礦區(qū)沉降監(jiān)測的InSAR質(zhì)量圖相位解纏方法

劉萬利 張秋昭 胡 江

（1.徐州工程學(xué)院數(shù)學(xué)與統(tǒng)計學(xué)院，江蘇 徐州 221018；2.中國礦業(yè)大學(xué)環(huán)境與測繪學(xué)院，江蘇 徐州 221116）

二維相位解纏是合成孔徑雷達干涉測量技術(shù)（Interferometric Synthetic Aperture Radar，InSAR）中的一個關(guān)鍵環(huán)節(jié)，解纏的成功與否直接影響到目標(biāo)高程信息的準確性，因此一直是InSAR技術(shù)研究的熱點和難點［1-2］。近年來，許多學(xué)者對相位解纏進行了深入研究，主要可分為兩大類：一類是基于積分的路徑跟蹤算法，另一類是最小范數(shù)算法［3］。質(zhì)量圖在這兩類算法中起著非常重要的作用，直接影響其解纏精度［4］。尤其對于路徑跟蹤算法，一個好的質(zhì)量指標(biāo)通?？梢员苊饣驕p少積分路徑穿過殘差點，可以使得高質(zhì)量區(qū)域獲得比較高的精度［5］。

基于此，許多學(xué)者就如何選擇合適的質(zhì)量指標(biāo)展開了相關(guān)研究［6-14］。GHIGLIA和PRITT在文獻［12］中比較了相干系數(shù)（Correlation Coefficient，CC）、偽相干系數(shù)（Pseudo Correlation Coefficient，PCC）、相位導(dǎo)數(shù)變化（Phase Derivative Variation，PDV）、最大相位梯度（Maximum Phase Gradient，MPG）4種質(zhì)量指標(biāo)，并指出相干系數(shù)是最有效的一個。OSMANOGLU等［13］描述了6種質(zhì)量指標(biāo)，試驗結(jié)果顯示Fisher距離的性能總體上優(yōu)于相干系數(shù)、相位導(dǎo)數(shù)變化、結(jié)合枝切法的相位導(dǎo)數(shù)變化、二階導(dǎo)數(shù)可靠度及線掃描指標(biāo)等5種質(zhì)量指標(biāo)。LIU等［15］提出一種基于灰度共生矩陣（Gray Level Co-occurrence Matrix，GLCM）的熵差質(zhì)量指標(biāo)，該指標(biāo)對分級級數(shù)要求比較高，需要事先確定級數(shù)值。劉萬利等［5］對幾種常見的質(zhì)量指標(biāo)做了比較深入的分析研究，并提出了一種改進的Fisher距離指標(biāo)，試驗結(jié)果表明，F(xiàn)isher和改進的Fisher信息質(zhì)量指標(biāo)總體上比較穩(wěn)定。

綜上分析，基于相干系數(shù)和相位導(dǎo)數(shù)變化的質(zhì)量指標(biāo)簡單易行，相干系數(shù)強調(diào)時間相干性，相位導(dǎo)數(shù)變化強調(diào)空間相似性［5］；并且上述質(zhì)量圖指標(biāo)大多單方面偏重時間相關(guān)、空間相關(guān)、噪聲等某一個指標(biāo)。InSAR技術(shù)在監(jiān)測山區(qū)礦區(qū)地表形變時，獲得的干涉圖通常具有噪聲大、條紋密且復(fù)雜等特點，單一質(zhì)量圖指標(biāo)用于山區(qū)礦區(qū)等復(fù)雜區(qū)域干涉圖相位解纏時，適應(yīng)性較差。為此，本研究提出一種融合相位導(dǎo)數(shù)變化和相關(guān)系數(shù)兩種質(zhì)量指標(biāo)的新的質(zhì)量圖指標(biāo)。該指標(biāo)包含了空間相關(guān)性、時間相關(guān)性等更多的質(zhì)量信息，因此更適用于地物覆蓋復(fù)雜地區(qū)、條紋比較密集干涉圖的相位解纏（尤其是山區(qū)煤礦區(qū)沉降監(jiān)測干涉圖）。采用模擬數(shù)據(jù)和實測數(shù)據(jù)試驗驗證了所提質(zhì)量指標(biāo)的可靠性和適用性。

