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

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

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

②(中國(guó)科學(xué)院大學(xué) 北京 100049)

1 引言

合成孔徑雷達(dá)干涉(Interferometry Synthetic Aperture Radar, InSAR)技術(shù)，是廣泛應(yīng)用在地理測(cè)量中的一種雷達(dá)遙感技術(shù)。利用同一場(chǎng)景的兩幅合成孔徑雷達(dá)(Synthetic Aperture Radar, SAR)復(fù)數(shù)據(jù)，提取相位信息來(lái)獲取地面的高程信息[1,2]?？梢詰?yīng)用于地形測(cè)繪、形變監(jiān)測(cè)、海洋應(yīng)用、以及地球動(dòng)力學(xué)應(yīng)用[3–5]。

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

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

2 高程重建算法

本文針對(duì)常規(guī)的雙頻干涉處理流程中的關(guān)鍵步驟進(jìn)行了改進(jìn)，改進(jìn)后的流程圖如圖1所示，主要包括頻段內(nèi)圖像配準(zhǔn)、去平地相位、干涉參數(shù)精確估計(jì)、頻段間干涉相位配準(zhǔn)、雙頻聯(lián)合相位解纏、解纏相位噪點(diǎn)剔除等步驟。下面通過(guò)比較分析，將給出各個(gè)步驟的處理方法，并在機(jī)載實(shí)測(cè)雙頻干涉處理中驗(yàn)證其有效性。

2.1 圖像配準(zhǔn)

圖像配準(zhǔn)是通過(guò)調(diào)整主、輔SAR復(fù)圖像，使得兩幅圖像上相同的點(diǎn)對(duì)應(yīng)地面場(chǎng)景上的同一點(diǎn)。配準(zhǔn)的步驟如下[20,21]：

步驟1 主輔圖像選取若干控制點(diǎn)，采用相關(guān)法計(jì)算控制點(diǎn)處粗偏移量，精度為1個(gè)像素。

計(jì)算控制點(diǎn)處的偏移量可以采用空域相關(guān)函數(shù)法的計(jì)算基于窗口的配準(zhǔn)方法，分別在主圖像s1和輔圖像s2以控制點(diǎn)為中心建立匹配窗和搜索窗。滑動(dòng)搜索窗求得窗口的相關(guān)函數(shù)值，相關(guān)函數(shù)最大的滑動(dòng)位置即為圖像的粗偏移量(u，v)；

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

步驟3 在求得控制點(diǎn)處精確偏移量后，可以利用二項(xiàng)式進(jìn)行擬合：

其中， (x，y)為主圖像的控制點(diǎn)坐標(biāo)，通過(guò)最小二乘法可以得擬合系數(shù)，再利用擬合系數(shù)得到輔圖像整幅圖像配準(zhǔn)后精確的位置，通過(guò)線(xiàn)性插值把輔圖像插值到求得的精確位置來(lái)實(shí)現(xiàn)精配準(zhǔn)[21]。

下面采用某機(jī)載雙頻干涉系統(tǒng)飛行的數(shù)據(jù)為例來(lái)分析圖像配準(zhǔn)技術(shù)。

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

2.2 去平地

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

2.3 干涉參數(shù)精確估計(jì)

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

2.4 頻段間干涉圖配準(zhǔn)

交叉頻段干涉相位配準(zhǔn)精度往往較低，本文采用的方法為首先進(jìn)行同頻段圖像配準(zhǔn)、濾波、去平地，通過(guò)利用各自主圖像的幅度信息擬合出多項(xiàng)式系數(shù)，進(jìn)而來(lái)實(shí)現(xiàn)頻段間干涉相位的配準(zhǔn)。

由于不同頻段的主圖像幅度存在較多的斑點(diǎn)噪聲，造成不同頻段的主圖像幅度圖的相似性較差。由圖5可知，NL-InSAR可以有效地對(duì)SAR圖像的幅度圖進(jìn)行去噪。圖7(a)為不進(jìn)行NL-InSAR幅度濾波的C波段主圖像幅度圖，圖7(b)為不進(jìn)行NLInSAR濾波的X波段主圖像幅度?？梢钥闯霾贿M(jìn)行NL-InSAR濾波的不同波段主圖像幅度圖較多的斑點(diǎn)噪聲，造成不同波段的主圖像幅度圖之間相似性較差。圖8(a)為采用NL-InSAR對(duì)幅度濾波后的C波段主圖像幅度，圖8(b)為采用NL-InSAR對(duì)幅度濾波后的X波段主圖像。NL-InSAR可以很好地去除SAR圖像幅度圖的斑點(diǎn)噪聲，同時(shí)較好地保持了圖像的邊緣細(xì)節(jié)信息，提高了不同波段主圖像幅度圖之間的相似性。因此本文采用經(jīng)過(guò)NL-InSAR濾波后的兩個(gè)頻段的主圖像幅度擬合出式(1)的多項(xiàng)式系數(shù)，從而實(shí)現(xiàn)對(duì)兩個(gè)頻段的干涉相位圖之間的配準(zhǔn)。針對(duì)實(shí)測(cè)數(shù)據(jù)不同頻段配準(zhǔn)前后的不同波段之間相關(guān)系數(shù)圖如圖9所示。圖9(a)所示的配準(zhǔn)前X波段和C波段的相關(guān)系數(shù)較低，均值為0.46，圖9(b)所示的不采用NL-InSAR對(duì)幅度進(jìn)行濾波的配準(zhǔn)后的相關(guān)系數(shù)相對(duì)較高，均值為0.90。圖9(c)所示為采用NL-InSAR對(duì)幅度進(jìn)行濾波的配準(zhǔn)好的相關(guān)系數(shù)最高，均值為0.94，驗(yàn)證了本文提出的不同波段干涉相位配準(zhǔn)方法的有效性。

