基于地形輔助的無(wú)人機(jī)載InSAR圖像分區(qū)配準(zhǔn)方法

謝 鑫 鄧云開*② 楊志軍 田衛(wèi)明④

①(北京理工大學(xué)信息與電子學(xué)院雷達(dá)技術(shù)研究所 北京 100081)

②(嵌入式實(shí)時(shí)信息處理技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室 北京 100081)

③(北京理工大學(xué)重慶創(chuàng)新中心 重慶 401120)

④(北京理工大學(xué)前沿技術(shù)研究院 濟(jì)南 250300)

1 引言

干涉合成孔徑雷達(dá)(Interferometric Synthetic Aperture Radar,InSAR)廣泛應(yīng)用于滑坡監(jiān)測(cè)和地形測(cè)繪等領(lǐng)域[1]。InSAR圖像配準(zhǔn)是指根據(jù)兩幅圖像同名點(diǎn)之間的坐標(biāo)映射關(guān)系，將輔圖像重采樣到與主圖像相同的像素網(wǎng)格，使配準(zhǔn)后兩幅圖像的同名點(diǎn)對(duì)位于同一分辨單元內(nèi)[2]。InSAR產(chǎn)生的復(fù)圖像必須經(jīng)過(guò)精確配準(zhǔn)才能保證圖像對(duì)之間具有良好的相干性。因此，配準(zhǔn)作為干涉測(cè)量中重要的一環(huán)，直接影響到干涉圖的質(zhì)量和地形測(cè)量的精度[3,4]。

針對(duì)InSAR圖像的配準(zhǔn)，其核心是偏移量的提取與估計(jì)。偏移量的提取與估計(jì)方法一般分為兩類，一類是利用較為精準(zhǔn)的外部軌道等信息，結(jié)合干涉測(cè)量幾何直接計(jì)算同名點(diǎn)對(duì)之間的偏移量[5]。這類方法通常以成像區(qū)域的平均高程作為整個(gè)區(qū)域的統(tǒng)一高程，即需要滿足平地假設(shè)[6]。相對(duì)于飛行高度，當(dāng)成像場(chǎng)景相對(duì)平坦，沒有大高程起伏時(shí)，配準(zhǔn)精度能達(dá)到亞像素級(jí)。因此，這類方法通常應(yīng)用于星載InSAR測(cè)量[7]。另一類方法則是基于In-SAR圖像中的信息來(lái)估計(jì)偏移量。這類方法主要采用最大互相關(guān)系數(shù)(Maximum Cross Correlation,MCC)[8]、最大譜函數(shù)(Maximum Spectrum Function,MSF)[9]或最小平均波動(dòng)函數(shù)(Average Fluctuation Function,AFF)[10]等作為測(cè)度函數(shù)，基于滑動(dòng)窗口搜索，實(shí)現(xiàn)同名點(diǎn)對(duì)之間的偏移量估計(jì)[11]。根據(jù)提取到的偏移量，基于所選像素點(diǎn)的坐標(biāo)位置構(gòu)建多項(xiàng)式擬合模型，通過(guò)最小二乘法估計(jì)模型參數(shù)，建立起從主圖像到輔圖像之間的映射關(guān)系[12]。在實(shí)際處理中，基于多項(xiàng)式擬合的方法需根據(jù)干涉圖像對(duì)之間偏移量具體情況選擇合適的多項(xiàng)式模型階數(shù)，但對(duì)于空變的偏移量易失效。

與傳統(tǒng)的InSAR配準(zhǔn)算法不同，其他利用SAR圖像信息的配準(zhǔn)方法不再使用測(cè)度函數(shù)來(lái)匹配同名點(diǎn)，而是基于SAR圖像中的點(diǎn)[13]、線[14]和邊緣結(jié)構(gòu)[15]等顯著幾何特征來(lái)估計(jì)配準(zhǔn)偏移量。然而，傳統(tǒng)的特征提取方法十分復(fù)雜，且難以獲得更全面和更有表現(xiàn)力的特征[16]。隨著深度學(xué)習(xí)的發(fā)展，卷積神經(jīng)網(wǎng)絡(luò)(Convolution Neural Network,CNN)逐步應(yīng)用于圖像配準(zhǔn)中的特征提取和匹配[17]。這些方法更適合于具有明顯變化的圖像處理場(chǎng)景，例如多角度圖像配準(zhǔn)[18]和異源圖像配準(zhǔn)[19]等，需要相對(duì)固定的成像場(chǎng)景[20]和大量的圖像數(shù)據(jù)集進(jìn)行訓(xùn)練，不適合InSAR圖像配準(zhǔn)。

