劉付剛， 張洪全， ?？■危?江曉林， 董 軍

(黑龍江科技大學(xué) 電子與信息工程學(xué)院，哈爾濱 150022)

基于D-InSAR的礦區(qū)地面沉降監(jiān)測(cè)

劉付剛， 張洪全， 常俊鑫， 江曉林， 董 軍

(黑龍江科技大學(xué) 電子與信息工程學(xué)院，哈爾濱 150022)

為獲得礦區(qū)地表的具體形變量和形變趨勢(shì)，根據(jù)礦區(qū)地表沉陷形變特征，分析了合成孔徑雷達(dá)差分干涉測(cè)量(D-InSAR)技術(shù)中SAR影像配準(zhǔn)、去平地效應(yīng)、干涉圖濾波以及相位解纏等關(guān)鍵問(wèn)題。結(jié)合實(shí)際礦區(qū)SAR影像數(shù)據(jù)，利用相應(yīng)技術(shù)的經(jīng)典處理算法，獲取了煤炭開(kāi)采引起的地面沉降數(shù)據(jù)。實(shí)驗(yàn)結(jié)果與實(shí)地勘察結(jié)論基本一致，證明D-InSAR技術(shù)在地面沉降監(jiān)測(cè)方面具有強(qiáng)大優(yōu)勢(shì)。該研究可為礦區(qū)生產(chǎn)、整體規(guī)劃與數(shù)字化礦山建設(shè)提供技術(shù)支持。

D-InSAR； 去平地效應(yīng)； 干涉濾波； 地面沉降監(jiān)測(cè)

0 引 言

煤炭的不合理開(kāi)采破壞了礦區(qū)原有的地質(zhì)結(jié)構(gòu)，極易引起地面形變和沉陷，導(dǎo)致人類(lèi)生存環(huán)境質(zhì)量的嚴(yán)重破環(huán)。常規(guī)的礦區(qū)地面沉降監(jiān)測(cè)主要有GPS定位和水準(zhǔn)測(cè)量，由于礦區(qū)監(jiān)測(cè)地域廣、地表形變漸變性和長(zhǎng)期性等特點(diǎn)，導(dǎo)致地表形變監(jiān)測(cè)和預(yù)測(cè)效果不佳，甚至失效。D-InSAR技術(shù)是20世紀(jì)90年代發(fā)展起來(lái)的一種可以高精度監(jiān)測(cè)大面積微小地面形變的新技術(shù)。在礦區(qū)地面監(jiān)測(cè)應(yīng)用方面，各發(fā)達(dá)國(guó)家科研人員利用D-InSAR技術(shù)取得了一系列的突出成果。Marco等于1995年利用ERS數(shù)據(jù)監(jiān)測(cè)美國(guó)Belridge油田的地面形變，并將監(jiān)測(cè)分析結(jié)果同傳統(tǒng)水準(zhǔn)法對(duì)比，在年均30～40 cm的下沉量中，最大高程誤差不足5 mm；此后，世界各國(guó)學(xué)者針對(duì)不同的煤礦開(kāi)采區(qū)，如德國(guó)圖林根州鈾礦開(kāi)采區(qū)、波蘭Upper Silesian煤田開(kāi)采區(qū)、澳大利亞某煤田開(kāi)采區(qū)等礦區(qū)的地面沉降進(jìn)行監(jiān)測(cè)分析研究[1]，得到了毫米級(jí)的垂直方向形變數(shù)據(jù)，進(jìn)一步證明了D-InSAR技術(shù)在礦區(qū)地面形變沉降監(jiān)測(cè)中的作用和良好的監(jiān)測(cè)效果。近年來(lái)，我國(guó)科研工作者利用D-InSAR技術(shù)開(kāi)展礦區(qū)地面沉降研究也取得了可喜的成果[2-5]。2017年3月，李振洪團(tuán)隊(duì)利用我國(guó)“高分三號(hào)”衛(wèi)星獲取了高質(zhì)量重軌干涉合成孔徑雷達(dá)實(shí)驗(yàn)數(shù)據(jù)，生成了我國(guó)第一幅衛(wèi)星干涉SAR影像，并從影像中提取到了亞厘米級(jí)地面沉降沉息，開(kāi)啟了我國(guó)自主雷達(dá)衛(wèi)星InSAR形變監(jiān)測(cè)的新進(jìn)程[6-7]。筆者針對(duì)合成孔徑雷達(dá)差分干涉測(cè)量中關(guān)鍵問(wèn)題進(jìn)行分析和討論，并結(jié)合實(shí)際礦區(qū)SAR影像數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，以驗(yàn)證D-InSAR在礦區(qū)地表形變監(jiān)測(cè)方面的應(yīng)用效果。

