樊仲藜，張 力，王慶棟，劉思婷,2，葉沅鑫

1. 中國(guó)測(cè)繪科學(xué)研究院，北京 100039; 2. 蘭州交通大學(xué)測(cè)繪與地理信息學(xué)院, 甘肅 蘭州 730070; 3. 西南交通大學(xué)地球科學(xué)與環(huán)境工程學(xué)院，四川 成都 611756

近年來(lái)，隨著SAR影像數(shù)據(jù)在地理國(guó)情監(jiān)測(cè)[1]、地表沉降變化檢測(cè)[2]等領(lǐng)域發(fā)揮的作用越來(lái)越顯著，其與光學(xué)影像截然不同的成像機(jī)理與特性引起學(xué)術(shù)界和工業(yè)界的廣泛關(guān)注。結(jié)合SAR影像與光學(xué)影像各自優(yōu)勢(shì)，實(shí)現(xiàn)對(duì)大量地表特性信息的互補(bǔ)，能夠更好地解決變化檢測(cè)、目標(biāo)識(shí)別、地物分類等任務(wù)。目前SAR影像和光學(xué)影像聯(lián)合應(yīng)用的關(guān)鍵前提在于影像匹配，影像匹配的本質(zhì)是在多幅影像之間尋找同名點(diǎn)的過(guò)程[3]，但是由于SAR影像獨(dú)特的成像機(jī)理，使其與光學(xué)影像之間存在顯著的非線性灰度差異，同時(shí)SAR影像還會(huì)受到嚴(yán)重的相干斑噪聲干擾[4]，這些影響因素使光學(xué)影像上的顯著特征點(diǎn)在SAR影像上極不穩(wěn)定，甚至失效，導(dǎo)致現(xiàn)有的影像匹配技術(shù)難以應(yīng)用于SAR影像和光學(xué)影像間的匹配，因此SAR影像和光學(xué)影像的匹配技術(shù)是目前遙感領(lǐng)域的難點(diǎn)，也是熱點(diǎn)研究問(wèn)題之一[5-7]。

影像匹配方法主要可分為基于特征的方法、基于區(qū)域的方法[8]和基于深度學(xué)習(xí)的方法[9-10]?；谔卣鞯姆椒ㄍㄟ^(guò)在輸入的影像之間提取共有特征，并檢測(cè)特征間的相似性來(lái)實(shí)現(xiàn)影像匹配。這些影像特征可以是點(diǎn)特征[11]、線特征[12]、面特征[13]或是局部不變性特征(SIFT[14]、SURF[15]等)。為了得到較好的匹配效果，檢測(cè)到的影像特征應(yīng)具有良好的穩(wěn)定性、可區(qū)分性和可重復(fù)性，但是在具有顯著非線性灰度差異的影像之間很難提取到具有高重復(fù)率的共有特征，從而導(dǎo)致匹配不到足夠數(shù)量的同名點(diǎn)，使得此類方法在SAR影像和光學(xué)影像匹配任務(wù)上失敗。

基于區(qū)域的方法也被稱為模板匹配方法，通過(guò)在影像間設(shè)置一定尺寸大小的影像窗口，并選擇某種相似性測(cè)度作為檢測(cè)影像相似性的準(zhǔn)則來(lái)實(shí)現(xiàn)影像匹配[16]。傳統(tǒng)的模板匹配方法通常是利用影像的灰度信息進(jìn)行匹配，代表性的方法有：互信息(mutual information,MI)[17]、歸一化互相關(guān)(normalized cross correlation,NCC)[18]等。由于NCC方法能夠?qū)τ跋耖g的線性灰度變化保持穩(wěn)健，因此被廣泛應(yīng)用于光學(xué)遙感影像匹配，但其對(duì)影像間的非線性灰度差異較為敏感。對(duì)于MI方法，其本質(zhì)是依賴于影像灰度的統(tǒng)計(jì)信息進(jìn)行匹配，與影像間灰度變化的方式相關(guān)性較弱，故理論上其對(duì)影像間的非線性灰度差異有一定的適應(yīng)性，然而MI方法忽略了鄰域像素的影響，使得計(jì)算結(jié)果容易陷入局部極值而導(dǎo)致誤匹配，并且其對(duì)模板窗口大小比較敏感，匹配過(guò)程計(jì)算量較大，因此MI方法在SAR影像和光學(xué)影像的匹配時(shí)效果不佳。