與機(jī)載InSAR相比，無(wú)人機(jī)體積小，飛行易受氣流擾動(dòng)的影響，航跡非線性[21]；載重有限，通常僅負(fù)載單天線，需要多次飛行才能實(shí)現(xiàn)干涉測(cè)量，航跡不平行[22]。非線性、不平行的飛行軌跡導(dǎo)致兩幅SAR圖像存在幾何畸變，偏移量變化復(fù)雜。另外，無(wú)人機(jī)飛行高度有限，測(cè)量區(qū)域的高程與飛行高度相差不大，星載InSAR或常規(guī)機(jī)載InSAR測(cè)量中的平地近似假設(shè)不再成立。在復(fù)雜地形條件下，干涉圖像對(duì)的不同同名像素之間的偏移量具有明顯的空變特性，特別是在長(zhǎng)基線模式下，平地的區(qū)域偏移量較小，而在有高程區(qū)域，偏移量可以達(dá)到數(shù)個(gè)像素[23]。當(dāng)成像區(qū)域內(nèi)高程起伏較大時(shí)，偏移量空變明顯，傳統(tǒng)的全局多項(xiàng)式擬合配準(zhǔn)方法不再適用[24,25]。

本文圍繞無(wú)人機(jī)載InSAR圖像配準(zhǔn)問(wèn)題開展了研究，分析了無(wú)人機(jī)載InSAR成像投影幾何模型，建立了適用于無(wú)人機(jī)載InSAR的配準(zhǔn)偏移量模型，論證了地形對(duì)無(wú)人機(jī)載InSAR圖像配準(zhǔn)的影響。為了解決復(fù)雜地形條件下的配準(zhǔn)難題，本文提出了一種基于圖像分區(qū)的精配準(zhǔn)算法。首先，基于無(wú)人機(jī)載InSAR偏移量模型，結(jié)合航跡信息，通過(guò)設(shè)定偏移量門限得到高程門限；其次，利用外部地形數(shù)據(jù)進(jìn)行圖像分區(qū)，并進(jìn)行連通域分析，剔除面積較小的區(qū)域；最后，構(gòu)建帶邊界約束的子塊多項(xiàng)式全局?jǐn)M合模型，引入Lagrange算子，利用最小二乘法聯(lián)合求解模型參數(shù)，并進(jìn)行圖像重采樣獲取精配準(zhǔn)后的圖像。利用P波段無(wú)人機(jī)載InSAR實(shí)測(cè)數(shù)據(jù)驗(yàn)證了該方法的有效性，實(shí)驗(yàn)數(shù)據(jù)結(jié)果表明配準(zhǔn)后的相干性大大提高，殘差點(diǎn)數(shù)量減少。在高程較大的區(qū)域，干涉相位的質(zhì)量顯著提高。

2 無(wú)人機(jī)載SAR成像幾何模型

常規(guī)線性孔徑SAR工作在正側(cè)視條帶模式下，天線波束指向始終與平臺(tái)運(yùn)動(dòng)方向垂直[26]。圖1所示為常規(guī)線性孔徑SAR成像幾何示意圖。常規(guī)線性孔徑SAR成像時(shí)，是將觀測(cè)區(qū)域的三維地形投影在二維成像平面上。以平臺(tái)運(yùn)動(dòng)方向作為x方向，豎直向上方向作為z方向，根據(jù)右手準(zhǔn)則，構(gòu)建三維成像直角坐標(biāo)系。其中，x-y平面為成像平面，雷達(dá)運(yùn)動(dòng)軌跡的中心位置S為坐標(biāo)原點(diǎn)，軌跡兩端分別為S1和S2[27]。

圖1 常規(guī)線性孔徑SAR成像幾何Fig.1 Conventional linear aperture SAR imaging geometry

其中，H代表無(wú)人機(jī)平臺(tái)相對(duì)于成像平面的高度。

在實(shí)際的無(wú)人機(jī)載SAR系統(tǒng)中，平臺(tái)易受到氣流擾動(dòng)等外界因素的影響，難以實(shí)現(xiàn)勻速直線運(yùn)動(dòng)，呈現(xiàn)出非線性的航跡特性。圖2所示為實(shí)際無(wú)人機(jī)載SAR成像幾何。此時(shí)，天線相位中心(Antenna Phase Center,APC)與地形點(diǎn)P(xp,yp,zp)之間的距離徙動(dòng)曲線可表示為

圖2 無(wú)人機(jī)載SAR成像幾何Fig.2 UAV SAR imaging geometry

其中，(Δx(t),Δy(t),Δz(t))代表偏離勻速直線運(yùn)動(dòng)軌跡的三維運(yùn)動(dòng)誤差。

實(shí)際數(shù)值求解時(shí)，通過(guò)提取無(wú)人機(jī)非線性航跡的線性分量，基于無(wú)人機(jī)載SAR成像幾何，根據(jù)式(1)計(jì)算出該地形點(diǎn)P(xp,yp,zp)在線性分量下的投影點(diǎn)，作為迭代的初始解。

利用仿真地形和實(shí)測(cè)航跡對(duì)投影點(diǎn)位置獲取方法進(jìn)行驗(yàn)證。實(shí)測(cè)航跡如圖3所示，在線性分量附近無(wú)規(guī)律波動(dòng)。仿真地形如圖4所示，仿真地形中存在地形起伏。

圖3 實(shí)測(cè)航跡運(yùn)動(dòng)分量Fig.3 Motion component of the actual trajectory

圖4 仿真地形高度Fig.4 Simulation topography elevation