1 D-InSAR成像基本原理

D-InSAR技術(shù)有別于傳統(tǒng)的基于點(diǎn)觀測(cè)的地面測(cè)量方法，該技術(shù)融合了合成孔徑雷達(dá)成像技術(shù)和干涉測(cè)量技術(shù)，利用單軌道雙天線(xiàn)模式(兩副天線(xiàn)同時(shí)觀測(cè))或單天線(xiàn)重復(fù)軌道模式(兩次近平行觀測(cè))獲得同一地區(qū)的兩景數(shù)據(jù)，通過(guò)獲取同一目標(biāo)對(duì)應(yīng)的兩個(gè)回波信號(hào)之間的相位差，結(jié)合軌道數(shù)據(jù)來(lái)獲取高精度、高分辨率的地面高程信息，精確測(cè)量地表某一點(diǎn)的三維空間位置及微小變化。星載重復(fù)軌道成像示意如圖1所示[8]。

圖1中b為基線(xiàn)距，基線(xiàn)距與水平方向的傾角為α,S1和S2分別表示主輔圖像傳感器，主圖像入射角為θ，H為主傳感器相對(duì)地面高度，R和R+ΔR分別為主輔圖像斜距，P為地面目標(biāo)點(diǎn)，其高程為h。設(shè)地面目標(biāo)點(diǎn)P兩次成像時(shí)的圖像分別為c1和c2,則可表示為

c1=A1ejφ(R),

(1)

c2=A2ejφ(R+ΔR),

(2)

式中:A1、A2——c1、c2的振幅；

φ(R)、φ(R+ΔR)——c1、c2的相位。

圖1 InSAR成像示意

忽略大氣效應(yīng)等因素，可以認(rèn)為相位信息由兩部分組成，一是與信號(hào)往返路徑相對(duì)應(yīng)的相位，二是散射特性引起的相位。

假設(shè)主輔圖像傳感器成像散射特性不變，對(duì)兩圖像信號(hào)進(jìn)行復(fù)共扼相乘可得到復(fù)干涉紋圖，即

(3)

式中：*——共軛； λ——雷達(dá)波的波長(zhǎng)。

由式(1)～(3)可以得到干涉圖的相位為

(4)

根據(jù)圖1所示，利用三角函數(shù)可得

(5)

z=H-Rcosθ,

(6)

展開(kāi)式(5)，由于ΔR遠(yuǎn)遠(yuǎn)小于R，故可得

ΔR≈bsin(θ-α),

(7)

將式(7)帶入(4)，可得

(8)

對(duì)干涉相位進(jìn)行解纏可以得到式(8)中的N?，F(xiàn)只考慮干涉圖相位的主值有

(9)

將式(9)帶入到(6)得

至此，可以求出目標(biāo)點(diǎn)的高度值h。由此可見(jiàn)，InSAR技術(shù)屬于主動(dòng)微波遙感系統(tǒng)，工作機(jī)理和方式與傳統(tǒng)光學(xué)不同，它不依賴(lài)于可見(jiàn)光, 能全天時(shí)全天候工作,精度相對(duì)較高。

2 研究區(qū)概況和數(shù)據(jù)來(lái)源

濟(jì)寧市北部礦井，煤田南北長(zhǎng)4 km，東西寬6 km。礦區(qū)內(nèi)地勢(shì)東高西低，地貌較為復(fù)雜，地面標(biāo)高均值是+38 m，氣候溫和，歷年最大積雪0.15 m，最大凍土0.31 m。兗石鐵路和新兗鐵路貫穿東西，交通便利。文中所用數(shù)據(jù)為歐空局ENVISAT衛(wèi)星提供的2008年12月至2009年7月的6景C波段ASAR數(shù)據(jù)，影像覆蓋濟(jì)寧市及北部煤田區(qū)域。影像數(shù)據(jù)所采用的雷達(dá)為ASAR型，軌道號(hào)為35571和35573～35577，軌跡為447，幀號(hào)為2889。

3 D-InSAR關(guān)鍵技術(shù)分析

