王志勇，孟慶穎

(山東科技大學(xué) 測(cè)繪科學(xué)與工程學(xué)院，山東 青島 266590)

1 引 言

隨著地下礦產(chǎn)資源的開采，礦區(qū)地表沉陷問題日益嚴(yán)重，在礦區(qū)形成大量的塌陷坑、塌陷槽，引起農(nóng)田破壞、房屋損毀等，嚴(yán)重影響人民的生命財(cái)產(chǎn)安全和經(jīng)濟(jì)的可持續(xù)發(fā)展。雷達(dá)干涉測(cè)量(Interferometric Sythetic Aperture Radar，InSAR)技術(shù)提供了一種全新的地表形變監(jiān)測(cè)方法，具有全天候、全天時(shí)、低成本等優(yōu)勢(shì)[1-6]。相位解纏是InSAR數(shù)據(jù)處理及應(yīng)用的一個(gè)關(guān)鍵技術(shù)和處理步驟，特別是在InSAR礦區(qū)沉降監(jiān)測(cè)數(shù)據(jù)處理中，它直接影響InSAR地表形變監(jiān)測(cè)的可靠性及精確性。

InSAR技術(shù)已經(jīng)初步被應(yīng)用于礦區(qū)地面沉降監(jiān)測(cè)中[7]，但標(biāo)志性的成果還很少，尚未在礦區(qū)沉降監(jiān)測(cè)中得到大規(guī)模應(yīng)用，這主要是由于礦區(qū)環(huán)境的復(fù)雜性以及InSAR技術(shù)的局限性決定的。礦區(qū)形變比較特殊[7]，主要表現(xiàn)為沉降速率不均勻、礦區(qū)形變量非常大(如在活躍期內(nèi)，一個(gè)月的沉降量可達(dá)半米，甚至更大)、沉降范圍分散；其次，沉降范圍很小(多為半徑幾百米的沉降漏斗)；礦區(qū)農(nóng)田多、地表植被覆蓋率高，雷達(dá)數(shù)據(jù)相干性很差，致使在干涉條紋圖中出現(xiàn)很多不連續(xù)的區(qū)域?；谶@些特點(diǎn)，部分InSAR相位解纏在礦區(qū)沉降監(jiān)測(cè)中不能正確解纏，致使InSAR監(jiān)測(cè)的可靠性存在很大問題。

由于不同的相位解纏算法都具有一定的適應(yīng)性，并不是所有的相位解纏算法都能解決礦區(qū)大形變沉降監(jiān)測(cè)時(shí)的相位解纏問題。因此，針對(duì)礦區(qū)形變的特殊性，對(duì)現(xiàn)有的相位解纏算法進(jìn)行對(duì)比分析，探討適合礦區(qū)沉降特點(diǎn)的InSAR相位解纏算法，對(duì)于礦區(qū)沉降監(jiān)測(cè)的可靠性及精度方面都具有十分重要的意義。

本文選取沉降特點(diǎn)非常明顯的濟(jì)寧市梁寶寺礦區(qū)為試驗(yàn)區(qū)域，采用真實(shí)的ALOS PALSAR(日本雷達(dá)衛(wèi)星，2006年發(fā)射)雷達(dá)數(shù)據(jù)，生成差分干涉圖。分別用經(jīng)典的相位解纏方法進(jìn)行礦區(qū)大形變區(qū)域的相位解纏實(shí)驗(yàn)，并對(duì)解纏結(jié)果進(jìn)行分析評(píng)價(jià)，找出適于礦區(qū)環(huán)境的、精度較高、可靠性較好的解纏方法。

2 相位解纏算法分析

所謂的相位解纏就是將相位由主值恢復(fù)為真實(shí)值的過程[2]。在InSAR數(shù)據(jù)處理中，針對(duì)相位解纏已經(jīng)進(jìn)行了大量的研究，提出了許多算法，主要可以概括為以下三大類[3，8-9，14]：路徑跟蹤類相位解纏算法、最小范數(shù)類相位解纏算法、基于最優(yōu)估計(jì)的相位解纏算法。