圖5所示為迭代解與真實(shí)解之差。真實(shí)解是指利用BP成像算法，對(duì)仿真地形進(jìn)行成像處理，將三維地形點(diǎn)投影至二維平面上，提取的各地形點(diǎn)的成像幅度峰值點(diǎn)的坐標(biāo)位置。迭代解是指結(jié)合航跡信息，利用式(4)對(duì)仿真地形上的各地形點(diǎn)進(jìn)行迭代處理，得到的投影點(diǎn)位置?？梢园l(fā)現(xiàn)，迭代出的近似解與真實(shí)投影點(diǎn)位置的差異在10-13m量級(jí)，完全可以忽略。因此，迭代得到的近似解可以作為地形點(diǎn)在二維平面上的投影點(diǎn)位置。

圖5 迭代解與真實(shí)投影點(diǎn)位置誤差Fig.5 Error between the iteration solution and the real projection point

與此同時(shí)，在該條航跡下，對(duì)比非線性航跡中的線性分量得到的初始解和迭代后得到的近似解，其二者差異在10-5m量級(jí)，如圖6所示。因此，在一定的約束條件下，可以用非線性航跡中線性分量下的投影點(diǎn)位置作為非線性航跡的真實(shí)投影點(diǎn)位置。

圖6 迭代解與初始解位置誤差Fig.6 Error between the iterative solution and initial solution

下面分析以線性分量下投影點(diǎn)位置代替非線性航跡投影點(diǎn)位置的邊界約束條件。實(shí)際上，在應(yīng)用式(4)進(jìn)行投影點(diǎn)位置獲取時(shí)，假設(shè)在成像過(guò)程中，地形點(diǎn)P(xp,yp,zp)在二維平面上的成像結(jié)果不散焦。此時(shí)，可利用非線性航跡中線性分量下的投影點(diǎn)位置作為非線性航跡的投影點(diǎn)位置。因此，可將非線性航跡下地形點(diǎn)P(xp,yp,zp)在二維平面上的成像結(jié)果不散焦的條件作為快速獲取投影點(diǎn)位置的邊界條件。

在進(jìn)行BP成像時(shí)，相干疊加時(shí)多普勒相位不能完全補(bǔ)償，存在相位誤差，嚴(yán)重時(shí)導(dǎo)致成像散焦。斜距誤差歷程可以表示為

忽略沿航跡方向的運(yùn)動(dòng)誤差 Δx(t)，將式(1)代入式(5)進(jìn)行近似化簡(jiǎn)可得

圖7 帶運(yùn)動(dòng)誤差的投影幾何Fig.7 Projection geometry with motion error

假設(shè)最大斜距誤差為 ΔRm，最大運(yùn)動(dòng)誤差為Δrm。通常認(rèn)為對(duì)于點(diǎn)目標(biāo)成像時(shí)，相位誤差大于π/4時(shí)會(huì)影響成像質(zhì)量。因此，可以假設(shè)滿足成像質(zhì)量的條件為

其中，λ為發(fā)射信號(hào)波長(zhǎng)。

當(dāng)航跡確定時(shí)，對(duì)于三維空間中任意目標(biāo)點(diǎn)，其投影點(diǎn)位置也確定，成像處理點(diǎn)目標(biāo)不散焦時(shí)，能量聚集于投影點(diǎn)位置。因此，當(dāng)式(10)成立時(shí)，目標(biāo)真實(shí)投影點(diǎn)位置等效于航跡線性分量對(duì)應(yīng)的投影點(diǎn)位置。

實(shí)際成像處理時(shí)，即使不考慮雷達(dá)位置的測(cè)量誤差，對(duì)于起伏地形成像時(shí)，波長(zhǎng)較小的高頻段SAR很難滿足式(10)所示的邊界條件，而工作在較低頻段的SAR則較容易滿足。圖8所示為地形與最大斜距誤差的關(guān)系。當(dāng)最大運(yùn)動(dòng)誤差和目標(biāo)地距一定時(shí)，最大斜距誤差隨目標(biāo)高度的增大而增大；當(dāng)目標(biāo)高度一定時(shí)，最大斜距誤差隨著地距的變大而減小，隨著最大運(yùn)動(dòng)誤差的增大而增大。當(dāng)運(yùn)動(dòng)誤差越小，目標(biāo)地距越大時(shí)，成像質(zhì)量受目標(biāo)高度的影響越小。圖8中λ為1 m。

圖8 不同最大運(yùn)動(dòng)誤差和目標(biāo)地距下最大斜距誤差與目標(biāo)高程之間的關(guān)系Fig.8 Relationship between the maximum slant range error and the target height with different the maximum motion error and the target ground distance

3 外部地形輔助分區(qū)的無(wú)人機(jī)載InSAR圖像配準(zhǔn)方法

3.1 無(wú)人機(jī)載InSAR偏移量建模

基于前文推導(dǎo)的無(wú)人機(jī)載SAR成像投影幾何模型與邊界條件可知，在較低頻段下，利用導(dǎo)航系統(tǒng)提供的雷達(dá)位置信息，成像投影點(diǎn)位置可以等效于非線性航跡中線性分量對(duì)應(yīng)的投影點(diǎn)位置。