1 相位解纏理論

基于濾波的相位解纏方法在展開密條紋高噪聲干涉圖相位時具有獨特的優(yōu)勢［5，15］，本研究采用容積卡爾曼濾波相位解纏方法［16］進行分析。

1.1 相位解纏模型

容積卡爾曼濾波相位解纏模型在繼承卡爾曼濾波相位解纏模型的基礎(chǔ)上，采用容積卡爾曼濾波器處理非線性方程，其狀態(tài)方程和觀測方程為

式中，y(k)和x(k)分別為像元k的觀測值和解纏相位；φ(k)為相位梯度估計值，通過最大似然估計器獲得；h[?]表示觀測函數(shù)；v(k)為k像元的觀測誤差向量；w(k)為梯度估計誤差，被認為高斯白噪聲且滿足如下條件：

式中，E[?]表示期望；D[?]表示方差；為梯度估計誤差方差，其計算參考文獻［5］。v1(k)和v2(k)分別為復(fù)觀測值實部和虛部的觀測誤差，被認為是高斯白噪聲且滿足如下條件：

式中，R(k)為觀測誤差協(xié)方差陣；，SkSNR為像元k的信噪比。

容積卡爾曼濾波相位解纏是沿著某一特定路徑進行的，而路徑的選擇是基于某一質(zhì)量指標(biāo)從高質(zhì)量像元到低質(zhì)量像元。質(zhì)量指標(biāo)選擇得當(dāng)與否直接關(guān)系到解纏效果。

1.2 相位解纏步驟

本研究相位解纏實施流程如圖1所示。具體步驟為：①根據(jù)某一質(zhì)量指標(biāo)生成質(zhì)量圖；②選取質(zhì)量圖中的一個高質(zhì)量點A1作為解纏始點；③檢測像素點A1的4個鄰接像素點；④利用容積卡爾曼濾波解纏模型估計4個鄰接像素點中最高質(zhì)量的像素點A2；⑤選取像素點A1、A2鄰接的質(zhì)量最高的未解纏像素點A3并解纏。重復(fù)以上步驟，直到相位圖中的全部像素點都被解纏完畢。

2 基于相干系數(shù)與相位導(dǎo)數(shù)變化的新質(zhì)量指標(biāo)

相干系數(shù)和相位導(dǎo)數(shù)變化是兩種較常用的衡量干涉圖像質(zhì)量的質(zhì)量指標(biāo)，分別介紹如下。

2.1 相干系數(shù)

相干系數(shù)是最直接、最常用的衡量［12］指標(biāo)，計算公式為

式中，s1=c1+n1，s2=c2+n2，c1、c2為圖像的信號部分，n1、n2為圖像的噪聲部分；[?]*為共軛運算。

對于實測干涉圖，CC是綜合多種失相干因素得出的結(jié)果，表征的是兩幅參與干涉的影像的相干性。通常情況下，反射能力強的地物對應(yīng)像元的CC值較高。本質(zhì)上是用相位的一階導(dǎo)數(shù)來描述相位質(zhì)量，對相位梯度比較敏感。因此，對于地形陡峭的區(qū)域，通過計算相位梯度，會將其標(biāo)識為低質(zhì)量數(shù)據(jù)。

2.2 相位導(dǎo)數(shù)變化

對于二維干涉圖數(shù)據(jù)，某像元的PDV值計算公式為［12］

由式（5）可以看出，PDV是以當(dāng)前像元為中心的一定窗口內(nèi)的像元在水平方向和豎直方向上梯度方差大小之和，反映的是目標(biāo)窗口內(nèi)梯度數(shù)值的離散程度，是梯度空間相似性好壞的評價指標(biāo)。本質(zhì)上是用相位的二階導(dǎo)數(shù)來描述相位質(zhì)量，對梯度變化比較敏感。因此在一些相干性較低的同質(zhì)地物覆蓋區(qū)，該指標(biāo)會表現(xiàn)出高質(zhì)量結(jié)果。