路徑跟蹤類算法主要采用不同的策略尋找最優(yōu)化的積分路徑，避開殘差點(diǎn)造成的誤差傳遞，以滿足相位梯度閉合路徑積分為零的條件，它一般是通過識(shí)別殘差點(diǎn)，設(shè)置正確的枝切線阻止積分路徑穿過；或者是在相位質(zhì)量圖的幫助下，從高質(zhì)量數(shù)據(jù)開始積分。主要包括：枝切法[10]、掩膜切線法[9]、Flynn的最小不連續(xù)法等。

最小范數(shù)類相位解纏算法是一種全局求解算法，以纏繞相位的離散偏微分與解纏相位的離散偏微分的差最小為準(zhǔn)則建立全局的數(shù)學(xué)模型來求解相位解纏的估計(jì)值。主要包括：最小LP范數(shù)法[2-3]、高斯-賽德爾迭代法、PCG預(yù)解共軛梯度法、FFT/DCT最小二乘法、多級(jí)格網(wǎng)法[2-3]。

基于最優(yōu)估計(jì)的相位解纏算法將最優(yōu)估計(jì)的算法應(yīng)用到相位解纏中，主要包括：網(wǎng)絡(luò)規(guī)劃法[11]、Kalman濾波法[12]、遺傳算法[13]等。

表1列出了幾種經(jīng)典的相位解纏算法的核心思想以及其優(yōu)缺點(diǎn)。

表1 相位解纏方法及特點(diǎn)分析

3 研究區(qū)域及InSAR相位解纏實(shí)驗(yàn)

3.1 研究區(qū)域及實(shí)驗(yàn)數(shù)據(jù)

選取了我國(guó)礦區(qū)地面沉降比較嚴(yán)重的濟(jì)寧礦區(qū)的梁寶寺煤礦作為實(shí)驗(yàn)區(qū)，選用2009年01月10日和2009年02月25日獲取的兩景ALOS PALSAR數(shù)據(jù)組成干涉對(duì)，其path號(hào)為449，F(xiàn)rame號(hào)為700，其工作模式為FBS(Fine Beam Single Polarization)模式，L波段(波長(zhǎng)23.6cm)，產(chǎn)品級(jí)別為L(zhǎng)evel1.1(即單視復(fù)影像數(shù)據(jù))，HH極化，視角為34.3°，方位向像元大小為3.148m，距離向像元大小為4.684m。干涉對(duì)的垂直基線為268.542m，時(shí)間間隔46d。

除此之外，為了去除地形相位的影響，本文采用了SRTM DEM數(shù)據(jù)。

3.2 InSAR技術(shù)生成差分干涉圖

采用雙軌法差分干涉測(cè)量[7]的方法生成差分干涉圖，并采用改進(jìn)的Goldstein濾波算法濾除了部分相位噪聲。為使實(shí)驗(yàn)更具針對(duì)性，對(duì)原始單視復(fù)影像(SLC)進(jìn)行區(qū)域裁剪，獲得單個(gè)礦區(qū)的差分干涉圖，采用1∶2(距離向：方位向)多視處理，生成的差分干涉圖的大小為：470行×481列。圖1(a)為梁寶寺礦區(qū)的濾波后的差分干涉圖，圖1(b)為相干圖。

圖1 差分干涉圖及相干圖

3.3 相位解纏對(duì)比實(shí)驗(yàn)

對(duì)濾波后的差分干涉相位圖分別采用枝切法、區(qū)域生長(zhǎng)法(Region Grow)、質(zhì)量引導(dǎo)掩膜切線法、最小不連續(xù)法、LP最小范數(shù)法、無權(quán)多級(jí)網(wǎng)絡(luò)法、加權(quán)多級(jí)網(wǎng)絡(luò)法、PCG預(yù)解共軛梯度法、最小費(fèi)用流法(minimum cost flow algorithm，MCF)進(jìn)行解纏實(shí)驗(yàn)，解纏結(jié)果如圖2所示。