受到深度學(xué)習(xí)技術(shù)的影響，一些基于深度神經(jīng)網(wǎng)絡(luò)進(jìn)行SAR影像和光學(xué)影像匹配的方法被提出。文獻(xiàn)[9]提出了一種基于圖像風(fēng)格遷移的匹配方法，通過(guò)改進(jìn)的VGG-19網(wǎng)絡(luò)模型將SAR影像和光學(xué)影像進(jìn)行風(fēng)格互換，再利用SIFT或SURF等算法對(duì)互換風(fēng)格后的影像進(jìn)行匹配，該方法能夠獲取較多的正確匹配點(diǎn)。文獻(xiàn)[10]改進(jìn)了一種用于檢測(cè)圖像相似性的Siamese網(wǎng)絡(luò)模型[19-20]，通過(guò)抽取特征以及去除池化層對(duì)原網(wǎng)絡(luò)進(jìn)行優(yōu)化，該方法取得了較高的正確匹配率和匹配精度?，F(xiàn)階段基于深度學(xué)習(xí)的方法存在一定的局限性，首先需要花費(fèi)大量的時(shí)間和精力來(lái)制作訓(xùn)練數(shù)據(jù)集以及訓(xùn)練模型，并且由于地表形態(tài)構(gòu)造的多樣性，使得訓(xùn)練數(shù)據(jù)集難以具有良好的通用性，導(dǎo)致對(duì)于不同類型的地表形態(tài)，訓(xùn)練模型缺乏良好的泛化能力；其次，此類方法的執(zhí)行效率較低，由于需要將遙感影像輸入到神經(jīng)網(wǎng)絡(luò)模型中進(jìn)行訓(xùn)練，這個(gè)過(guò)程需要對(duì)海量樣本數(shù)據(jù)進(jìn)行迭代計(jì)算，對(duì)計(jì)算環(huán)境的軟硬件配置提出挑戰(zhàn)，并且非常耗時(shí)，因此這類方法難以達(dá)到實(shí)際應(yīng)用的要求，缺乏時(shí)效性。

除上述的各種方法之外，近年來(lái)，一些研究者發(fā)現(xiàn)影像的幾何結(jié)構(gòu)和形狀特征能夠有效抵抗影像間的非線性灰度差異[21-22]，因而將其引入到SAR影像和光學(xué)影像的匹配當(dāng)中，并取得了較好的匹配效果。這類方法不直接利用影像的灰度信息進(jìn)行匹配，而是首先利用如影像的梯度、局部自相似性(local self-similarity,LSS)、相位一致性(phase congruency,PC)等能有效地表達(dá)影像結(jié)構(gòu)的信息來(lái)構(gòu)建特征描述符，隨后建立描述符之間的相似性測(cè)度，并采用模板匹配的策略對(duì)影像進(jìn)行匹配。文獻(xiàn)[23]利用具有對(duì)比度和光照不變性的相位一致性特征對(duì)方向梯度直方圖(histograms of oriented gradients,HOG)[24]進(jìn)行了改進(jìn)，提出了相位一致性方向直方圖(histogram of orientated phase congruency,HOPC)，在光學(xué)和SAR影像的匹配任務(wù)上取得了較好的效果。但是，HOG和HOPC是在一個(gè)稀疏的采樣格網(wǎng)(非逐像素)內(nèi)進(jìn)行特征構(gòu)建，是一種相對(duì)稀疏的特征表達(dá)方式，難以精確地捕獲影像的細(xì)節(jié)結(jié)構(gòu)信息，另外其計(jì)算效率較低。文獻(xiàn)[25]利用方向梯度信息構(gòu)建了一種稠密結(jié)構(gòu)特征描述符(channel features of orientated gradients,CFOG),通過(guò)逐像素地提取鄰域結(jié)構(gòu)特征，增強(qiáng)了描述符對(duì)于結(jié)構(gòu)信息的細(xì)節(jié)表達(dá)能力，同時(shí)在頻率域中完成影像匹配，顯著提高了匹配性能及計(jì)算效率。不過(guò)，CFOG是利用水平和垂直方向的梯度內(nèi)插出的方向梯度，其計(jì)算精度還可進(jìn)一步提高。針對(duì)這一問(wèn)題，本文提出一種更精確的方向梯度計(jì)算方法，并以此構(gòu)建特征描述符進(jìn)行SAR影像和光學(xué)影像的匹配。

本文提出的SOFM方法使用影像的梯度信息構(gòu)建出一種表達(dá)影像結(jié)構(gòu)的特征描述符，即角度加權(quán)方向梯度AWOG,隨后選擇描述符之間的差值的平方和(SSD)作為用于匹配的相似性測(cè)度，并利用快速傅里葉變換將匹配過(guò)程轉(zhuǎn)換到頻率域進(jìn)行。設(shè)計(jì)了一套基于本文方法的自動(dòng)匹配流程，并通過(guò)設(shè)計(jì)多組試驗(yàn)來(lái)驗(yàn)證本文方法的有效性和匹配性能。

1 結(jié)構(gòu)特性表達(dá)的影像相似性

本節(jié)將給出利用影像梯度信息構(gòu)建AWOG描述符的詳細(xì)步驟，以及給出在頻率域中表達(dá)的影像匹配函數(shù)的詳細(xì)推導(dǎo)過(guò)程。

1.1 特征描述符

AWOG描述符的提出受到HOG[23]和CFOG[25]的啟發(fā)。HOG描述符最初被用于處理行人檢測(cè)任務(wù)，由于對(duì)影像結(jié)構(gòu)具有良好的表達(dá)能力，其應(yīng)用范圍被拓展到目標(biāo)識(shí)別[26]、影像分類[27]和影像檢索[28]等領(lǐng)域。CFOG是在HOG的基礎(chǔ)上發(fā)展的一種逐像素的結(jié)構(gòu)特征描述符，其使用方向梯度信息來(lái)表達(dá)影像結(jié)構(gòu)特性，能夠有效抵抗影像間的光照和對(duì)比度變化，以及非線性灰度差異。因此，本文同樣選擇利用影像梯度幅值和方向信息來(lái)構(gòu)建AWOG特征描述符。