圖9所示為無(wú)人機(jī)載InSAR干涉測(cè)量投影幾何。對(duì)于地形點(diǎn)P(xp,yp,zp)，分別用兩次航過(guò)對(duì)其進(jìn)行成像處理。兩次成像的合成孔徑中心分別為C1(xc1,yc1,zc1)和C2(xc2,yc2,zc2)，兩次成像的投影點(diǎn)位置分別為

圖9 無(wú)人機(jī)載InSAR干涉測(cè)量投影幾何Fig.9 UAV InSAR interferometric projection geometry

以其中一次航過(guò)為例，提取非線性航跡中的線性分量，假設(shè)其單位方向向量為 (m0,n0,p0)，直線上一點(diǎn)為 (x1,y1,z1)，直線方程可表示為

其中，t為直線方程的參數(shù)。

其中，t0為C1在直線方程中對(duì)應(yīng)的參數(shù)。

因此，C1(xc1,yc1,zc1)可表示為

弧的半徑r可表示為

其中，

將式(15)和式(17)代入式(16)化簡(jiǎn)可得

當(dāng)孔徑為線性且平行于飛行方向時(shí)，即xp=xc1,yc1=0,zc1=H時(shí)，式(18)與式(1)一致。

進(jìn)一步地，對(duì)于地形點(diǎn)P(xp,yp,zp)兩次航過(guò)的偏移量模型可以表示為

由式(19)可以看出，無(wú)人機(jī)飛行高度受限，與實(shí)際地形相差不大，星載InSAR或常規(guī)機(jī)載In-SAR測(cè)量中的平地近似假設(shè)不再成立。干涉圖像對(duì)的偏移量與地形相關(guān)，具有很強(qiáng)的空變性。特別是在長(zhǎng)基線模式下，平地的區(qū)域偏移量較小，而在有高程區(qū)域，偏移量可以達(dá)到數(shù)個(gè)像素。當(dāng)成像區(qū)域內(nèi)高程起伏較大時(shí)，偏移量出現(xiàn)明顯的空變現(xiàn)象，傳統(tǒng)的配準(zhǔn)方法失效。因此，基于多項(xiàng)式的全局偏移量擬合模型無(wú)法準(zhǔn)確描述真實(shí)的偏移量情況，需要對(duì)不同地形區(qū)域分別處理，才能獲取較好的配準(zhǔn)效果。

3.2 基于多級(jí)高程門限的圖像分區(qū)

基于BP成像算法，對(duì)兩次SAR回波數(shù)據(jù)在同一成像坐標(biāo)系下進(jìn)行成像處理[28]。由式(19)可知，當(dāng)航跡固定時(shí)，任意三維空間內(nèi)的目標(biāo)，其二維配準(zhǔn)偏移量也隨之確定。那么已知目標(biāo)的三維空間位置即可計(jì)算出該點(diǎn)的二維配準(zhǔn)偏移量。

首先，提取非線性航跡中的線性航跡分量，利用式(14)獲取各像素點(diǎn)合成孔徑中心位置 (xc,yc,zc)。然后，根據(jù)式(19)，對(duì)于任一像素點(diǎn)，其二維配準(zhǔn)偏移量與目標(biāo)高程一一對(duì)應(yīng)，設(shè)置偏移量門限，代入成像平面中各像素點(diǎn)的位置，從而獲取各像素點(diǎn)的高程門限。通常認(rèn)為，當(dāng)偏移量小于1個(gè)像素時(shí)，使用多項(xiàng)式擬合的方法可以較好地實(shí)現(xiàn)全局偏移量的擬合。因此，需設(shè)置多級(jí)高程門限進(jìn)行圖像分區(qū)，而后對(duì)子區(qū)域內(nèi)進(jìn)行多項(xiàng)式擬合，才能更準(zhǔn)確地表示區(qū)域內(nèi)的偏移量情況。

最后，利用地形數(shù)據(jù)與獲取的多級(jí)高程門限相比，將成像平面進(jìn)行區(qū)域劃分，劃分為低于門限的接近成像平面的高程區(qū)域和高于門限的其他高程區(qū)域。對(duì)分塊后的結(jié)果進(jìn)行連通域分析，剔除面積較小的區(qū)域，得到最終的圖像分區(qū)結(jié)果。

實(shí)際處理過(guò)程中，在外部地形的精度不足時(shí)，不能保證分區(qū)的完全準(zhǔn)確。本文對(duì)圖像分區(qū)時(shí)采用亞像素級(jí)的偏移量作為分區(qū)門限，實(shí)際各區(qū)不可避免地會(huì)存在部分像素點(diǎn)的偏移量超過(guò)門限，但是本文方法核心思想是將全圖的偏移量擬合問(wèn)題，轉(zhuǎn)化為若干個(gè)子區(qū)域的偏移量在約束條件下的聯(lián)合擬合求解問(wèn)題，分區(qū)誤差對(duì)局部偏移量的多項(xiàng)式擬合影響有限。

3.3 帶約束的子塊多項(xiàng)式擬合