D-InSAR的關(guān)鍵技術(shù)包括圖像配準(zhǔn)、產(chǎn)生干涉圖、濾波、平地效應(yīng)消除、相位解纏和差分干涉處理等。D-InSAR的經(jīng)典方法為三軌法，采用三景雷達(dá)圖像，以其中之一作為主圖像，另外兩景為從圖像，產(chǎn)生干涉圖。圖2給出了三軌法差分干涉測(cè)量數(shù)據(jù)流程。

圖2 三軌法差分干涉測(cè)量數(shù)據(jù)流程

3.1 SAR圖像配準(zhǔn)

兩幅圖像在空間和灰度上的映射即為圖像配準(zhǔn)。在選取公共主影像之后，為得到干涉條紋圖，必須使得主從兩幅SAR圖像上的點(diǎn)對(duì)應(yīng)地面上的同一點(diǎn)，即將從影像配準(zhǔn)到主影像的影像坐標(biāo)系內(nèi)。由于地形對(duì)相位的影響非常敏感，加之時(shí)間基線(xiàn)和空間基線(xiàn)的影響，很容易導(dǎo)致一些影像由于相干性過(guò)低而無(wú)法得到較高的圖像配準(zhǔn)質(zhì)量。

圖像對(duì)的配準(zhǔn)包括粗配準(zhǔn)和精配準(zhǔn)兩個(gè)步驟。粗配準(zhǔn)就是利用衛(wèi)星軌道參數(shù)或人工選取少量的控制點(diǎn)計(jì)算主、輔圖像之間在方位向和距離向上的坐標(biāo)映射關(guān)系，其配準(zhǔn)精度約為30個(gè)像元。粗配準(zhǔn)常采用基于衛(wèi)星DEOS精密軌道參數(shù)的配準(zhǔn)方法[9]。精配準(zhǔn)又分為像元級(jí)配準(zhǔn)和亞像元級(jí)配準(zhǔn)。兩者的基本步驟和方法大體相同，一般來(lái)說(shuō)，亞像元級(jí)圖像配準(zhǔn)精度要優(yōu)于1/8個(gè)像元才能符合SAR干涉處理的精度要求，此時(shí)所造成的失相干很小(約為4%)。

3.2 去平地效應(yīng)

影像配準(zhǔn)后對(duì)應(yīng)的點(diǎn)共軛相乘，可以獲得干涉相位圖。干涉相位圖并不能直觀地體現(xiàn)地形高低的變化，這是因?yàn)楦缮嫦辔粓D的相位信息里包含著平地(高度不變)所對(duì)應(yīng)的相位量。這種由高度不變的平地引起干涉相位在距離向和方位向呈周期性變化的現(xiàn)象稱(chēng)為平地效應(yīng)，它在時(shí)域內(nèi)顯示沿方位向平行的明暗相間的密集條紋。平地效應(yīng)所引起的條紋和地形高低變化所引起的條紋相互疊加，掩蓋了地形高低變化對(duì)應(yīng)的相位信息，增加了相位解纏難度，必須進(jìn)行去除平地效應(yīng)處理。平地相位去除后，干涉圖的相位信息就能反映真實(shí)的地表高程相位和形變相位。針對(duì)ASAR數(shù)據(jù)，衛(wèi)星搭載的平臺(tái)飛行軌道穩(wěn)定性較ERS-I衛(wèi)星要差，飛行軌道之間存在著一定的夾角，其方位向也存在著平地效應(yīng)。采用精密軌道參數(shù)和成像區(qū)域中心點(diǎn)位置去除平地效應(yīng)干涉圖如圖3所示。

a 去除前

b 去除后

由圖3可見(jiàn)，去平地效應(yīng)之后的干涉紋圖條紋密度較去除前稀疏了很多，并且條紋的清晰度有明顯改善，去平地效應(yīng)后，有利于噪聲濾波和相位解纏的進(jìn)行。

3.3 干涉圖濾波

考慮到SAR圖像噪聲(通常定義為滑動(dòng)窗口內(nèi)的異常最大值或最小值)對(duì)圖像重構(gòu)的影響，還需要對(duì)成像后干涉圖進(jìn)行濾波處理。目前常用的濾波方法分為空間域?yàn)V波和頻率域?yàn)V波[10]。空間域?yàn)V波方法是利用圖像像素的空間相關(guān)性對(duì)相干斑進(jìn)行濾波，如多視平均處理方法、均值濾波處理方法、中值濾波處理方法等。頻率域的方法主要是利用小波變換，比較著名的有基于小波變換和多尺度分析的濾波方法。理想的SAR圖像濾波算法應(yīng)該能夠在抑制斑點(diǎn)噪聲的同時(shí)保持圖像的邊緣和紋理信息。采用小波變換的濾波方法得到的干涉圖如圖4所示。