圖2 不同相位解纏方法解纏結(jié)果

4 相位解纏結(jié)果分析及評(píng)價(jià)

為了對(duì)比各種相位解纏算法在礦區(qū)沉降監(jiān)測(cè)中的解纏結(jié)果，從剖面圖、統(tǒng)計(jì)分析、運(yùn)行時(shí)間、形變量等幾個(gè)方面進(jìn)行對(duì)比分析。

4.1 直接分析

根據(jù)圖2，分析不同相位解纏方法得到的解纏結(jié)果發(fā)現(xiàn)，基于路徑跟蹤的相位解纏算法整體效果較差，只有Flynn最小不連續(xù)法解纏效果相對(duì)較好，這是因?yàn)橹η蟹ㄔ谙辔徊贿B續(xù)點(diǎn)較多且比較密集的時(shí)候，無法進(jìn)行解纏，區(qū)域生長(zhǎng)法需要從質(zhì)量較好的區(qū)域開始，質(zhì)量引導(dǎo)掩膜切線需要提供質(zhì)量較好的品質(zhì)圖像，礦區(qū)環(huán)境復(fù)雜，相干性較差，相位噪聲嚴(yán)重，很難滿足上述要求，只有Flynn最小不連續(xù)法可以在沒有相位質(zhì)量圖指導(dǎo)的情況下較好地完成相位解纏；基于最小范數(shù)的相位解纏算法總體解纏效果不錯(cuò)，但是加權(quán)多級(jí)網(wǎng)絡(luò)由于權(quán)重選擇導(dǎo)致其在噪聲較大的區(qū)域無法實(shí)現(xiàn)正確解纏；基于網(wǎng)絡(luò)規(guī)劃的相位解纏算法，最小費(fèi)用流法(MCF)出現(xiàn)了較多不可解纏的區(qū)域，解纏效果較差。

4.2 剖面圖

實(shí)驗(yàn)選取了其中的一行數(shù)據(jù)(位于差分干涉圖的第210行)，對(duì)比分析每種相位解纏算法的結(jié)果圖的剖面圖，所生成的剖面圖如圖3所示。

根據(jù)圖3，通過分析不同解纏方法下的剖面圖，發(fā)現(xiàn)枝切法、最小不連續(xù)法 、LP最小范數(shù)、加權(quán)多級(jí)網(wǎng)絡(luò)剖面線走勢(shì)大致相同；區(qū)域生長(zhǎng)法、最小費(fèi)用流(MCF)剖面線出現(xiàn)大量不連續(xù)點(diǎn)；質(zhì)量引導(dǎo)掩膜切線出現(xiàn)突變值，形成兩個(gè)波峰，不符合整體走勢(shì)。無權(quán)多網(wǎng)絡(luò)、PCG預(yù)解共軛梯度雖然整體光滑，但是與其他剖面線有一定差別，可能在解纏過程中引入了誤差，導(dǎo)致解纏結(jié)果出現(xiàn)問題。

圖3 不同解纏算法的剖面圖比較

4.3 定量統(tǒng)計(jì)比較及運(yùn)行效率對(duì)比

對(duì)各相位解纏算法得到的解纏結(jié)果進(jìn)行了統(tǒng)計(jì)分析，分別統(tǒng)計(jì)了其最大值、最小值、均值及標(biāo)準(zhǔn)差等，表2列出了各種相位解纏方法的統(tǒng)計(jì)結(jié)果。