無(wú)人機(jī)載InSAR的二維配準(zhǔn)偏移量與地形有關(guān)。對(duì)于整個(gè)場(chǎng)景區(qū)域，其二維配準(zhǔn)偏移量應(yīng)該是連續(xù)的。進(jìn)行圖像分區(qū)后，如果直接采用常規(guī)的多項(xiàng)式擬合方法來(lái)獲取各個(gè)子區(qū)域內(nèi)的偏移量，由于邊界處的擬合偏移量往往過(guò)大或過(guò)小，易導(dǎo)致各子區(qū)域邊界連接處偏移量出現(xiàn)較大“跳變”，因此本文提出了帶約束的子塊多項(xiàng)式擬合方案。

以一階多項(xiàng)式擬合為例，假設(shè)圖像分為N個(gè)子區(qū)域，并且子區(qū)域內(nèi)M個(gè)參考點(diǎn)的二維坐標(biāo)和二維配準(zhǔn)偏移量已知，可以得到

其中，(xm,n,ym,n)代表第n個(gè)子區(qū)域第m個(gè)參考點(diǎn)在雷達(dá)圖像中的二維坐標(biāo)，(Δxm,n,Δym,n)代表第n個(gè)子區(qū)域第m個(gè)參考點(diǎn)的二維配準(zhǔn)偏移量。an1,bn1,cn1,an2,bn2,cn2分別為待估計(jì)參數(shù)。

對(duì)于區(qū)域邊緣的像素，需要附加約束條件，使二維配準(zhǔn)偏移量連續(xù)，即需滿足：

其中，(xm,nq,ym,nq)代表第n個(gè)子區(qū)域和第q個(gè)子區(qū)域分界處的第m個(gè)參考點(diǎn)的二維坐標(biāo)，(Δxm,nq,Δym,nq)代表第n個(gè)子區(qū)域和第q個(gè)子區(qū)域分界處的第m個(gè)參考點(diǎn)的二維配準(zhǔn)偏移量。an1,bn1,cn1,an2,bn2,cn2,aq1,bq1,cq1,aq2,bq2,cq2分別為待估計(jì)參數(shù)。

參考點(diǎn)選擇時(shí)，在子區(qū)域內(nèi)選擇主圖像中信噪比較高的點(diǎn)；同時(shí)，選擇所有子區(qū)域邊界處的點(diǎn)作為參考點(diǎn)。

如圖10所示，以主圖像的參考點(diǎn) (i,j)為中心點(diǎn)，確定一定大小的匹配窗口，在輔圖像上相同位置為中心選取比匹配窗口大的搜索窗口。沿著行向和列向，在搜索窗中順序移動(dòng)匹配窗，同時(shí)計(jì)算每個(gè)位置處的相干系數(shù)，選取相干系數(shù)最大的位置作為同名點(diǎn)的位置，并得到同名點(diǎn)在主輔圖像中的二維偏移量。

圖10 匹配窗和搜索窗示意圖Fig.10 Schematic diagram of match window and search window

根據(jù)式(20)和式(21)，構(gòu)建目標(biāo)函數(shù)聯(lián)合求解，目標(biāo)函數(shù)可表示為

引入Lagrange算子，求解目標(biāo)函數(shù)。引入Lagrange算子后的目標(biāo)函數(shù)可表示為

對(duì)于式(23)，可利用最小二乘法聯(lián)合求解模型參數(shù)。根據(jù)圖像中所有點(diǎn)的二維坐標(biāo)，即可以實(shí)現(xiàn)全局偏移量的估計(jì)。最后，對(duì)輔圖像進(jìn)行插值，即可獲得精配準(zhǔn)后的圖像。圖11所示為基于地形輔助的無(wú)人機(jī)載InSAR圖像分區(qū)配準(zhǔn)方法流程圖。

圖11 基于地形輔助的無(wú)人機(jī)載InSAR圖像分區(qū)配準(zhǔn)方法流程圖Fig.11 Flowchart of topography-assisted UAV InSAR registration method with image partition

4 實(shí)測(cè)數(shù)據(jù)驗(yàn)證

4.1 P波段無(wú)人機(jī)載SAR系統(tǒng)

為驗(yàn)證本文提出的基于圖像分區(qū)的無(wú)人機(jī)載InSAR圖像配準(zhǔn)方法的有效性，采用北京理工大學(xué)自研的P波段無(wú)人機(jī)載SAR原理樣機(jī)開展驗(yàn)證實(shí)驗(yàn)。系統(tǒng)照片如圖12所示，系統(tǒng)參數(shù)如表1所示。

表1 系統(tǒng)參數(shù)Tab.1 System parameters

圖12 P波段無(wú)人機(jī)載SAR原理樣機(jī)Fig.12 P-band UAV SAR system

4.2 實(shí)測(cè)數(shù)據(jù)分析

4.2.1 北京平谷金海湖機(jī)場(chǎng)實(shí)驗(yàn)

4.2.1.1 實(shí)驗(yàn)信息

利用P波段無(wú)人機(jī)載SAR在北京平谷區(qū)金海湖機(jī)場(chǎng)(40.19°N,117.23°E)附近進(jìn)行了無(wú)人機(jī)載SAR干涉測(cè)量實(shí)驗(yàn)。圖13所示為實(shí)驗(yàn)場(chǎng)景，藍(lán)色實(shí)線表示無(wú)人機(jī)飛行軌跡，紅色實(shí)線圈出的區(qū)域?yàn)槔走_(dá)觀測(cè)區(qū)域。觀測(cè)區(qū)域主要由一些山丘組成。航跡總長(zhǎng)480 m，無(wú)人機(jī)平臺(tái)飛行速度為8 m/s，單航過(guò)飛行時(shí)長(zhǎng)為60 s。