如圖1所示，對(duì)于一幅給定的影像，首先使用x方向的濾波器dx=[-1,0,1]和y方向的濾波器dy=[-1,0,1]T對(duì)影像進(jìn)行濾波，得到影像在x軸和y軸方向的方向梯度gx和gy。

圖1 影像梯度方向OG和梯度幅值IG的計(jì)算Fig.1 Computational process of orientation and magnitude of image gradient

隨后根據(jù)式(1)和式(2)計(jì)算影像的梯度幅值IG和梯度方向OG。由于不同類型的傳感器所獲取的影像之間存在輻射強(qiáng)度反轉(zhuǎn)的現(xiàn)象[29]，因此為削弱SAR影像和光學(xué)影像間因輻射強(qiáng)度反轉(zhuǎn)所造成的影響，以提高描述符的穩(wěn)健性，對(duì)梯度方向OG進(jìn)行一步歸化操作[25,30]，根據(jù)式(2)計(jì)算得到的梯度方向OG其值域?yàn)閇0°,360°) ，將屬于區(qū)間[180°,360°)的值統(tǒng)一減去180°，OG的值域被歸化到[0°,180°)

(1)

(2)

如圖2所示，按角度間隔為22.5°將OG的值域8等分，等分的區(qū)間被9個(gè)角度方向所隔開(kāi)，由于特征描述符將在這9個(gè)角度方向上統(tǒng)計(jì)得到，故將它們稱為9個(gè)特征方向，按角度從小到大為9個(gè)特征方向賦予索引號(hào)0～8。

圖2 右特征方向索引表IRF的構(gòu)建Fig.2 Construction of index table IRF

對(duì)于影像中任意像元p(x,y)，其梯度方向值OG(x,y)存在于某兩個(gè)特征方向所構(gòu)成的角度區(qū)間內(nèi)(極少數(shù)恰好位于某特征方向上，但不影響后續(xù)計(jì)算)，將該角度區(qū)間中角度值較大的特征方向稱為像元的左特征方向LF，角度值較小的特征方向稱為像元的右特征方向RF，根據(jù)式(3)進(jìn)行計(jì)算

(3)

可得到像元右特征方向的索引號(hào)IRF(x,y)，由此對(duì)整幅影像進(jìn)行計(jì)算，得到影像的右特征方向索引表IRF，該索引表將用于加速后續(xù)的計(jì)算過(guò)程。

如圖3所示，對(duì)于任意像元p(x,y)，其梯度方向值OG(x,y)存在于其左特征方向LF和右特征方向RF所構(gòu)成的角度區(qū)間內(nèi)，因此可根據(jù)梯度方向值OG(x,y)與左右特征方向的角度差，將該像元的梯度幅值IG(x,y)加權(quán)分配到其左右特征方向上，OG(x,y)與其左右特征方向中某一個(gè)的角度差越小，說(shuō)明OG(x,y)距離該特征方向越近，從而IG(x,y)向該特征方向分配的值越大；反之，OG(x,y)與該特征方向的角度差越大，說(shuō)明OG(x,y)距離該特征方向越遠(yuǎn)，IG(x,y)向該特征方向分配的值越小。

圖3 加權(quán)梯度幅值WIG1、WIG2的計(jì)算過(guò)程Fig.3 Computational process of weighted magnitude of image gradient

基于上述加權(quán)原則，根據(jù)式(4)計(jì)算像素的梯度方向值OG(x,y)與右特征方向RF的角度差θ1(x,y)以及與左特征方向LF的角度差θ2(x,y)

(4)

隨后根據(jù)θ1(x,y)和θ2(x,y)對(duì)梯度幅值IG(x,y)進(jìn)行加權(quán)分配，由此對(duì)整幅影像進(jìn)行處理，得到影像的加權(quán)梯度幅值WIG1和WIG2。

逐像元為中心在WIG1和WIG2上同時(shí)開(kāi)辟尺寸為3×3像素的統(tǒng)計(jì)窗口，對(duì)兩個(gè)窗口內(nèi)9個(gè)特征方向上的加權(quán)梯度幅值進(jìn)行統(tǒng)計(jì)，利用右特征方向索引表IRF加速統(tǒng)計(jì)過(guò)程，聯(lián)合兩個(gè)統(tǒng)計(jì)窗口，將統(tǒng)計(jì)結(jié)果以向量的形式輸出，該向量即為統(tǒng)計(jì)窗口中心像元的特征向量。使用索引表進(jìn)行統(tǒng)計(jì)相比于不使用索引表的統(tǒng)計(jì)方式，在統(tǒng)計(jì)窗口內(nèi)的每個(gè)像元處均可減少兩次乘法、兩次除法和一次減法運(yùn)算而只保留兩次加法運(yùn)算，因此，使用索引表進(jìn)行統(tǒng)計(jì)將極大地提高特征描述符的提取速度。

將每個(gè)像元處的特征向量沿垂直于影像平面的z軸方向排列便形成初步的影像描述符。隨后使用z軸方向的濾波器dz=[1,3,1]T對(duì)初步描述符進(jìn)行濾波，z軸上的濾波操作對(duì)不同特征方向的特征值起到平滑作用，該操作能夠減小影像間局部幾何和強(qiáng)度變形引起的畸變所造成的影響[25]。對(duì)z軸濾波后的描述符進(jìn)行歸一化便得到最終的AWOG特征描述符。

綜上，AWOG特征描述符的構(gòu)建過(guò)程如圖4所示。