結(jié)果顯示加權(quán)多級(jí)網(wǎng)絡(luò)和LP最小范數(shù)法的平均值和標(biāo)準(zhǔn)差較大，其他方法大致接近。在運(yùn)行速度上，枝切法運(yùn)行時(shí)間最短，LP最小范數(shù)運(yùn)行時(shí)間最長(zhǎng)。枝切法運(yùn)行速度最快，但是穩(wěn)定性較差；區(qū)域生長(zhǎng)法穩(wěn)定性較好，但是運(yùn)行時(shí)間過長(zhǎng)；無權(quán)多級(jí)網(wǎng)絡(luò)運(yùn)行較快且穩(wěn)定性較好；最小不連續(xù)法在解纏效果、標(biāo)準(zhǔn)差、運(yùn)行時(shí)間、直方圖形狀和剖面線走勢(shì)上比較均衡，沒有突出優(yōu)勢(shì)，但是整體較好；LP最小范數(shù)穩(wěn)定性一般、運(yùn)行最慢，但是解纏效果較好。

本文采用的9種相位解纏算法均運(yùn)行在Windows XP環(huán)境下，電腦配置為Intel(R) Core(TM)2 Duo CPU2.80GHz、4G內(nèi)存。各種解纏算法運(yùn)行時(shí)間如表2所示。從表2中可以看出，不同的相位解纏算法的運(yùn)行時(shí)間有很大差異，其中枝切法最快，而LP最小范數(shù)法由于是一種全局最優(yōu)解，其運(yùn)行時(shí)間相對(duì)較長(zhǎng)，是所有算法中最慢的一種方法。

4.4 相位解纏對(duì)形變結(jié)果的比較

以上的統(tǒng)計(jì)及分析還不能確定相位解纏算法在礦區(qū)大形變量解纏時(shí)的優(yōu)劣，還需要考慮解纏算法可靠性的問題。

在差分干涉圖中，一個(gè)條紋代表了半個(gè)波長(zhǎng)的形變量，因此，根據(jù)條紋的數(shù)量可以得到最大形變量的大小，通過不同相位解纏算法得到的礦區(qū)最大形變量與條紋數(shù)目直接反應(yīng)的最大形變量的差值應(yīng)該在較小的范圍內(nèi)，表3列出了各相位解纏方法得到的最大形變量與理論值之間的差值。

表2 相位解纏算法統(tǒng)計(jì)比較

表3 各種解纏方法最大形變量與理論值的差異

從表3中可以看出，盡管所有的相位解纏算法都能恢復(fù)相位的主值并且得到形變量，但在得到的最大形變量方面卻存在較大的差異。并且有的算法與理論值相差較大，在這些算法中，區(qū)域生長(zhǎng)法、無權(quán)多網(wǎng)絡(luò)算法、PCG預(yù)解共軛梯度方法得到的形變量與理論值相差較小。

通過分析各種解纏方法解纏后的形變區(qū)域的直方圖，發(fā)現(xiàn)區(qū)域生長(zhǎng)法、最小費(fèi)用流法(MCF)的直方圖較平滑，分布較好；枝切法、質(zhì)量引導(dǎo)掩膜切線、最小不連續(xù)法分布雖然較好，但是零值較高；LP最小范數(shù)、加權(quán)多級(jí)網(wǎng)絡(luò)、PCG預(yù)解共軛梯度直方圖分布不好且有毛刺；無權(quán)多網(wǎng)絡(luò)出現(xiàn)兩個(gè)峰值。

5 結(jié)束語