2.5 多通道聯(lián)合相位解纏

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

假設(shè)所有像素點(diǎn)之間是相互獨(dú)立的，那么：

對(duì)于N0個(gè)波段相位聯(lián)合解纏情況下，相應(yīng)似然函數(shù)為：

其中單個(gè)波段的觀(guān)測(cè)相位的似然函數(shù)為：

我們可以使用馬爾科夫隨機(jī)場(chǎng)來(lái)描述圖像的局部特性，根據(jù)MarKov-Gibbs等價(jià)性可以用Gibbsian分布來(lái)描述馬爾科夫先驗(yàn)分布[8]：

變分能量項(xiàng)最大后驗(yàn)估計(jì)使用的是一種在圖像領(lǐng)域常用的整體變分(TV)模型，其能量函數(shù)定義為：

式(7)中，β是一個(gè)固定的常量，用來(lái)調(diào)節(jié)真實(shí)相位的光滑性與觀(guān)測(cè)相位和真實(shí)相位之間相似性的平衡的常數(shù)，Np為像素點(diǎn)p的鄰域的集合，wp，q表征鄰域連通性，若Np為4鄰域模型，其值為1。

把式(7)代入式(6)結(jié)合式(2)–式(5)并取對(duì)數(shù)且添加負(fù)號(hào)，可以得到變分能量項(xiàng)最大后驗(yàn)估計(jì)的在TV模型下的相位估計(jì)公式為：

上述優(yōu)化模型可以通過(guò)圖割理論求解[11]，通過(guò)求解上述優(yōu)化模型可以求解得到各個(gè)波段干涉相位的纏繞量和聯(lián)合解纏繞相位。

對(duì)添加相同噪聲的仿真數(shù)據(jù)，采用上述3種方法分別進(jìn)行聯(lián)合解纏結(jié)果如圖10所示，用來(lái)仿真干涉相位的原始相位相減的殘差圖如圖11所示。通過(guò)上述仿真分析，進(jìn)一步驗(yàn)證了變分能量項(xiàng)最大后驗(yàn)估計(jì)用于雙頻聯(lián)合解纏的魯棒性。

2.6 聚類(lèi)均值濾波

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

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

如圖12所示，將兩種標(biāo)記結(jié)果相與，得到用于加權(quán)均值濾波的標(biāo)記結(jié)果，對(duì)聯(lián)合解纏結(jié)果進(jìn)行如下加權(quán)均值濾波：

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

下面采用實(shí)際飛行數(shù)據(jù)對(duì)本文方法進(jìn)行驗(yàn)證，同時(shí)對(duì)比了單波段和雙波段的高程重建效果。實(shí)測(cè)數(shù)據(jù)是由中國(guó)科學(xué)院電子學(xué)研究所航天微波遙感系統(tǒng)部提供的某場(chǎng)景雙波段(C波段中心頻率5.4 GHz和X波段中心頻率9.6 GHz)的機(jī)載雙頻干涉數(shù)據(jù)。

實(shí)測(cè)數(shù)據(jù)包括同一場(chǎng)景的C波段和X波段主、輔SAR圖像復(fù)數(shù)據(jù)，處理數(shù)據(jù)的場(chǎng)景如圖16所示，為陜西澄城地區(qū)的某山區(qū)，其中地形起伏比較大區(qū)域?yàn)閳D中右上角的高架，不同波段配準(zhǔn)之后的C波段干涉圖和X波段干涉圖分別如圖17所示。

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

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

由于本文所用的實(shí)測(cè)數(shù)據(jù)場(chǎng)景中沒(méi)有準(zhǔn)確的地本文采用的方法是將谷歌光學(xué)地圖面控制點(diǎn)，本文采用的方法是將谷歌光學(xué)地圖作為主圖像和SAR圖像的幅度圖進(jìn)行配準(zhǔn)，得到配準(zhǔn)的多項(xiàng)式系數(shù)，利用擬合的多項(xiàng)式系數(shù)實(shí)現(xiàn)對(duì)DEM進(jìn)行配準(zhǔn)。在光學(xué)圖像上選取5個(gè)控制點(diǎn)對(duì)DEM進(jìn)行精度分析。單波段反演高程的均方根誤差(Root Mean Square Error, RMSE)為4.5826 m，且存在部分地區(qū)高程反演失敗，而雙頻聯(lián)合反演的RMSE為1.3476 m，可以看出本文提出雙頻聯(lián)合高程重建方法精度優(yōu)于單波段高程重建。

4 結(jié)束語(yǔ)

本文提出了一種有效的雙頻干涉處理方法，分析并給出了算法中各個(gè)步驟的處理方案。針對(duì)雙頻聯(lián)合處理中的雙頻聯(lián)合配準(zhǔn)、雙頻解纏繞等關(guān)鍵技術(shù)難題，給出了有效的解決方案，提高了雙頻聯(lián)合反演DEM的精度和魯棒性。采用常規(guī)單波段方法和本文方法對(duì)機(jī)載實(shí)際飛行數(shù)據(jù)處理，處理結(jié)果表明，單波段使用MCF進(jìn)行相位解纏，無(wú)法實(shí)現(xiàn)大起伏地形區(qū)域的高程重建；而采用本文提出的處理流程能夠準(zhǔn)確地反演大起伏地形區(qū)域的高程，驗(yàn)證了本文所提出的雙頻干涉處理方法的有效性。