圖13 金海湖機(jī)場(chǎng)實(shí)驗(yàn)場(chǎng)景Fig.13 Jinhai Lake airport experiment scene

兩次航過(guò)的航跡情況如圖14所示，可以看到水平向和垂直向偏離直線運(yùn)動(dòng)的分量很小，最大運(yùn)動(dòng)誤差不超過(guò)0.2 m。圖15為插值前后的外部數(shù)字高程模型(Digital Elevation Model,DEM)。原始DEM來(lái)自SRTM官方數(shù)據(jù)，空間分辨率為90 m，需將DEM結(jié)果插值到成像像素分辨單元下，這里插值后空間分辨率為1 m。結(jié)合圖15，其相位誤差分量滿足式(10)所示條件，成像結(jié)果可以聚焦。

圖14 金海湖場(chǎng)景兩次航過(guò)航跡運(yùn)動(dòng)分量Fig.14 Motion components of two trajectories in Jinhai Lake

圖15 金海湖場(chǎng)景SRTM-DEM結(jié)果Fig.15 DEM results from SRTM of Jinhai Lake

圖16為金海湖場(chǎng)景中P波段無(wú)人機(jī)載SAR成像結(jié)果，其中以最大幅值為參考對(duì)圖像進(jìn)行了歸一化處理。利用高精度的位置測(cè)量系統(tǒng)，基于BP成像算法能夠很好地補(bǔ)償多普勒相位，實(shí)現(xiàn)高質(zhì)量聚焦成像[29]。

圖16 金海湖場(chǎng)景成像結(jié)果Fig.16 The imaging result of Jinhai Lake

4.2.1.2 配準(zhǔn)結(jié)果分析

對(duì)于垂直基線為12 m的兩個(gè)航過(guò)的雷達(dá)圖像，獲取的相干系數(shù)圖與濾波后[30]的干涉相位圖如圖17所示。可以看出，高密度干涉條紋所在區(qū)域的相干性較低。這些嚴(yán)重失相干區(qū)域的高程較大，即有地形的區(qū)域。

圖17 金海湖場(chǎng)景12 m基線配準(zhǔn)前相干系數(shù)圖和干涉相位圖Fig.17 Coherence coefficient diagram and interferogram with 12 m baseline before registration of Jinhai Lake

結(jié)合兩個(gè)航過(guò)的航跡信息，根據(jù)式(19)，設(shè)置偏移量門限為0.5個(gè)像素單元，可以得到如圖18(a)所示的門限結(jié)果。可以看到，距離向越近的目標(biāo)，斜距越短，對(duì)高程也越敏感，在很低的高程即會(huì)出現(xiàn)較大的配準(zhǔn)偏移量。圖18(b)表示偏移量大于1.5個(gè)像素單元時(shí)所需的高程門限，可以看到此時(shí)的二級(jí)高程門限已經(jīng)高于圖15所示的DEM結(jié)果，此時(shí)無(wú)需使用二級(jí)高程門限。以高程門限將所示的外部地形結(jié)果進(jìn)行分區(qū)，并進(jìn)行連通域處理，剔除面積較小的部分，得到圖19所示的圖像分區(qū)結(jié)果，可以看到，地形區(qū)域均有效地被提取出來(lái)。

圖18 金海湖場(chǎng)景高程門限Fig.18 Threshold results of Jinhai Lake

圖19 金海湖場(chǎng)景圖像分區(qū)結(jié)果Fig.19 The segmentation result of Jinhai Lake

圖20所示為圖像分區(qū)結(jié)果與相關(guān)系數(shù)、干涉相位的對(duì)比結(jié)果?？梢钥闯?，由于地形導(dǎo)致的圖像失配，被提取出來(lái)的地形區(qū)域普遍存在相干性低、干涉相位質(zhì)量差的問(wèn)題。

圖20 金海湖場(chǎng)景地形分區(qū)結(jié)果Fig.20 Comparison between the segmentation results of Jinhai Lake

圖像分區(qū)后，將各區(qū)域內(nèi)高信噪比和區(qū)域邊界處的像素點(diǎn)作為參考點(diǎn)，利用基于窗口滑動(dòng)的配準(zhǔn)偏移量提取方法對(duì)所有參考點(diǎn)的二維配準(zhǔn)偏移量進(jìn)行提取，并進(jìn)一步篩選出高相干性的參考點(diǎn)。圖21為精配準(zhǔn)結(jié)果，可以看到在地形起伏區(qū)域干涉相位質(zhì)量大幅提升。

圖21 金海湖場(chǎng)景本文所提方法的精配準(zhǔn)結(jié)果Fig.21 Fine registration results of the proposed method of Jinhai Lake

4.2.1.3 比較與評(píng)估

為了進(jìn)一步說(shuō)明本文所提方法的有效性，對(duì)比了幾種常見的InSAR配準(zhǔn)方法[31]，分別是平均波動(dòng)函數(shù)、最大頻譜法和最大互相關(guān)法。利用3種方法配準(zhǔn)后的相干系數(shù)圖和干涉相位圖如圖22所示。