a 濾波前

b 濾波后

對(duì)比圖4a、b可以看出，采用小波變換的濾波方法濾波具有一定的效果，一些噪聲明顯濾除，濾波后的干涉圖效果明顯改善。利用工具Cursor Location/Value，查看圖4的像元值，可以看到相干性系數(shù)分布在0～1，圖像中越白的點(diǎn)，相干性系數(shù)值越大，說(shuō)明該區(qū)域的相干性越高。圖5給出了截取濾波前后兩幅圖中第250行數(shù)據(jù)的剖面。由圖5可見(jiàn)，經(jīng)過(guò)濾波去噪處理后圖像曲線(xiàn)變化平穩(wěn)，數(shù)據(jù)處理過(guò)程中濾掉了大量的噪聲數(shù)據(jù)。

a 濾波前

b 濾波后

3.4 相位解纏

去平地效應(yīng)和濾波處理之后的干涉圖中，相位值并非真實(shí)的絕對(duì)相位值，而是以非線(xiàn)性的方式混合纏繞了2π的整數(shù)倍，如式(4)所示，假設(shè)二維矩陣φi,j表示纏繞相位，為了解纏該二維矩陣，對(duì)矩陣中的每個(gè)點(diǎn)(i,j)加減2π的整數(shù)倍，可以得到一個(gè)連續(xù)的解纏相位函數(shù)φi,j,即

φi,j=φi,j+2πki,j，i[0,M-1]，j[0,N-1],

其中,Ki,j為整數(shù)，-π≤φi，j<π。

相位解纏的準(zhǔn)確性直接關(guān)系到地表形變量的測(cè)量精度。目前，相位解纏方法主要分為三類(lèi)：路徑跟蹤法、最小范數(shù)法基于網(wǎng)絡(luò)規(guī)劃的相位解纏算法。利用加權(quán)最小二乘法解纏干涉圖如圖6。

a 解纏前

b 解纏后

圖6中解纏后相位圖中干涉條紋明顯比原始數(shù)據(jù)少，其相位值是一個(gè)比原始纏繞相位值更加平滑的估計(jì)值，解纏后的相位值更接近真實(shí)的相位值。

3.5 實(shí)驗(yàn)結(jié)果分析

綜合上述實(shí)驗(yàn)數(shù)據(jù)及實(shí)驗(yàn)過(guò)程，根據(jù)圖4 的干涉圖可以得出礦區(qū)局部地區(qū)的地表形變結(jié)果。在2008年12月至2009年7月之間，煤田北部東西長(zhǎng)0.8 km，南北長(zhǎng)1.5 km塌陷較為突出，最大塌陷處達(dá)3.0 cm，其他地區(qū)的沉降量處于[0.7,2.2]區(qū)間內(nèi)，沉降趨勢(shì)較緩慢。最大沉降地東北方向1.4 km處有兩個(gè)村莊，地下煤層的開(kāi)采引起的地質(zhì)形變將會(huì)影像到村莊的建筑物，需要加強(qiáng)注漿充填開(kāi)采的力度。經(jīng)過(guò)實(shí)地勘察和調(diào)查資料對(duì)干涉測(cè)量結(jié)果進(jìn)行驗(yàn)證，發(fā)現(xiàn)分析結(jié)果與實(shí)際情況基本一致。由此可見(jiàn)，D-InSAR技術(shù)監(jiān)測(cè)礦區(qū)地表沉降效果較好，與傳統(tǒng)的地表形變監(jiān)測(cè)方法相比較，該技術(shù)具有覆蓋范圍廣、監(jiān)測(cè)成本低、監(jiān)測(cè)效率高的明顯優(yōu)勢(shì)。

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