2.3 PDV和CC融合質(zhì)量指標(biāo)

由2.1及2.2節(jié)可以看出CC可以很好地反映像元的個體特征，但不能很好地反映空間相似性，對于一些存在孤立高相干點的干涉圖，容易造成解纏路徑穿過低質(zhì)量區(qū)，去解纏那些高質(zhì)量點然后再穿過低質(zhì)量區(qū)去尋找次高質(zhì)量點，從而造成誤差過早累積。而PDV會將一些空間上勻質(zhì)但相干性不高的區(qū)域誤判為高質(zhì)量區(qū)?；诖耍狙芯刻岢隽艘环N融合PDV和CC的質(zhì)量指標(biāo)。運算步驟如下：

（1）數(shù)據(jù)標(biāo)準化。由式（4）和式（5）獲得的CC和PDV值量綱和數(shù)量級不同，首先對其進行標(biāo)準化處理可得：

由式（7）可知，PDV_CC值越小，質(zhì)量越高，該指標(biāo)融合了像元相位的空間相似性和時間相干性兩個特征，是一種信息比較全面的質(zhì)量指標(biāo)。

3 試驗與分析

我國多數(shù)煤礦區(qū)地表下沉速度快、植被覆蓋多，尤其是山區(qū)礦區(qū)地形起伏大，因此InSAR干涉圖通常呈現(xiàn)出條紋密集復(fù)雜、噪聲大、殘差點多等特點。如果質(zhì)量指標(biāo)選用不當(dāng)，會造成解纏路徑穿過大量殘差點，進而導(dǎo)致沉降監(jiān)測精度不理想。本研究采用一組多山場景模擬數(shù)據(jù)（干涉圖具有條紋密集、噪聲大的特點）和礦區(qū)實測數(shù)據(jù)對所提出的質(zhì)量指標(biāo)進行試驗分析。為詳細比較試驗結(jié)果，整個相位解纏過程只有質(zhì)量指標(biāo)不同，解纏算法采用容積卡爾曼濾波相位解纏方法。采用解纏誤差d和偏差χ2兩個指標(biāo)進行定量評價，兩者計算公式為

式中，φu為解纏相位值；φr為相位先驗值；σ為標(biāo)準差；N為樣本數(shù)。

3.1 模擬數(shù)據(jù)試驗

多山場景地形模擬數(shù)據(jù)如圖2（a）所示。為了詳細分析地形起伏對相位解纏的影響，模擬數(shù)據(jù)的垂直基線設(shè)為150 m，像元個數(shù)為250×250。數(shù)據(jù)具有地勢起伏大、陡峭度較高等特點。此場景的模擬成像加噪干涉圖如圖2（c）所示。本試驗的含噪干涉圖體現(xiàn)的是幾何失相干程度，幾何失相干嚴重的區(qū)域噪聲強度大，由干涉圖可以看出多個區(qū)域條紋比較密集。

圖3給出了模擬數(shù)據(jù)的3種質(zhì)量圖。從大致的紋理來看，CC質(zhì)量圖與其余兩種質(zhì)量圖差別較大，高低質(zhì)量區(qū)區(qū)分比較明顯。這是因為模擬數(shù)據(jù)CC值的計算僅體現(xiàn)了幾何失相干，條紋密集的區(qū)域質(zhì)量低，稀疏的區(qū)域質(zhì)量高。PDV_CC與PDV的大致紋理比較接近，但顏色有所區(qū)別，前者顏色介于CC和PDV之間。這是由于PDV_CC是融合了CC和PDV兩種質(zhì)量指標(biāo)所致。另外，由圖3還可以看出在條紋稀疏區(qū)，PDV_CC和PDV也有少量低質(zhì)量散點，是因為含噪干涉圖的噪聲添加具有隨機特性，而這兩者的計算都是基于含噪干涉相位的。