對(duì)梁寶寺礦區(qū)的差分干涉圖，采用多種方法進(jìn)行了解纏實(shí)驗(yàn)，并對(duì)解纏結(jié)果進(jìn)行了定性和定量分析，通過分析發(fā)現(xiàn)：在礦區(qū)環(huán)境中，最小費(fèi)用流法沒有顯著優(yōu)勢(shì)，反而出現(xiàn)了較多不可解纏的區(qū)域，而最小不連續(xù)法、LP最小范數(shù)、無權(quán)多級(jí)網(wǎng)絡(luò)和PCG預(yù)解共軛梯度在解纏效果相對(duì)不錯(cuò)。其中PCG在解纏過程中引入誤差導(dǎo)致直方圖和剖面線與整體不一致，LP最小范數(shù)運(yùn)行時(shí)間過長(zhǎng)。綜合分析得出最適于用于礦區(qū)的相位解纏方法是最小不連續(xù)法和無權(quán)多級(jí)網(wǎng)絡(luò)法，它們解纏效率較高、精度較高，穩(wěn)定性和可靠性較好，在礦區(qū)解纏實(shí)驗(yàn)中，受噪聲影響較小，能較好地完成相位解纏。本文較全面地分析了在礦區(qū)環(huán)境中各種相位解纏方法的優(yōu)缺點(diǎn)，為礦區(qū)相位解纏方法的選擇提供了依據(jù)，對(duì)礦區(qū)沉降監(jiān)測(cè)精度的提高具有重要意義。

致謝：本研究還得到了山東省泰山學(xué)者建設(shè)工程專項(xiàng)以及山東科技大學(xué)科研創(chuàng)新團(tuán)隊(duì)支持計(jì)劃項(xiàng)目(2011KYTD103)的資助，在此表示感謝。

參考文獻(xiàn)：

[1] HANSSEN R.Radar interferometry-data interpretation and analysis[M].Kluwer Academic Publisher，2001.

[2] 李平湘，楊杰.雷達(dá)干涉測(cè)量原理與應(yīng)用[M].北京：測(cè)繪出版社，2006.

[3] 王超，張紅.星載合成孔徑雷達(dá)干涉測(cè)量[M].北京：科學(xué)出版社，2002.

[4] 張永紅，張繼賢，龔文瑜.基于SAR干涉點(diǎn)目標(biāo)分析技術(shù)的城市地表形變監(jiān)測(cè)[J].測(cè)繪學(xué)報(bào)，2009，38(6)：482-487.

[5] 程璞，許才軍，王華.InSAR相位解纏算法研究[J].大地測(cè)量與地球動(dòng)力學(xué)，2007，27(3)：50-55.

[6] 黃國(guó)滿，張繼賢，趙爭(zhēng)，等.機(jī)載干涉SAR測(cè)繪制圖應(yīng)用系統(tǒng)研究[J].測(cè)繪學(xué)報(bào)，2008，37(3)：277-279.

[7] 王志勇，張繼賢，黃國(guó)滿.基于InSAR技術(shù)的濟(jì)寧礦區(qū)沉降精細(xì)化監(jiān)測(cè)與分析[J].中國(guó)礦業(yè)大學(xué)學(xué)報(bào)，2014，43(1)：169-174.

[8] 靳國(guó)旺，徐國(guó)華，余懋勛，等.基于瞬時(shí)頻率估計(jì)的InSAR相位解纏[J].測(cè)繪科學(xué)技術(shù)學(xué)報(bào)，2009，26(1)：33-35.

[9] GHIGLIA D C.Two-dimensional phase unwrapping：Theory，algorithms and software[M].New York：John Wiley & Sons，Inc，1998.

[10] GOLDSTEIN R M，ZEBKER H A，WERNER C L.Satellite radar interferometry：Two-dimensional phase unwrapping[J].Radio Science，1988，23(4)：713-720.

[11] CHEN C W.Statistical-cost network-flow approaches to two-dimensional phase unwrapping for radar interferometry[D].California：Standford University，2001.

[12] 劉國(guó)林，郝華東，陶秋香.卡爾曼濾波相位解纏及其與其他方法的對(duì)比分析[J].武漢大學(xué)學(xué)報(bào)(信息科學(xué)版)，2010，35(10)：1174-1178.

[13] 趙爭(zhēng)，張繼賢，張過.遺傳算法在InSAR相位解纏中的應(yīng)用[J].測(cè)繪科學(xué)，2002，27(3)：37-39.

[14] 劉志敏，張景發(fā)，羅毅，等.InSAR相位解纏算法的實(shí)驗(yàn)對(duì)比研究[J].遙感信息，2012，27(2)：71-76.