一種有效的機載雙頻干涉地形高程重建方法

劉華有 鄭明潔 張 衡* 王 宇 秦小芳

①(中國科學院電子學研究所 北京 100190)

②(中國科學院大學 北京 100049)

1 引言

合成孔徑雷達干涉(Interferometry Synthetic Aperture Radar, InSAR)技術，是廣泛應用在地理測量中的一種雷達遙感技術。利用同一場景的兩幅合成孔徑雷達(Synthetic Aperture Radar, SAR)復數據，提取相位信息來獲取地面的高程信息[1,2]?？梢詰糜诘匦螠y繪、形變監測、海洋應用、以及地球動力學應用[3–5]。

傳統的單波段InSAR系統只需兩幅SAR圖像形成一幅干涉圖，即可實現高程反演，其處理流程[1]主要包括復圖像配準、干涉相位去平地、干涉相位濾波、相位解纏、數字高程模型(Digital Elevation Model, DEM)重建等步驟。雖然單波段InSAR系統數據獲取相對容易，處理流程相對簡單[3]，但是單波段相位解纏必須滿足Iton假設[6]，無法實現地形起伏比較大的高程反演。90年代以后多通道InSAR技術快速發展，多通道InSAR技術主要包括多頻率干涉和多基線干涉兩種。1994年，Xu Wei等人[7]提出了中國剩余定理(Chinese Reminder Theorem,CRT)、投影法和線性組合法3種多通道聯合相位解纏算法。2000年，Pascazio 和Schirinzi[8]提出將大帶寬的SAR數據進行分割來獲取多頻率干涉數據，并且使用最大似然估計(Maximum Likelihood,ML)來實現多頻率聯合相位解纏。2003年，Hoge等人[9]提出了一種基于頻域奇異值分解(Singular Value Decomposition, SVD)的配準方法，可以很高效的實現圖像亞像素級配準。2004年，Ferraiuolo等人[10]提出了一種高斯能量項的最大后驗估計(Maximum A Posterior, MAP)的算法來多通道聯合高程重建，該方法使用馬爾科夫隨機場來建模地形高程的先驗分布，并通過仿真的多頻率干涉驗證了算法的有效性。2009年，Ferraioli等人[11]提出了一種基于圖割法的變分能量項最大后驗估計(Maximum Posteriori Estimation of Total Variational Energy Term, TV-MAP)的多通道聯合高程重建方法。Ferraiuolo和 Meglio等人[12]深入分析最大似然估計和最大后驗估計等兩種多通道聯合相位解纏算法，證明了最大后驗估計法魯棒性優于最大似然估計。2011年，于翰雯等人[13]提出了一種基于聚類分析的多通道聯合相位解纏繞算法，該方法將具有相同模糊矢量的像素點聚成一類，隨后逐類進行聯合解纏。Deledalle等人[14]提出了一種非局部參數估計(Nonlocal Interferogram Estimation, NL-InSAR)的方法，可以實現對干涉相位、相干系數、幅度圖的精確估計。2013年，袁志輝[15]提出了一種基于閉合形式魯棒性中國剩余定理(CRT)的多通道聯合高程重建算法。2016年，曾濤[16]等人提出了結合最大似然估計相位梯度和枝切法的多頻率聯合相位解纏算法。2017年，斯奇等人[17]提出了一種基于最大后驗框架的聚類分析多基線干涉SAR高度重建算法。

雙頻InSAR技術利用不同頻段的兩幅干涉圖進行聯合高程反演，可以有效地提取地形復雜區域的高程信息。為實現高精度、強魯棒性的高程反演，雙頻干涉數據處理需要解決不同頻段之間的非相干配準、干涉相位降噪、聯合相位解纏等一系列問題。針對這些問題，結合現有圖像配準、相位濾波、多通道InSAR技術等研究進展[18,19]，本文提出了一種有效的雙頻干涉SAR地表高程重建方法。該方法對常規處理流程中的關鍵步驟進行了改進，首先采用NL-InSAR技術對幅度圖、相干系數、干涉相位通過迭代進行精確估計，利用各個波段的幅度信息來實現不同波段的干涉相位的配準。然后采用聚類分析技術對聯合解纏相位標記有效點和噪點，并利用有效點對聯合解纏相位進行均值濾波。通過這些改進，提高了高程重建效果與高程反演精度。目前該方法已經應用在實際機載雙頻InSAR系統中，取得了較好的結果。

2 高程重建算法

本文針對常規的雙頻干涉處理流程中的關鍵步驟進行了改進，改進后的流程圖如圖1所示，主要包括頻段內圖像配準、去平地相位、干涉參數精確估計、頻段間干涉相位配準、雙頻聯合相位解纏、解纏相位噪點剔除等步驟。下面通過比較分析，將給出各個步驟的處理方法，并在機載實測雙頻干涉處理中驗證其有效性。

2.1 圖像配準

圖像配準是通過調整主、輔SAR復圖像，使得兩幅圖像上相同的點對應地面場景上的同一點。配準的步驟如下[20,21]：