圖22 金海湖場(chǎng)景常見InSAR方法配準(zhǔn)后的相干系數(shù)圖和干涉相位圖Fig.22 Coherence coefficient diagrams and interferograms after registration of common InSAR methods of Jinhai Lake

從圖22可以看出，常見的InSAR配準(zhǔn)方法雖然能夠提升一定的相干系數(shù)和干涉圖質(zhì)量，但對(duì)于起伏地形區(qū)域提升有限。相比常見的InSAR配準(zhǔn)方法，本文所提的配準(zhǔn)方法，以地形分區(qū)為先導(dǎo)，對(duì)地形分區(qū)后的不同高程子區(qū)域進(jìn)行分區(qū)配準(zhǔn)，保證了各子區(qū)域的干涉圖質(zhì)量。如圖23所示，本文所提方法較常規(guī)方法有一定程度的提升，平均相干系數(shù)與殘差點(diǎn)數(shù)量如表2所示，可以看到殘差點(diǎn)數(shù)量大幅減少，證明了本文所提方法的有效性。圖24為所提方法反演出的地形。

表2 金海湖場(chǎng)景干涉圖評(píng)價(jià)Tab.2 Evaluation of interferograms of Jinhai Lake

圖23 金海湖場(chǎng)景不同方法相干系數(shù)分布曲線Fig.23 Coherence coefficient distribution with different methods of Jinhai Lake

圖24 所提方法反演的地形結(jié)果Fig.24 Topography result of the proposed method

4.2.2 重慶奉節(jié)老林溝實(shí)驗(yàn)

4.2.2.1 實(shí)驗(yàn)信息

圖25所示為重慶奉節(jié)縣老林溝(30.89°N,109.48°E)實(shí)驗(yàn)場(chǎng)景，藍(lán)色實(shí)線表示無(wú)人機(jī)飛行軌跡，紅色實(shí)線圈出的區(qū)域?yàn)槔走_(dá)觀測(cè)區(qū)域。觀測(cè)區(qū)域主要由一舌頭形山體組成，如黃色虛線所示。航跡總長(zhǎng)480 m，無(wú)人機(jī)平臺(tái)飛行速度為8 m/s，單航過(guò)飛行時(shí)長(zhǎng)為60 s。

圖25 老林溝實(shí)驗(yàn)場(chǎng)景Fig.25 Laolin Gou experiment scene

兩次航過(guò)的航跡情況如圖26所示，可以看到水平向和垂直向偏離直線運(yùn)動(dòng)的分量同樣很小，最大運(yùn)動(dòng)誤差不超過(guò)0.2 m。如圖27所示為插值前后的DEM，插值后空間分辨率仍為1 m。結(jié)合圖27所示的DEM，其相位誤差分量滿足式(10)所示條件，成像結(jié)果可以聚焦。圖28為P波段無(wú)人機(jī)載SAR成像結(jié)果，同樣以最大幅值為參考對(duì)圖像進(jìn)行了歸一化處理。從圖28可以看到，圖像的左邊區(qū)域有部分陰影區(qū)域，信噪比較低。

圖26 老林溝場(chǎng)景兩次航過(guò)航跡運(yùn)動(dòng)分量Fig.26 Motion components of two trajectories in Laolin Gou

圖27 老林溝場(chǎng)景SRTM-DEM結(jié)果Fig.27 DEM results from SRTM of Laolin Gou

圖28 老林溝場(chǎng)景成像結(jié)果Fig.28 The imaging result of Laolin Gou

4.2.2.2 配準(zhǔn)結(jié)果分析

對(duì)于垂直基線為13 m的兩個(gè)航過(guò)的雷達(dá)圖像，獲取的相干系數(shù)圖與干涉相位圖如圖29所示?？梢钥闯?，由于地形相對(duì)較高，除了接近平地的區(qū)域有一定的相干性和干涉條紋外，其余區(qū)域內(nèi)幾乎完全失相干，看不到干涉條紋。

圖29 老林溝場(chǎng)景13 m基線配準(zhǔn)前相干系數(shù)圖和干涉相位圖Fig.29 Coherence coefficient diagram and interferogram with 13 m baseline before registration of Laolin Gou

結(jié)合兩個(gè)航過(guò)的航跡信息，根據(jù)式(19)，設(shè)置偏移量門限為0.5個(gè)像素單元，可以得到如圖30(a)所示的門限結(jié)果。圖30(b)表示偏移量大于1.5個(gè)像素單元時(shí)所需的高程門限，可以看到此時(shí)的二級(jí)高程門限已經(jīng)包含了絕大部分的DEM結(jié)果，無(wú)需構(gòu)建三級(jí)高程門限。以高程門限將所示的外部地形結(jié)果進(jìn)行分區(qū)，并進(jìn)行連通域處理剔除面積較小的部分，得到圖31所示的圖像分區(qū)結(jié)果。對(duì)照?qǐng)D29所示的干涉相位圖和相干系數(shù)圖，可以看出低于一級(jí)門限的區(qū)域尚能保持一定的相干性并有一定程度的干涉條紋，而高于一級(jí)門限區(qū)域幾乎失相干，說(shuō)明地形區(qū)域均有效地被提取出來(lái)。