圖4 AWOG特征描述符的構(gòu)建Fig.4 Construction of AWOG

1.2 相似性測(cè)度

由于AWOG描述符具有多層級(jí)結(jié)構(gòu)，相比于傳統(tǒng)的單層描述符來(lái)說(shuō)其數(shù)據(jù)量較大，在空間域中進(jìn)行匹配較為耗時(shí)，因此選擇在頻率域中進(jìn)行匹配[31]，在頻率域中能夠顯著提高匹配時(shí)的計(jì)算速度[25]，選擇描述符之間差值的平方和作為同名點(diǎn)匹配時(shí)的相似性準(zhǔn)則，并利用快速傅里葉變換推導(dǎo)出在頻率域中表達(dá)的影像匹配函數(shù)。

差值的平方與表達(dá)的AWOG描述符之間的相似性測(cè)度表示為

d(v)=∑x[A1(x)-A2(x-v)]2

(5)

式中，x表示某像元的像素坐標(biāo)；A1(x)和A2(x-v)分別表示參考影像描述符和輸入影像描述符被模板窗口覆蓋的區(qū)塊；變量v表示輸入影像描述符的搜索窗口中某像元相對(duì)于搜索窗口中心像元的偏移向量；d(v)為兩個(gè)描述符區(qū)塊之間的差值的平方和函數(shù)。

為了得到參考影像中的某像元在輸入影像上的最佳匹配位置，需對(duì)式(5)求最小值，由此，同名點(diǎn)的搜索問(wèn)題轉(zhuǎn)化為計(jì)算使式(5)達(dá)到最小值的偏移向量vi。進(jìn)而給出如下影像匹配函數(shù)

(6)

對(duì)于模板匹配的過(guò)程，需要在輸入影像搜索窗口中的每個(gè)像元位置上，計(jì)算與模板窗口間的相似性值，該計(jì)算過(guò)程非常耗時(shí)，因此，為了提高計(jì)算效率，采用快速傅里葉變換將該過(guò)程轉(zhuǎn)換到頻率域進(jìn)行加速。

將式(6)中的平方項(xiàng)展開(kāi)

2∑xA1(x)·A2(x-v)}

(7)

對(duì)于式(7)，等式右側(cè)第1項(xiàng)的計(jì)算值是一個(gè)常數(shù)，因此計(jì)算時(shí)可將其忽略。后兩項(xiàng)可以視為描述符之間的相關(guān)計(jì)算，由于空間域中的相關(guān)對(duì)應(yīng)頻率域中的相乘，因此可利用快速傅里葉變換(二維或三維)將后兩項(xiàng)轉(zhuǎn)換到頻率域中計(jì)算。影像的匹配函數(shù)轉(zhuǎn)化為

2F-1[F(A2(x-v))F*(A1(x))]}

(8)

式中，F(xiàn)和F-1表示快速傅里葉變換的正變換和逆變換；F*表示對(duì)F求復(fù)共軛。

對(duì)于式(8),等式右側(cè)第1項(xiàng)表達(dá)輸入影像描述符與自身的相關(guān)計(jì)算，對(duì)于進(jìn)行了歸一化的特征描述符如AWOG、FLSS[32]、FHOG[32]等，右側(cè)第1項(xiàng)的計(jì)算值幾乎為一個(gè)常數(shù)，同時(shí)其在影像不同位置處的計(jì)算值的變化率對(duì)于式(8)右側(cè)第2項(xiàng)的計(jì)算結(jié)果幾乎沒(méi)有影響，因此參照文獻(xiàn)[32]，選擇忽略式(8)右側(cè)第1項(xiàng)而只保留第2項(xiàng)。因此同名點(diǎn)匹配被轉(zhuǎn)化為計(jì)算對(duì)式(8)右側(cè)第2項(xiàng)反符號(hào)求得最大值時(shí)對(duì)應(yīng)的偏移向量vi。最終得到如下影像匹配函數(shù)

(9)

式(9)表達(dá)了兩幅影像的描述符在頻率域中的互相關(guān)，該計(jì)算會(huì)得到一幅相關(guān)圖，因此通過(guò)確定最大反應(yīng)值所在位置，便可以得到同名匹配點(diǎn)。

為驗(yàn)證上述影像匹配函數(shù)的有效性，利用x方向的相似性曲線將其與MI[17]和NCC[18]進(jìn)行比較，選擇尺寸為41×41像素的模板窗口在一組高分辨率城區(qū)影像上進(jìn)行測(cè)試，測(cè)試影像如圖5所示。

圖5 相似性曲線測(cè)試影像Fig.5 Images for testing similarity curve

圖6展示SOFM、MI[17]和NCC[18]等3種方法在X方向上的相似性曲線，通過(guò)對(duì)比可以發(fā)現(xiàn)，SOFM方法的相似性曲線比較光滑并成功檢測(cè)到了正確匹配位置，相比之下，MI和NCC方法的相似性曲線則不夠光滑，且峰值位置與正確匹配位置存在一定偏差。通過(guò)相似性曲線的對(duì)比，初步證明了本文方法的有效性，關(guān)于本文方法匹配性能的詳細(xì)分析見(jiàn)第3節(jié)。

2 SAR影像和光學(xué)影像模板匹配流程