步驟1 主輔圖像選取若干控制點，采用相關法計算控制點處粗偏移量，精度為1個像素。

計算控制點處的偏移量可以采用空域相關函數法的計算基于窗口的配準方法，分別在主圖像s1和輔圖像s2以控制點為中心建立匹配窗和搜索窗?；瑒铀阉鞔扒蟮么翱诘南嚓P函數值，相關函數最大的滑動位置即為圖像的粗偏移量(u，v)；

步驟2 利用相關法求得圖像粗偏移量，并修正輔圖像控制點的位置，為得到亞像素級的精配準，可以通過空域相干系數法[21]、傅里葉變換法[21]、頻域SVD法[9]求得控制點處的精確偏移量；

步驟3 在求得控制點處精確偏移量后，可以利用二項式進行擬合：

其中， (x，y)為主圖像的控制點坐標，通過最小二乘法可以得擬合系數，再利用擬合系數得到輔圖像整幅圖像配準后精確的位置，通過線性插值把輔圖像插值到求得的精確位置來實現精配準[21]。

下面采用某機載雙頻干涉系統飛行的數據為例來分析圖像配準技術。

空域相關法、傅里葉變換法、頻域SVD分解法這3種方法都具有較高的配準精度，本文分別采用了這3種方法對上述圖像進行了配準，配準結果如圖2所示。由于頻域SVD分解法運算效率最高，因此機載雙頻干涉高重建過程中采用頻域SVD法配準。

2.2 去平地

配準后的干涉相位由于平地相位的存在，使干涉相位纏繞量較大，干涉條紋過密，會影響后續的相位解纏。因此在相位解纏之前，需要去除平地相位，以便降低聯合相位解纏的難度。去平地方法主要有頻譜法和軌道法兩種，頻譜法是將干涉復圖像變換到頻域，找出頻譜中最大值對應的頻率并搬移到零頻。軌道法是根據軌道參數計算平地相位。兩種方法的處理結果如圖3所示。由于在實測數據處理中頻譜法并不能有效地去除平地相位，因此，本文采用軌道法去平地。

2.3 干涉參數精確估計

為了提高不同頻段配準和聯合解纏的精度，采用NL-InSAR技術分別對各個頻段的干涉相位、相干系數、幅度進行精確估計[14]，再進行不同波段配準，可以有效提高不同頻段之間干涉相位的配準精度，并且實現了對干涉相位去噪的同時較好地保持了干涉相位的紋理信息。圖4給出了NLInSAR濾波前后的干涉相位對比情況。從圖中可以看出，NL-InSAR在去除干涉相位存在的斑點噪聲。NL-InSAR濾波前的幅度圖和濾波后的幅度圖如圖5所示，對比表明NL-InSAR在保持SAR圖像邊緣細節信息的前提下能較好地去除SAR圖像的斑點噪聲。NL-InSAR估計前的相干系數如圖6所示，可以看出，NL-InSAR估計后的相干系數的可靠性更強，可以有效輔助后續雙頻聯合解纏。

2.4 頻段間干涉圖配準

交叉頻段干涉相位配準精度往往較低，本文采用的方法為首先進行同頻段圖像配準、濾波、去平地，通過利用各自主圖像的幅度信息擬合出多項式系數，進而來實現頻段間干涉相位的配準。

由于不同頻段的主圖像幅度存在較多的斑點噪聲，造成不同頻段的主圖像幅度圖的相似性較差。由圖5可知，NL-InSAR可以有效地對SAR圖像的幅度圖進行去噪。圖7(a)為不進行NL-InSAR幅度濾波的C波段主圖像幅度圖，圖7(b)為不進行NLInSAR濾波的X波段主圖像幅度。可以看出不進行NL-InSAR濾波的不同波段主圖像幅度圖較多的斑點噪聲，造成不同波段的主圖像幅度圖之間相似性較差。圖8(a)為采用NL-InSAR對幅度濾波后的C波段主圖像幅度，圖8(b)為采用NL-InSAR對幅度濾波后的X波段主圖像。NL-InSAR可以很好地去除SAR圖像幅度圖的斑點噪聲，同時較好地保持了圖像的邊緣細節信息，提高了不同波段主圖像幅度圖之間的相似性。因此本文采用經過NL-InSAR濾波后的兩個頻段的主圖像幅度擬合出式(1)的多項式系數，從而實現對兩個頻段的干涉相位圖之間的配準。針對實測數據不同頻段配準前后的不同波段之間相關系數圖如圖9所示。圖9(a)所示的配準前X波段和C波段的相關系數較低，均值為0.46，圖9(b)所示的不采用NL-InSAR對幅度進行濾波的配準后的相關系數相對較高，均值為0.90。圖9(c)所示為采用NL-InSAR對幅度進行濾波的配準好的相關系數最高，均值為0.94，驗證了本文提出的不同波段干涉相位配準方法的有效性。

2.5 多通道聯合相位解纏