圖4給出了3種質(zhì)量指標(biāo)引導(dǎo)解纏路徑的相位解纏結(jié)果圖。相對于真實相位圖（圖2（b）），可以看出CC引導(dǎo)的解纏失?。▓D4（a）），PDV（圖4（b））和PDV_CC（圖4（c））解纏結(jié)果總體保持了原樣，但在條紋非常密集的區(qū)域解纏結(jié)果表現(xiàn)出了不連續(xù)現(xiàn)象，即解纏相位誤差較大。這是因為這些區(qū)域噪聲強度大、圖像質(zhì)量低、條紋密集且變化較快，進而導(dǎo)致條紋梯度的估計誤差較大。進一步觀察可以發(fā)現(xiàn)，在條紋密集的區(qū)域，PDV_CC較PDV解纏結(jié)果更連續(xù)光滑（主要體現(xiàn)在橢圓和矩形區(qū)域）。說明兼顧空間相似性和時間相干性的PDV_CC質(zhì)量指標(biāo)優(yōu)于僅考慮空間相似性的PDV質(zhì)量指標(biāo)。

為了進一步討論所提質(zhì)量指標(biāo)的可靠性，圖5給出了3種指標(biāo)的解纏誤差圖及其對應(yīng)的誤差直方圖。從誤差圖可以明顯看出CC解纏誤差圖右側(cè)部分出現(xiàn)了大面積的大誤差區(qū)域，PDV和PDV_CC解纏誤差圖只是在條紋密集噪聲較大的區(qū)域存在較大誤差，這說明CC指標(biāo)最不可靠。對比圖5（b）和圖5（c），可以看出PDV_CC方法解纏圖中橢圓和矩形區(qū)域的不連續(xù)區(qū)域明顯小于PDV方法，這說明PDV_CC質(zhì)量指標(biāo)優(yōu)于PDV。從誤差直方圖也可以看出PDV_CC的誤差分布圖形中心峰（落在0附近部分）高度最高且寬度最窄，說明該指標(biāo)引導(dǎo)的解纏結(jié)果誤差方差最小。因此PDV_CC的解纏結(jié)果優(yōu)于其他兩種指標(biāo)。

表1給出了模擬數(shù)據(jù)3種質(zhì)量指標(biāo)性能的定量評價指標(biāo)取值，分別為誤差均值、誤差方差以及偏差均值。從解纏誤差和偏差表達式可知，以上3個定量評價指標(biāo)的值越接近0說明解纏結(jié)果越好。因此從表1可以看出CC指標(biāo)最不可靠，PDV_CC優(yōu)于CC，是一種比較魯棒可靠的質(zhì)量指標(biāo)。

分析偏差值隨解纏路徑的變化情況可以很好地解釋誤差從哪里產(chǎn)生、誤差隨路徑如何傳播等問題。一個好的質(zhì)量指標(biāo)應(yīng)該是較大的偏差值出現(xiàn)得越晚越好，這樣可以避免誤差較早累積。模擬數(shù)據(jù)解纏結(jié)果的偏差值隨解纏路徑的走向曲線圖如圖6所示。可以看出，CC的偏差值最大，PDV次之，PDV_CC最??；CC方法的大誤差出現(xiàn)最早（在3 700個像元左右），PDV次之，PDV_CC方法的大誤差出現(xiàn)最晚（在56 000個像元左右）。這說明PDV_CC方法出現(xiàn)誤差累積現(xiàn)象較晚，有更多的像元解纏結(jié)果較可靠。

3.2 礦區(qū)實測數(shù)據(jù)試驗

研究區(qū)為神東礦區(qū)大柳塔煤礦52304工作面，像元個數(shù)為300×250，區(qū)域內(nèi)地形略有緩坡，有荒草沙石覆蓋。植被類型有低矮灌木、沙生植被、荒草等。干涉圖由2013年4月2日及2013年4月24日兩景TerraSAR-X 影像生成（圖7（a）），外部DEM采用課題組通過收集的航拍地形圖資料并結(jié)合GPS控制點生成的Relief-DEM。干涉圖的生成采用Gamma軟件。由于地表覆蓋導(dǎo)致的失相干較嚴重，故干涉圖采用的是2視處理結(jié)果。3種質(zhì)量指標(biāo)對應(yīng)的質(zhì)量圖見圖7（b）至圖7（d）。由于試驗采用的是沉降監(jiān)測干涉圖，整個區(qū)域的真實沉降值難以獲得，因此解纏結(jié)果主要從連續(xù)性角度進行分析。