(1) 將SAR影像作為參考影像，光學(xué)影像作為輸入影像，為兩影像構(gòu)建影像金字塔。若影像間分辨率差異不超過(guò)2倍，一般為影像金字塔設(shè)置4個(gè)層級(jí)，層級(jí)間的縮放比為2；若影像間分辨率差異較大，則根據(jù)分辨率差異倍數(shù)為分辨率較高的影像增加一定的金字塔影像層數(shù)，并通過(guò)建立查找表的方式在分辨率差異較小的層級(jí)間進(jìn)行匹配。由低層級(jí)影像構(gòu)建高層級(jí)影像時(shí)，先對(duì)低層級(jí)影像進(jìn)行高斯濾波，后按照雙線性內(nèi)插的采樣方式生成高層級(jí)影像。隨后由頂層金字塔影像進(jìn)入后續(xù)處理步驟。

圖6 SOFM、MI和NCC方法的相似性曲線Fig.6 Similarity curves of SOFM, MI and NCC

(2) 為保證所提取影像特征點(diǎn)的均勻性，選擇使用分塊Harris算子，計(jì)算參考影像中每個(gè)像元處的Harris特征值，隨后將影像劃分為n×n個(gè)互不重疊的格網(wǎng)區(qū)域，提取每個(gè)格網(wǎng)內(nèi)特征值較大的k個(gè)點(diǎn)作為特征點(diǎn)，整幅影像將得到k×n×n個(gè)特征點(diǎn)。關(guān)于劃分的格網(wǎng)數(shù)n和格網(wǎng)內(nèi)特征點(diǎn)數(shù)k的選擇，一般經(jīng)驗(yàn)性的將k設(shè)置為1～5個(gè)點(diǎn)，再根據(jù)期望得到的總點(diǎn)數(shù)計(jì)算出格網(wǎng)數(shù)n。

(3) 考慮到影像之間存在旋轉(zhuǎn)與尺度差異，而AWOG描述符不具備旋轉(zhuǎn)與尺度不變性，因此在匹配之前對(duì)影像進(jìn)行image-reshaping操作(圖7)來(lái)消除影像間的旋轉(zhuǎn)和尺度差異，目前所允許的影像間的縮放倍數(shù)約為0.2～5。處理步驟為：①以特征點(diǎn)為中心在參考影像上確定一定尺寸的搜索區(qū)域；②給定搜索區(qū)域的近似高程面(區(qū)域平均高程面或已有概略DEM)，利用參考影像的定向參數(shù)將搜索區(qū)域投影到近似高程面上，得到搜索區(qū)域?qū)?yīng)的物方面元；③利用影像成像模型(如RFM模型)將物方面元反投影到輸入影像上得到相應(yīng)的搜索區(qū)域范圍；④建立參考影像的搜索區(qū)域和輸入影像的搜索區(qū)域之間的透視變換或仿射變換模型，并根據(jù)所建立的變換模型對(duì)輸入影像的搜索區(qū)域進(jìn)行重采樣，基本上消除其與參考影像之間的旋轉(zhuǎn)和尺度差異。

(4) 根據(jù)本文提出的SOFM算法，對(duì)經(jīng)過(guò)步驟(3)處理的參考影像和輸入影像的影像塊提取AWOG描述符，并根據(jù)給出的影像匹配函數(shù)完成同名點(diǎn)匹配。

(5) 重復(fù)步驟(3)—步驟(4)，直至遍歷所有特征點(diǎn)，得到初始的同名點(diǎn)集合，由于匹配過(guò)程在提取的影像塊上進(jìn)行，此時(shí)同名點(diǎn)坐標(biāo)由局部像素坐標(biāo)系所表達(dá)，因此通過(guò)坐標(biāo)轉(zhuǎn)換將同名點(diǎn)坐標(biāo)表達(dá)在原始輸入影像的像素坐標(biāo)系上，并輸出初始同名點(diǎn)集合。

圖7 影像匹配時(shí)的image-reshaping操作 Fig.7 Image-reshaping processing used for image matching

(6) 構(gòu)建二維透視變換表達(dá)的RANSAC框架，設(shè)置一定的粗差閾值及迭代次數(shù)，將步驟(5)輸出的初始同名點(diǎn)集合輸入該框架，迭代剔除殘差大于閾值的點(diǎn)對(duì)，輸出內(nèi)點(diǎn)率最高的同名點(diǎn)集合。

(7) 構(gòu)建二維透視變換表達(dá)的整體最小二乘模型，將步驟(6)剔除粗差后的同名點(diǎn)集合作為輸入，計(jì)算單個(gè)點(diǎn)的殘差和整體的標(biāo)準(zhǔn)差，迭代刪除殘差較大的同名點(diǎn)對(duì)，直至標(biāo)準(zhǔn)差小于給定的閾值，將剩余的同名點(diǎn)集合作為匹配結(jié)果，并由此計(jì)算當(dāng)前層級(jí)金字塔影像的仿射變換參數(shù)。

(8) 重復(fù)步驟(2)—步驟(7)，并將計(jì)算出的仿射變換參數(shù)向低層級(jí)影像傳遞，從而實(shí)現(xiàn)由粗到精的影像匹配策略，直至完成金字塔底層影像的匹配，并記錄最終輸出結(jié)果。

3 試驗(yàn)與分析