目前，多頻段聯合相位解纏算法大致有基于統計信號處理和非統計信號處理兩類，本文仿真對比了非統計信號處理類的魯棒性中國剩余定理(CRT)[15]，統計信號處理類的最大似然估計(ML)[8]和變分最大后驗估計(TV-MAP)[11]。在實測數據處理中，雙頻聯合相位解纏算法對算法魯棒性要求較高，本文選用魯棒性較強的變分最大后驗估計作為雙頻聯合解纏的處理方案。最大后驗估計相位解纏就利用觀測相位Ψ并且加入真實相位的先驗分布P(Φ)來估計真實相位Φ，即：

假設所有像素點之間是相互獨立的，那么：

對于N0個波段相位聯合解纏情況下，相應似然函數為：

其中單個波段的觀測相位的似然函數為：

我們可以使用馬爾科夫隨機場來描述圖像的局部特性，根據MarKov-Gibbs等價性可以用Gibbsian分布來描述馬爾科夫先驗分布[8]：

變分能量項最大后驗估計使用的是一種在圖像領域常用的整體變分(TV)模型，其能量函數定義為：

式(7)中，β是一個固定的常量，用來調節真實相位的光滑性與觀測相位和真實相位之間相似性的平衡的常數，Np為像素點p的鄰域的集合，wp，q表征鄰域連通性，若Np為4鄰域模型，其值為1。

把式(7)代入式(6)結合式(2)–式(5)并取對數且添加負號，可以得到變分能量項最大后驗估計的在TV模型下的相位估計公式為：

上述優化模型可以通過圖割理論求解[11]，通過求解上述優化模型可以求解得到各個波段干涉相位的纏繞量和聯合解纏繞相位。

對添加相同噪聲的仿真數據，采用上述3種方法分別進行聯合解纏結果如圖10所示，用來仿真干涉相位的原始相位相減的殘差圖如圖11所示。通過上述仿真分析，進一步驗證了變分能量項最大后驗估計用于雙頻聯合解纏的魯棒性。

2.6 聚類均值濾波

模糊矢量[13,22]定義為：

其中，Kc為C波段干涉相位的纏繞量，Kx為X波段干涉相位的纏繞量。將圖像各個點的模糊矢量求和得到：

如圖12所示，將兩種標記結果相與，得到用于加權均值濾波的標記結果，對聯合解纏結果進行如下加權均值濾波：

3 實測數據分析

下面采用實際飛行數據對本文方法進行驗證，同時對比了單波段和雙波段的高程重建效果。實測數據是由中國科學院電子學研究所航天微波遙感系統部提供的某場景雙波段(C波段中心頻率5.4 GHz和X波段中心頻率9.6 GHz)的機載雙頻干涉數據。

實測數據包括同一場景的C波段和X波段主、輔SAR圖像復數據，處理數據的場景如圖16所示，為陜西澄城地區的某山區，其中地形起伏比較大區域為圖中右上角的高架，不同波段配準之后的C波段干涉圖和X波段干涉圖分別如圖17所示。

單波段使用最小代價流(MCF)解纏，處理結果如圖18所示。使用最大似然估計方法雙頻聯合解纏結果如圖19所示，由于波段數較少，最大似然估計魯棒性較差，完全無法實現正確的相位解纏。按照本文方法，雙頻聯合解纏結果如圖20所示。對圖中B區域紅框內地形起伏比較大的高架區域解纏結果放大對比如圖21所示，從圖中可以看出在地形起伏較大的區域，單波段也無法實現正確解纏，而使用本文提出的雙頻聯合解纏算法可以正確解纏。

根據軌道參數反演的DEM如圖22所示，圖中下半部分為上半部分中沿虛線的地形走勢，通過對比可知，無論是X波段還是C波段，都無法正確反演高架橋高程，只有聯合X波段和C波段，才能正確反演。這也更加明確地說明了，在地形起伏比較大的區域，只有雙頻聯合才可以有效地反演高程信息。

由于本文所用的實測數據場景中沒有準確的地本文采用的方法是將谷歌光學地圖面控制點，本文采用的方法是將谷歌光學地圖作為主圖像和SAR圖像的幅度圖進行配準，得到配準的多項式系數，利用擬合的多項式系數實現對DEM進行配準。在光學圖像上選取5個控制點對DEM進行精度分析。單波段反演高程的均方根誤差(Root Mean Square Error, RMSE)為4.5826 m，且存在部分地區高程反演失敗，而雙頻聯合反演的RMSE為1.3476 m，可以看出本文提出雙頻聯合高程重建方法精度優于單波段高程重建。

4 結束語

本文提出了一種有效的雙頻干涉處理方法，分析并給出了算法中各個步驟的處理方案。針對雙頻聯合處理中的雙頻聯合配準、雙頻解纏繞等關鍵技術難題，給出了有效的解決方案，提高了雙頻聯合反演DEM的精度和魯棒性。采用常規單波段方法和本文方法對機載實際飛行數據處理，處理結果表明，單波段使用MCF進行相位解纏，無法實現大起伏地形區域的高程重建；而采用本文提出的處理流程能夠準確地反演大起伏地形區域的高程，驗證了本文所提出的雙頻干涉處理方法的有效性。