文中采用歐空局ENVISAT衛(wèi)星提供的C波段ASAR數(shù)據(jù)，利用D-InSAR技術(shù)中的典型算法對(duì)濟(jì)寧北部煤田區(qū)域進(jìn)行了地表形變監(jiān)測(cè)，并將實(shí)驗(yàn)結(jié)果通過(guò)實(shí)地勘察和資料查證進(jìn)行驗(yàn)證，得出了與實(shí)際結(jié)果一致的結(jié)論。通過(guò)分析可知，D-InSAR技術(shù)在礦區(qū)地表沉降形變監(jiān)測(cè)應(yīng)用中完全可行，并且具有明顯的技術(shù)優(yōu)勢(shì)，利用該技術(shù)能夠?yàn)榈V產(chǎn)資源的開(kāi)采提供科學(xué)合理的決策參考，同時(shí)可以科學(xué)預(yù)測(cè)可能出現(xiàn)的地質(zhì)災(zāi)害，保證和諧礦區(qū)、綠色礦區(qū)的可持續(xù)發(fā)展。可以預(yù)見(jiàn)，未來(lái)隨著SAR衛(wèi)星的不斷專(zhuān)業(yè)化、InSAR處理軟件和算法的日趨成熟，D-InSAR技術(shù)的實(shí)際應(yīng)用水平將會(huì)越來(lái)越高。

[1] Carnec C，Massonnet D，King C. Two examples of the use of SAR interferometry on displacement of small spatial extent[J].Geophysical Research Letters, 1996，23(4): 3579-3582.

[2] 楊澤發(fā), 易輝偉, 朱建軍, 等. 基于InSAR時(shí)序形變的礦區(qū)全盆地沉降時(shí)空演化規(guī)律分析[J]. 中國(guó)有色金屬學(xué)報(bào), 2016, 26(7): 1515-1522.

[3] 陸燕燕, 柯長(zhǎng)青, 陳德良, 等. PS-InSAR在沛縣礦區(qū)地表沉降監(jiān)測(cè)中的應(yīng)用[J]. 地理空間信息, 2016, 14(5): 96-99.

[4] Saghafi A, Pinetown K, Javanmard H. Developing a new method to identify the source of gas emissions into longwall and goaf from sSurrounding strata[C]// Coal Operators’ Conference 2015. Wollongong, Australia: The University of Wollongong, The Australasian Institute of Mining and Metallurgy and Mine Managers Association of Australia, 2015.

[5] 汪寶存, 郭凌飛, 王軍見(jiàn), 等. 礦區(qū)地表形變InSAR監(jiān)測(cè)——以永城市為例[J]. 測(cè)繪與空間地理信息, 2016, 39(6): 24-27.

[6] 丁赤飚, 劉佳音, 雷 斌, 等. 高分三號(hào)SAR衛(wèi)星系統(tǒng)級(jí)幾何定位精度初探[J]. 雷達(dá)學(xué)報(bào), 2017, 6(1): 11-16.

[7] 劉淮宇. 高分專(zhuān)項(xiàng)應(yīng)用推廣全面啟動(dòng)[J]. 太空探索, 2016(2): 7-7.

[9] Rudzinski L, Mirek K, Mirek J. Source mechanism analysis of strong mining induced seismic event and its influence on ground deformation observed by InSAR technique[C]//EGU General Assembly 2016, Vienna, Austria： European Geophysical Society and the European Union of Geoscienes, 2016，18: 6759.

[10] Chaussard E, Kerosky S. Characterization of black sand mining activities and their environmental impacts in the philippines using remote sensing[J]. Remote Sensing, 2016, 8(2): 100-116.

(編校 王 冬)

D-InSAR-based monitoring of land subsidence in mining area

LiuFugang,ZhangHongquan,ChangJunxin,JiangXiaolin,DongJun

(School of Electronics & Information Engineering, Heilongjiang University of Science & Technology, Harbin 150022, China)

This paper provides an insight into the specific mining surface deformation and deformation trend. The research building on the characteristics behind subsidence deformation in mines, works towards the following objectives：analyzing the key issues, such as SAR image registration, flat-earth removal, interferogram filtering, and phase unwrapping, all of which are involved in synthetic aperture radar differential interferometry(D-InSAR)；combined with SAR image data derived from actual mining area, employing the classical processing algorithm of corresponding technology to carry out tests and thereby derive the data on land subsidence caused by coal mining; and providing the results basically consistent with the field investigation results to verify the strong advantage D-InSAR technology has in land subsidence monitoring. The research could provide technical support for mining area production, overall planning and digital mine construction.

D-InSAR; flat-earth removal; interference filter; ground subsidence monitoring

2017-04-03

黑龍江省教育廳科學(xué)技術(shù)研究項(xiàng)目(12541716)

劉付剛(1981-)，男，山東省德州人，副教授，博士，研究方向：信號(hào)檢測(cè)處理、InSAR遙感遙測(cè)，E-mail：liufugang_36@163.com。

10.3969/j.issn.2095-7262.2017.03.013

TD173.4

2095-7262(2017)03-0265-05

A