為驗(yàn)證所提出方法的有效性，本文試驗(yàn)將僅考慮影像間的平移差異，對(duì)于影像間的旋轉(zhuǎn)和尺度差異可通過(guò)第2節(jié)中的image-reshaping操作進(jìn)行消除，使得影像間只剩余平移參數(shù)未知，而平移參數(shù)與影像自身的無(wú)控定位精度有關(guān)。因此，為簡(jiǎn)化試驗(yàn)步驟，試驗(yàn)中所使用的參考影像均為經(jīng)過(guò)多視處理后的SAR影像的數(shù)字正射影像圖(以下簡(jiǎn)稱DOM),輸入影像均為經(jīng)過(guò)正射糾正后的光學(xué)影像的DOM。試驗(yàn)過(guò)程中使用的計(jì)算機(jī)配置為Intel(R)Core(TM)i9-10885H CPU @2.4 GHz,64 GB。

3.1 試驗(yàn)數(shù)據(jù)

選擇5組不同分辨率且覆蓋不同地面場(chǎng)景的SAR影像和光學(xué)影像的DOM作為試驗(yàn)數(shù)據(jù)，如圖8所示。具體的數(shù)據(jù)描述見(jiàn)表1。

表1 試驗(yàn)數(shù)據(jù)描述

圖8為試驗(yàn)影像的縮略圖，可以看出每組影像間均存在顯著的非線性灰度差異。

3.2 SOFM方法對(duì)模板窗口大小的適應(yīng)性分析

使用表1中的第5組試驗(yàn)數(shù)據(jù)，按照第2節(jié)的匹配流程對(duì)影像進(jìn)行匹配，固定搜索窗口為31×31像素，改變模板窗口的尺寸(41×41像素～101×101像素),記錄最終輸出的點(diǎn)數(shù)以及匹配到的總點(diǎn)數(shù)，試驗(yàn)結(jié)果見(jiàn)表2。

表2 不同模板窗口下的影像匹配結(jié)果

分析表2中的結(jié)果可知，隨著設(shè)置的模板窗口尺寸的增大，最終輸出的點(diǎn)數(shù)逐漸增多且趨于穩(wěn)定，輸出點(diǎn)數(shù)占匹配總點(diǎn)數(shù)的比率逐漸上升，當(dāng)模板窗口尺寸設(shè)置為61×61像素時(shí)，輸出點(diǎn)的數(shù)量以其占總點(diǎn)數(shù)的比率均已達(dá)到較高水平，因此在應(yīng)用中將模板窗口的尺寸固定為61×61像素即可。

3.3 SOFM方法對(duì)搜索窗口大小的適應(yīng)性分析

使用表1中的第5組試驗(yàn)數(shù)據(jù)，對(duì)試驗(yàn)影像建立影像金字塔(4級(jí))，相鄰層級(jí)間的縮放倍數(shù)為2，選擇61×61像素的模板窗口，對(duì)頂層金字塔影像進(jìn)行整個(gè)影像范圍內(nèi)的搜索，計(jì)算影像間的仿射變換參數(shù)并向低層級(jí)影像傳遞，并對(duì)低層級(jí)影像重新提取特征點(diǎn)，在低層級(jí)影像中設(shè)置21×21像素的搜索窗口，記錄最終輸出的點(diǎn)數(shù)以及匹配到的總點(diǎn)數(shù)，試驗(yàn)結(jié)果見(jiàn)表3。

表3 影像金字塔匹配結(jié)果

圖8 試驗(yàn)影像Fig.8 The images used in the experiments

對(duì)表3進(jìn)行分析可知，在影像匹配時(shí)，通過(guò)使用對(duì)影像建立影像金字塔的方式，既能加快影像匹配的速度，又能穩(wěn)定地獲取數(shù)量足夠多的同名點(diǎn)。

3.4 匹配性能對(duì)比分析

為驗(yàn)證SOFM方法的匹配性能，使用表1中的第1～4組試驗(yàn)數(shù)據(jù)，選擇CFOG[25]、HOPC[23]、DLSC[33]、MI[17]、NCC[18]等方法進(jìn)行對(duì)比分析與評(píng)價(jià)，對(duì)比的內(nèi)容包括：正確匹配率(正確匹配率=正確匹配點(diǎn)數(shù)/總匹配點(diǎn)數(shù))、計(jì)算效率和均方根誤差(RMSE)等3個(gè)方面。試驗(yàn)步驟為：①在參考影像上提取200個(gè)Harris特征點(diǎn)；②利用影像金字塔為每個(gè)特征點(diǎn)在輸入影像上確定一個(gè)21×21像素的搜索窗口，采用不同尺寸的模板窗口(31×31像素～91×91像素)對(duì)影像進(jìn)行匹配并記錄不同方法的性能表現(xiàn)。

圖9展示了參與對(duì)比的6種匹配方法在4組試驗(yàn)影像上的正確匹配率(容差為1.5像素)。

通過(guò)對(duì)圖9進(jìn)行分析可以得到以下結(jié)論。

(1) NCC方法總體上在各組試驗(yàn)中的正確匹配率最低，甚至在第3組試驗(yàn)中出現(xiàn)了近乎匹配失敗的情況，究其原因，NCC方法直接利用影像的灰度信息進(jìn)行匹配，其能夠?qū)τ跋耖g的線性灰度變化具有良好的適應(yīng)性，但對(duì)影像間的非線性灰度差異非常敏感，因此NCC方法在進(jìn)行光學(xué)影像間的匹配時(shí)表現(xiàn)良好，但在匹配SAR影像和光學(xué)影像時(shí)效果較差。