由圖7（b）至圖7（d）可以看出，CC質(zhì)量圖的高質(zhì)量區(qū)比較小，而PDV和PDV_CC顯現(xiàn)出大面積高質(zhì)量區(qū)，這是因為后兩個指標(biāo)考慮了空間相似性，而研究區(qū)地面覆蓋空間相似性較高。PDV_CC與PDV質(zhì)量圖相似度較高，但比后者顏色稍淡，與公式所體現(xiàn)的特征相符。由圖7（e）至圖7（g）解纏相位圖可以看出，CC質(zhì)量指標(biāo)解纏結(jié)果在形變區(qū)右下角出現(xiàn)了較大的相位突變，相位連續(xù)性最差。在此區(qū)域，PDV和PDV_CC兩種方法解纏效果相當(dāng)，均優(yōu)于CC解纏方法。結(jié)合干涉圖和質(zhì)量圖可以看出，此區(qū)域噪聲大并且失相干比較嚴重，說明在低質(zhì)量區(qū)PDV和PDV_CC指標(biāo)魯棒性較好。在形變區(qū)右上方，信噪比較高，CC和PDV_CC取得了比較連續(xù)的解纏結(jié)果，而PDV方法出現(xiàn)了相位跳變。這說明PDV指標(biāo)在引導(dǎo)解纏路徑時出現(xiàn)了從低質(zhì)量到高質(zhì)量的情況，導(dǎo)致誤差提前大量累積從而造成高質(zhì)量區(qū)解纏結(jié)果不理想。從3種質(zhì)量指標(biāo)引導(dǎo)的解纏順序圖（圖7（h）至圖7（j），尤其是圖7（j）橢圓區(qū)域）可以看出，解纏結(jié)果的不同是由于質(zhì)量指標(biāo)引導(dǎo)的路徑不同導(dǎo)致的誤差累積不同所致。由此試驗可以看出本研究提出的PDV_CC質(zhì)量指標(biāo)魯棒性較好。

綜合以上試驗結(jié)果可以看出，無論是模擬數(shù)據(jù)還是實測數(shù)據(jù)，PDV_CC都取得了優(yōu)于PDV和CC的解纏結(jié)果。這說明融合指標(biāo)可靠性好，是一種比較適合引導(dǎo)復(fù)雜地物覆蓋、沉降較快工作面干涉圖相位解纏的指標(biāo)。此外，對于模擬數(shù)據(jù)，CC指標(biāo)解纏結(jié)果極其不可靠，是因為模擬的地形在小范圍內(nèi)起伏變化較快，存在高質(zhì)量孤島；CC指標(biāo)引導(dǎo)的路徑出現(xiàn)了從低質(zhì)量到高質(zhì)量情形，從而造成了誤差累積傳播。對于實測數(shù)據(jù)，在沉陷區(qū)PDV_CC和PDV取得了相似的解纏結(jié)果，較符合實際沉降情況，但在非沉陷區(qū)PDV出現(xiàn)了解纏相位不連續(xù)現(xiàn)象，而CC和PDV_CC結(jié)果相似且較符合實際情況。兩種試驗結(jié)果有所差別，主要是因為模擬數(shù)據(jù)的噪聲產(chǎn)生原因比較單一，主要由幾何失相干引起，而實測數(shù)據(jù)噪聲原因比較復(fù)雜。

4 結(jié)論

（1）提出了一種融合PDV和CC的新質(zhì)量指標(biāo)，該指標(biāo)不僅考慮了時間相干性還兼顧了空間相似性，更適用于噪聲大、條紋密集的復(fù)雜干涉圖相位解纏。

（2）高山地形模擬數(shù)據(jù)和礦區(qū)實測數(shù)據(jù)試驗均表明，新的質(zhì)量指標(biāo)與單一質(zhì)量指標(biāo)相比，具有較好的魯棒性和實用性。

（3）相位梯度估計噪聲與質(zhì)量指標(biāo)的相互影響機制是下一步研究的方向。