圖30 老林溝場(chǎng)景高程門限Fig.30 Threshold results of Laolin Gou

圖31 老林溝場(chǎng)景圖像分區(qū)結(jié)果Fig.31 Segmentation result of Laolin Gou

圖像分區(qū)后，利用所提方法進(jìn)行圖像配準(zhǔn)。圖32所示為精配準(zhǔn)結(jié)果，可以看到在地形起伏區(qū)域干涉相位質(zhì)量大幅提升。

圖32 老林溝場(chǎng)景本文所提方法的精配準(zhǔn)結(jié)果Fig.32 Fine registration results of the proposed method of Laolin Gou

4.2.2.3 比較與評(píng)估

同樣利用4.2.1.3節(jié)所述的3種常見InSAR配準(zhǔn)方法進(jìn)行對(duì)比，得到的相干系數(shù)圖和干涉相位圖如圖33所示。

圖33 老林溝場(chǎng)景常見InSAR方法配準(zhǔn)后的相干系數(shù)圖和干涉相位圖Fig.33 Coherence coefficient diagrams and interferograms after registration of common InSAR methods of Laolin Gou

從圖33可以看出，在整體信噪比和相干系數(shù)較低的情況下，AFF和MSF方法受到噪聲的影響嚴(yán)重，靈敏度大幅降低；MCC方法雖然能夠提升一定的相干系數(shù)和干涉圖質(zhì)量，但對(duì)于起伏地形區(qū)域提升有限。相比常規(guī)方法，本文所提的配準(zhǔn)方法更好地提升了干涉圖質(zhì)量。如圖34所示，本文所提方法較常規(guī)方法有一定程度的提升，平均相干系數(shù)與殘差點(diǎn)數(shù)量如表3所示，可以看到殘差點(diǎn)數(shù)量大幅減少，證明了本文所提方法的先進(jìn)性。

表3 老林溝場(chǎng)景干涉圖評(píng)價(jià)Tab.3 Evaluation of interferograms of Laolin Gou

圖34 老林溝場(chǎng)景不同方法相干系數(shù)分布曲線Fig.34 Coherence coefficient distribution with different methods of Laolin Gou

5 結(jié)語(yǔ)

本文面向無(wú)人機(jī)平臺(tái)干涉測(cè)量，重點(diǎn)研究了無(wú)人機(jī)平臺(tái)本身給圖像配準(zhǔn)帶來(lái)的影響和解決方案。對(duì)于復(fù)雜地形，無(wú)人機(jī)載InSAR不同同名像素之間的偏移量大且具有明顯的空變特性，常規(guī)的多項(xiàng)式擬合的變換方法難以取得良好的全局配準(zhǔn)效果。針對(duì)這一問(wèn)題，本文提出了一種基于圖像分區(qū)的無(wú)人機(jī)載InSAR圖像配準(zhǔn)方法，利用外部粗地形，并結(jié)合航跡信息對(duì)圖像進(jìn)行分區(qū)處理，將全圖的偏移量擬合問(wèn)題，轉(zhuǎn)化為若干個(gè)子區(qū)域的偏移量在約束條件下的聯(lián)合擬合求解問(wèn)題，實(shí)現(xiàn)圖像對(duì)之間的精配準(zhǔn)。

該方法首先基于航跡信息生成高程門限，再利用外部粗地形對(duì)測(cè)量區(qū)域進(jìn)行圖像分區(qū)處理。之后對(duì)各子區(qū)域內(nèi)的偏移量進(jìn)行多項(xiàng)式擬合，并對(duì)區(qū)域之間邊界處的偏移量施加約束條件，構(gòu)建帶約束的多項(xiàng)式擬合目標(biāo)函數(shù)，引入Lagrange算子，利用最小二乘法聯(lián)合求解模型參數(shù)，并進(jìn)行圖像重采樣獲取精配準(zhǔn)后的圖像。基于P波段無(wú)人機(jī)載InSAR獲取的實(shí)測(cè)數(shù)據(jù)，驗(yàn)證了該方法的有效性。

本文還有部分問(wèn)題尚未解決。首先，本文所用數(shù)據(jù)集十分有限，無(wú)法保證所提方法能夠適應(yīng)不同復(fù)雜地形條件。其次，我們暫未獲取到對(duì)應(yīng)場(chǎng)景的較為準(zhǔn)確的DEM，因此未對(duì)生成地形的精度進(jìn)行評(píng)估。最后，對(duì)于所提方法自身，由于實(shí)驗(yàn)條件所限，也未對(duì)配準(zhǔn)精度直接進(jìn)行定量評(píng)估。后續(xù)，課題組將繼續(xù)致力于無(wú)人機(jī)載InSAR測(cè)量領(lǐng)域的關(guān)鍵技術(shù)研究，獲取盡可能多的地形場(chǎng)景驗(yàn)證所提方法的魯棒性和適應(yīng)性。

利益沖突所有作者均聲明不存在利益沖突

Conflict of InterestsThe authors declare that there is no conflict of interests