(2) MI方法在多組試驗(yàn)中的表現(xiàn)雖然優(yōu)于NCC方法，但其總體上的正確匹配率仍然較低，MI方法本質(zhì)是依賴影像灰度的統(tǒng)計(jì)信息進(jìn)行匹配，故理論上來(lái)講其匹配效果應(yīng)與影像間灰度變化的方式相關(guān)性較弱，但MI方法在匹配過(guò)程中不考慮鄰域像素的影響，導(dǎo)致計(jì)算結(jié)果容易陷入局部極值而產(chǎn)生誤匹配，因此MI方法不適用于SAR影像和光學(xué)影像的匹配。

圖9 SOFM、CFOG、HOPC、DLSC、MI和NCC在不同模板尺寸下的正確匹配率Fig.9 The correct matching ratios of SOFM, CFOG, HOPC, DLSC, MI and NCC with different template sizes

(3) HOPC和DLSC方法的表現(xiàn)較好，總體來(lái)看，HOPC在各組試驗(yàn)影像上的正確匹配率均高于DLSC。這兩種方法屬于基于影像結(jié)構(gòu)特性的匹配方法，HOPC和DLSC方法分別利用了影像的相位一致性(PC)和局部自相似性(LSS)信息來(lái)表達(dá)影像的幾何構(gòu)造，但是由圖9可知，相比于SOFM這兩種方法的匹配性能受模板尺寸變化的影響較大，在模板尺寸較大時(shí)，能夠得到較高的正確匹配率；但當(dāng)設(shè)置的模板尺寸逐漸減小時(shí)，這兩種方法的正確匹配率會(huì)出現(xiàn)較明顯的降低。

(4) 相比于上述的幾種方法，SOFM和COFG在各試驗(yàn)像對(duì)上的正確匹配率更高，但是整體而言，SOFM方法在每組試驗(yàn)中的表現(xiàn)均優(yōu)于CFOG。SOFM的正確匹配率受模板尺寸變化的影響更小，性能更加穩(wěn)定，在模板尺寸較小時(shí)仍能保證較高的正確匹配率。當(dāng)設(shè)置模板尺寸為31×31像素時(shí)，SOFM方法在4組試驗(yàn)上的正確匹配率平均要比CFOG方法高17.8%。分析其原因，可能是由于CFOG方法是利用水平和垂直方向的梯度插值得到多個(gè)固定方向的方向梯度，并對(duì)方向梯度進(jìn)行高斯卷積來(lái)得到目標(biāo)像素的特征向量。而這種計(jì)算方式會(huì)引入一定程度的模糊，降低特征向量之間的區(qū)分性，因此當(dāng)設(shè)置模板尺寸較小或影像質(zhì)量較差時(shí)，會(huì)導(dǎo)致其匹配效果出現(xiàn)一定程度的下降。而SOFM方法則是以方向間的角度差為權(quán)重，將像素的梯度值加權(quán)分配到與梯度方向相關(guān)性最強(qiáng)的兩個(gè)固定方向上，同時(shí)結(jié)合索引表進(jìn)行鄰域統(tǒng)計(jì)得到目標(biāo)像素的特征向量，增強(qiáng)了特征向量的獨(dú)特性，因此其穩(wěn)健性更好，匹配效果更優(yōu)。

圖10展示了各方法在4組試驗(yàn)中的平均運(yùn)行時(shí)間。

圖10 SOFM、CFOG、HOPC、DLSC、MI和NCC在不同模板尺寸下的平均運(yùn)行時(shí)間Fig.10 The average running times of SOFM, CFOG, HOPC, DLSC, MI and NCC with different template sizes

分析圖10可得如下結(jié)論。

(1) 總體上MI方法的計(jì)算效率最低，其次較低的是DLSC方法，NCC方法的計(jì)算效率雖然較高，但是其匹配效果較差，HOPC方法的計(jì)算效率略高于NCC。CFOG方法的計(jì)算效率最高，SOFM方法的計(jì)算效率略低于CFOG方法。

(2) SOFM和CFOG方法的時(shí)間消耗隨著模板尺寸的增大只有微小增加，而其他4種方法的時(shí)間消耗會(huì)隨著模板尺寸的增大而迅速增加。當(dāng)模板尺寸較大時(shí)(47×47像素及以上)，HOPC和DLSC方法能夠得到較高的正確匹配率，但同時(shí)也會(huì)有更長(zhǎng)的時(shí)間消耗。

(3) 當(dāng)設(shè)置模板尺寸91×91像素時(shí)，SOFM、CFOG、HOPC、DLSC、MI、NCC等方法的時(shí)間消耗分別為2.51、1.31、11.18、30.17、29.44和14.99 s，通過(guò)數(shù)據(jù)對(duì)比可以看出，除CFOG比SOFM方法快1.2 s之外，NCC方法的時(shí)間消耗約是SOFM方法的6倍，最慢的DLSC方法的時(shí)間消耗約是SOFM方法的12倍。因此，相較之下，SOFM方法在保證較高的正確匹配率的前提下，達(dá)到了快速匹配的效果。

圖11展示了當(dāng)模板尺寸為91×91像素時(shí)各方法在4組數(shù)據(jù)上的RMSE。由圖11可知SOFM方法在各組試驗(yàn)影像上均具有最小的RMSE，這說(shuō)明相較于其余幾種方法，SOFM方法具有最高的匹配精度，能夠在進(jìn)行SAR影像和光學(xué)影像的匹配時(shí)更具優(yōu)勢(shì)。

圖11 SOFM、CFOG、DLSC、HOPC、MI和NCC在各試驗(yàn)影像上的RMSEFig.11 The RMSEs of correct matches of SOFM, CFOG, HOPC, DLSC, MI and NCC

圖12展示了SOFM方法在各組試驗(yàn)影像上的匹配結(jié)果(模板尺寸61×61像素)，可以看出匹配到的同名點(diǎn)數(shù)量充足且分布均勻，觀察局部放大圖能夠看到同名點(diǎn)的對(duì)應(yīng)位置關(guān)系是正確的。

通過(guò)以上的各組試驗(yàn)和對(duì)比分析，證明了本文提出的SOFM方法在正確匹配率、匹配穩(wěn)定性、計(jì)算效率和匹配精度等方面均更具優(yōu)勢(shì)，能夠快速、高效地完成SAR影像與光學(xué)影像之間的自動(dòng)匹配任務(wù)。

4 討 論

針對(duì)當(dāng)前遙感領(lǐng)域較為困難的SAR影像和光學(xué)影像間的匹配問(wèn)題，本文提出了一種基于影像結(jié)構(gòu)特性的快速匹配方法SOFM。該方法利用影像的梯度信息構(gòu)建了一種特征描述符—角度加權(quán)方向梯度AWOG，然后利用描述符之間的差值的平方和建立了用于模板匹配的相似性測(cè)度，并推導(dǎo)出了利用快速傅里葉變換在頻率域中表達(dá)的影像匹配函數(shù)，在頻率域中進(jìn)行影像匹配能夠大幅提高計(jì)算效率，隨后給出了使用本文方法的模板匹配流程。通過(guò)與現(xiàn)有的影像匹配方法CFOG、HOPC、DLSC、MI和NCC在多個(gè)方面進(jìn)行對(duì)比分析，結(jié)果顯示：SOFM方法能夠較好地抵抗SAR影像與光學(xué)影像之間的非線性灰度差異，其匹配性能良好且更具穩(wěn)定性，同時(shí)在計(jì)算效率和匹配精度方面也有著優(yōu)異的表現(xiàn)。SOFM方法能夠快速地完成SAR影像與光學(xué)影像之間的自動(dòng)匹配任務(wù)，獲得精度高且分布均勻的同名匹配點(diǎn)，為后續(xù)的相關(guān)應(yīng)用提供了可靠保障。

圖12 SOFM的匹配結(jié)果Fig.12 The image matching results of SOFM

考慮到SAR影像生成的正射影像圖(DOM)的平面精度基本只與傳感器的定軌精度有關(guān)，與瞬時(shí)成像姿態(tài)角幾乎沒(méi)有相關(guān)性，而目前衛(wèi)星定軌精度已經(jīng)達(dá)到足夠高的水平，因此SAR影像的DOM能夠保證較高的平面精度[34]。筆者下一步的工作將使用由天繪二號(hào)衛(wèi)星數(shù)據(jù)生成的正射影像圖(DOM)及數(shù)字高程模型(DEM)作為控制信息，將光學(xué)原始影像與之進(jìn)行匹配，匹配得到的同名點(diǎn)將作為控制點(diǎn)來(lái)精化光學(xué)原始影像的RPC參數(shù)，從而可形成一套完全基于國(guó)產(chǎn)SAR衛(wèi)星數(shù)據(jù)和光學(xué)衛(wèi)星數(shù)據(jù)的無(wú)控測(cè)圖方法流程，以期在無(wú)法實(shí)地測(cè)得控制點(diǎn)的境外地區(qū)實(shí)現(xiàn)高精度的無(wú)控測(cè)圖。

試驗(yàn)過(guò)程中也發(fā)現(xiàn)本文方法目前存在一定不足：①SOFM方法中的AWOG描述符不具有旋轉(zhuǎn)和尺度不變性，但可利用影像的定向參數(shù)基本消除影像間的旋轉(zhuǎn)和尺度差異(見(jiàn)第2節(jié))。試驗(yàn)結(jié)果表明，對(duì)于地形起伏較小的區(qū)域如平地、丘陵地甚至是山體高度變化不大的山區(qū)，本文方法都能得到較好的匹配結(jié)果，但是對(duì)于高程變化劇烈的山區(qū)如我國(guó)橫斷山區(qū)，方法便有可能失效，因此尚需要進(jìn)一步研究能夠抵抗局部影像間旋轉(zhuǎn)和尺度差異的方法，以更好地解決復(fù)雜地形區(qū)域的SAR影像與光學(xué)影像之間的匹配。②SOFM方法依賴于影像中蘊(yùn)含的結(jié)構(gòu)特征進(jìn)行匹配，如果影像的結(jié)構(gòu)特征不夠豐富，如影像中含有大面積的森林區(qū)域，則方法得到的有效特征點(diǎn)數(shù)量會(huì)有一定程度的減少，同時(shí)其分布情況也會(huì)在一定程度上變差，因此，在進(jìn)一步的研究中將考慮引入紋理特征來(lái)輔助進(jìn)行匹配。