雙軌D-InSAR在礦區(qū)地面沉降監(jiān)測(cè)中的應(yīng)用*

劉同文，高騰飛，于廣婷，徐進(jìn)達(dá)

(1.山東省地質(zhì)測(cè)繪院，山東 濟(jì)南 250002； 2.山東科技大學(xué)，山東 青島 266590)

雙軌D-InSAR在礦區(qū)地面沉降監(jiān)測(cè)中的應(yīng)用*

劉同文1，高騰飛2，于廣婷1，徐進(jìn)達(dá)1

(1.山東省地質(zhì)測(cè)繪院，山東 濟(jì)南 250002； 2.山東科技大學(xué)，山東 青島 266590)

合成孔徑雷達(dá)差分干涉測(cè)量技術(shù)(Differential Interferometry Synthetic Aperture Radar，D-InSAR)已成為監(jiān)測(cè)礦區(qū)地表沉降的有力工具之一。文章選取2景2004-12-26—2005-01-30覆蓋濟(jì)寧礦區(qū)的C波段的ENVISAT ASAR影像作為干涉對(duì)，首先對(duì)整景影像進(jìn)行差分干涉處理，得到D-InSAR 一系列結(jié)果圖，然后選取沉降明顯的葛亭礦區(qū)進(jìn)行精細(xì)化研究，結(jié)合GIS軟件疊加數(shù)字正射影像圖 DOM、開(kāi)采平面圖和礦區(qū)邊界圖等進(jìn)行空間分析。結(jié)果表明：雙軌D-InSAR技術(shù)可以對(duì)礦區(qū)的沉降分布進(jìn)行有效的探測(cè)與定位，且能對(duì)礦區(qū)的沉降程度進(jìn)行監(jiān)測(cè)。

礦區(qū)；地面沉降；合成孔徑雷達(dá)差分干涉測(cè)量；GAMMA

Abstract：D-InSAR is one of the useful techniques for monitoring the ground subsidence.Two C-land ENVISAT ASAR images from December 26th，2004 to January 30th，2005 were selected as interferometry pair.First of all，the total image were processed and a series of maps were obtained，then the DOM，mining plan and mining boundary map in Geting mining area were performed for the spatial analysis in GIS.The results show that the technology of double-track D-InSAR can be used to effectively detect and locate the subsidence distribution in the mining area，and monitor the degree of subsidence in the mining area.

Keywords：mining area；ground subsidence；D-InSAR；GAMMA

0 引言

目前，國(guó)內(nèi)煤礦開(kāi)采的地面沉降監(jiān)測(cè)方法主要有傳統(tǒng)大地測(cè)量、精密水準(zhǔn)測(cè)量、近景攝影測(cè)量、GPS、三維激光掃描等，地表移動(dòng)變形觀測(cè)站的設(shè)置以剖面線(xiàn)形式為主，監(jiān)測(cè)數(shù)據(jù)為離散監(jiān)測(cè)點(diǎn)的形變信息，只能對(duì)監(jiān)測(cè)點(diǎn)繪制出的沉降剖面曲線(xiàn)圖進(jìn)行分析[1]。如果要對(duì)沉降區(qū)地面沉降情況進(jìn)行全面把握和分析處理，必須布設(shè)大量的監(jiān)測(cè)點(diǎn)，通過(guò)擬合、插值等手段，繪制出沉降等值線(xiàn)圖，進(jìn)而為合理安排礦區(qū)開(kāi)采，預(yù)測(cè)、評(píng)估沉降損失等提供依據(jù)，顯然傳統(tǒng)的測(cè)量方法和手段無(wú)能為力[2]。合成孔徑雷達(dá)差分干涉測(cè)量(D-InSAR)是在InSAR的基礎(chǔ)上發(fā)展起來(lái)的，大量研究與應(yīng)用實(shí)例表明，其能夠應(yīng)用于礦區(qū)長(zhǎng)期緩慢的地表形變監(jiān)測(cè)。2005年，劉國(guó)林等[3]對(duì)合成孔徑雷達(dá)干涉測(cè)量與全球定位系統(tǒng)數(shù)據(jù)融合監(jiān)測(cè)礦區(qū)地表沉降的可行性進(jìn)行了分析。2007 年，董玉森等[4]收集了1992年12月至1998年6月的JERS-1 L波段雷達(dá)數(shù)據(jù)，利用D-InSAR 技術(shù)進(jìn)行礦區(qū)地表沉降監(jiān)測(cè)研究。2012年，陶秋香等[5]用ALOS PALSAR和ENVISAT ASAR數(shù)據(jù)分析討論了L和C波段雷達(dá)干涉數(shù)據(jù)的礦區(qū)地面沉降監(jiān)測(cè)能力，試驗(yàn)結(jié)果表明L波段雷達(dá)干涉數(shù)據(jù)更容易監(jiān)測(cè)到沉降梯度和沉降量較大的礦區(qū)地面沉降。2013年，劉一霖[6]針對(duì)礦區(qū)地表大量級(jí)突變沉陷形變?cè)趥鹘y(tǒng)InSAR監(jiān)測(cè)中的問(wèn)題及對(duì)礦震災(zāi)害的監(jiān)測(cè)，進(jìn)行了一系列的實(shí)驗(yàn)分析。2014 年，陳炳乾等[7]利用 D-InSAR對(duì)陜西省大柳塔礦進(jìn)行開(kāi)采沉陷監(jiān)測(cè)，并將監(jiān)測(cè)結(jié)果與支持向量回歸算法(SVR)相結(jié)合進(jìn)行開(kāi)采沉陷動(dòng)態(tài)預(yù)計(jì)；2016年，郭山川等[8]對(duì)黃土高原礦區(qū)使用D-InSAR技術(shù)監(jiān)測(cè)其復(fù)雜劇烈的動(dòng)態(tài)的地表形變。

基于此，本文闡述和分析了雙軌D-InSAR獲取地表形變信息的基本原理及基于GAMMA雷達(dá)數(shù)據(jù)處理軟件的數(shù)據(jù)處理流程，選取2景2004-12-26—2005-01-30覆蓋濟(jì)寧礦區(qū)的C波段的ENVISAT ASAR作為干涉對(duì)，外部DEM選為SRTM3 DEM，基于GAMMA軟件進(jìn)行雙軌差分干涉處理，將得到的沉降圖經(jīng)過(guò)地理編碼后導(dǎo)入 ArcGIS 軟件中與其它空間信息進(jìn)行疊加來(lái)分析礦區(qū)地面沉降。首先對(duì)試驗(yàn)區(qū)所在的整景圖像進(jìn)行差分干涉處理，得到D-InSAR 一系列結(jié)果圖(差分干涉圖、增強(qiáng)后的差分干涉圖、相位解纏圖、沉降圖等)。然后選取沉降較為明顯的葛亭礦區(qū)在該時(shí)間段內(nèi)的開(kāi)采沉降進(jìn)行精細(xì)化監(jiān)測(cè)研究，得到葛亭礦區(qū)的沉降量、沉降位置和分布等信息，并結(jié)合GIS軟件疊加數(shù)字正射影像圖 DOM、開(kāi)采平面圖和礦區(qū)邊界圖等進(jìn)行空間分析，以沉降剖面圖進(jìn)行顯示，直觀、清晰地顯示出礦區(qū)地面沉降的位置和范圍。

1 雙軌D-InSAR基本原理及數(shù)據(jù)處理流程

根據(jù)地形信息去除方法的不同，D-InSAR 測(cè)量地表形變的方法主要包括：雙軌法、三軌法和四軌法。本文在已獲得外部 DEM 的基礎(chǔ)上，使用雙軌法對(duì)所選數(shù)據(jù)進(jìn)行差分干涉處理。其基本思想是利用試驗(yàn)區(qū)地表變化前后的兩幅SAR影像經(jīng)過(guò)配準(zhǔn)、重采樣和干涉處理生成干涉條紋圖，再利用已有的DEM數(shù)據(jù)模擬地形相位，從干涉圖中去除模擬的地形相位得到差分干涉圖，對(duì)差分干涉圖進(jìn)行濾波、相位解纏處理，然后進(jìn)行相高轉(zhuǎn)換，即可得到地表形變圖。圖1所示為雙軌法幾何原理示意圖[9]。

圖1 雙軌法幾何原理示意圖Fig.1 Geometric principle of double-track method

圖中，S1和S2分別表示主輔影像傳感器位置，P表示地面目標(biāo)點(diǎn)，假設(shè)SAR衛(wèi)星兩次過(guò)境時(shí)，地面在SAR衛(wèi)星視線(xiàn)向(Line of Sihgt，LOS)上發(fā)生Δr的形變，則形變相位為：

(1)

從式中可以看出干涉相位對(duì)地形變化非常敏感，測(cè)量精度可達(dá)到波長(zhǎng)級(jí)。對(duì)于ENVISAT ASAR數(shù)據(jù)，當(dāng)?shù)乇戆l(fā)生2.8 cm的位移時(shí)，形變相位就會(huì)發(fā)生一個(gè)2π的變化。

基于GAMMA軟件的雙軌D-InSAR數(shù)據(jù)處理流程圖，如圖2所示。

圖2 D-InSAR數(shù)據(jù)處理流程圖Fig.2 D-InSAR data processing flowchart

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

濟(jì)寧礦區(qū)煤炭資源豐富，開(kāi)采煤層厚度大，導(dǎo)致地質(zhì)災(zāi)害分布廣、威脅大，因此對(duì)濟(jì)寧礦區(qū)的地表形變進(jìn)行監(jiān)測(cè)十分必要。本文選取了2景2004-12-26、2005-01-30拍攝的覆蓋濟(jì)寧礦區(qū)的ENVISAT ASAR數(shù)據(jù)，外部參考DEM選為SRTM3 DEM，基于GAMMA軟件進(jìn)行了差分干涉處理，系統(tǒng)研究和分析D-InSAR技術(shù)在礦區(qū)地面沉降監(jiān)測(cè)中的應(yīng)用，為礦區(qū)的安全開(kāi)采和塌陷區(qū)環(huán)境的綜合治理提供科學(xué)依據(jù)。研究區(qū)地理位置如圖3所示，方框內(nèi)為影像覆蓋范圍。

圖3 研究區(qū)地理位置Fig.3 The geographical location in research area

2.1 ENVISAT ASAR影像

ENVISAT衛(wèi)星是歐空局的對(duì)地觀測(cè)衛(wèi)星系列之一，2002年3月l日由阿里亞納5號(hào)火箭發(fā)射的一顆先進(jìn)的太陽(yáng)同步極軌地球環(huán)境監(jiān)測(cè)衛(wèi)星。ENVISAT衛(wèi)星上所載最大設(shè)備是先進(jìn)的合成孔徑雷達(dá)(ASAR)，具有多極化、多入射角、大幅寬等特性，可生成海洋、海岸、極地冰冠和陸地的高質(zhì)量高分辨率影像[10]。

考慮到時(shí)間基線(xiàn)和空間基線(xiàn)對(duì)影像相干性的影響，成像期間植被的影響和形變結(jié)果的可靠性等，最終選定的干涉對(duì)為表1中的影像對(duì)。

表1 ENVISAT ASAR影像對(duì)主要參數(shù)

2.2 外部參考DEM

外部 DEM 的精度對(duì)差分干涉處理結(jié)果的精度高低具有很大的影響，本文試驗(yàn)數(shù)據(jù)處理采用SRTM3 DEM，其平面參考系為WGS-84，垂直參考系為WGS-84 EGM96大地水準(zhǔn)面，絕對(duì)精度為±16 m，相對(duì)高程精度為10 m。濟(jì)寧礦區(qū)位于華北平原，地勢(shì)平坦，SRTM DEM在平原地區(qū)的高程精度比其它地貌區(qū)高，絕對(duì)高程誤差小于5 m。因此，試驗(yàn)中由DEM誤差引入的相位誤差可忽略[11]。圖4所示為選取的外部參考DEM。

3 試驗(yàn)數(shù)據(jù)處理與分析

根據(jù)圖2所示的雙軌D-InSAR數(shù)據(jù)處理流程，基于GAMMA軟件，以2004-12-26獲取的影像為主影像，以2005-01-30獲取的影像為輔影像，對(duì)兩幅SLC影像進(jìn)行配準(zhǔn)、重采樣后得到干涉相位圖，此時(shí)干涉相位包括由參考橢球引入的參考相位，由地形起伏引起的地形相位，由大氣擾動(dòng)引起的大氣延遲相位，由兩次SAR影像成像期間地表形變引起的形變相位及由各種噪聲引起的噪聲相位。為消除干涉相位中的地形相位的影響，使用SRTM3 DEM作為外部參考DEM來(lái)反演地形相位，從

圖4 外部參考DEMFig.4 The external reference DEM

而將地形相位剔除，生成差分干涉圖(如圖5(a)所示)，然后進(jìn)行Goldstein自適應(yīng)濾波處理(如圖5(b)所示)，生成相干系數(shù)圖(如圖5(c)所示，黃色代表高相干性區(qū)域)，使用最小費(fèi)用流和不規(guī)則三角網(wǎng)的方法進(jìn)行相位解纏(如圖5(d)所示)，生成相應(yīng)的形變圖(如圖5(e)所示)。

圖5 D-InSAR數(shù)據(jù)處理的系列結(jié)果圖Fig.5 A series of result maps after D-InSAR data processing

通過(guò)圖5(e)形變圖可以發(fā)現(xiàn)，試驗(yàn)區(qū)存在多個(gè)沉降漏斗，其中A為葛亭礦區(qū)，這些大部分都是由于地下采煤引起的地面沉降所形成的(圖6(a)所示為試驗(yàn)區(qū)的煤礦地理位置)。這充分說(shuō)明了利用雙軌D-InSAR技術(shù)可正確地得到礦區(qū)沉降的整體分布情況，對(duì)礦區(qū)沉陷的位置進(jìn)行探測(cè)與定位。

圖6 試驗(yàn)區(qū)煤礦地理位置及葛亭煤礦放大圖Fig.6 The geographical location in resarch coal mining area and the enlarged drawing of Geting coal mine

2004年12月26日—2005年1月30日監(jiān)測(cè)周期內(nèi)，沉降范圍最廣、沉降量最大的區(qū)域位于葛亭煤礦附近。其中，最大沉降量達(dá)到100 mm，葛亭煤礦是濟(jì)寧地區(qū)沉降情況較為突出的煤礦之一。另外，沉降量大于20 mm的地區(qū)主要集中運(yùn)河煤礦沉降區(qū)和岱莊煤礦沉降區(qū)。這也充分說(shuō)明了D-InSAR技術(shù)可以對(duì)礦區(qū)的沉降程度進(jìn)行監(jiān)測(cè)。

為進(jìn)一步分析由于煤礦開(kāi)采造成的地面沉降情況，對(duì)沉降較嚴(yán)重的葛亭煤礦做進(jìn)一步的精細(xì)化沉降分析。葛亭煤礦的地理坐標(biāo)為116°28′～116°32′E，35°29′～35°32′N(xiāo)，礦井范圍南北長(zhǎng)4～4.5 km，東西寬4～6 km。礦區(qū)內(nèi)地形平坦，地勢(shì)東北略高，西南稍低[11]。葛亭煤礦疊加礦區(qū)邊界、礦區(qū)開(kāi)采平面圖放大圖如圖6(b)所示。

將地理編碼后葛亭礦區(qū)的形變圖導(dǎo)入 ArcGIS 軟件中結(jié)合其它空間信息分析礦區(qū)地面沉降情況，對(duì)礦區(qū)的地面沉降進(jìn)行定性和定量分析。通過(guò)對(duì)D-InSAR 獲得的沉降圖疊加DOM、礦區(qū)開(kāi)采平面和礦區(qū)邊界圖等空間信息，可以驗(yàn)證沉降圖中沉降位置和沉降范圍與地下開(kāi)采活動(dòng)的位置和范圍是否一致，從而可以進(jìn)一步驗(yàn)證 D-InSAR 技術(shù)監(jiān)測(cè)煤礦區(qū)地面沉降的有效性。圖7所示為葛亭礦區(qū)形變圖疊加礦區(qū)開(kāi)采平面、礦區(qū)邊界圖。由圖中可以看出：

1)葛亭礦區(qū)存在不同程度的沉降。在這些沉降中，一些輕微的沉降是由地下水開(kāi)采和數(shù)據(jù)處理過(guò)程中不可避免的誤差引起的，但沉降量比較大的地方是由采礦引起的，與礦區(qū)的實(shí)際情況相符；

2)葛亭礦區(qū)在2004年12月26日—2005年1月30日監(jiān)測(cè)周期內(nèi)最大沉降量為100 mm，在該時(shí)間段內(nèi)存在多個(gè)開(kāi)采工作面，但沉降區(qū)域最為明顯的為2304工作面。2304工作面的回采時(shí)間為2004年10月～2005年04月。

圖7 葛亭礦區(qū)形變圖疊加礦區(qū)開(kāi)采平面、礦區(qū)邊界圖Fig.7 Deformation map in Geting mining area being superimposed in mining plan and boundary map

圖8所示為葛亭煤礦2304工作面沉降圖疊加開(kāi)采平面圖。

圖8 葛亭煤礦2304工作面沉降圖疊加開(kāi)采平面圖Fig.8 Subsidence map of in Geting mining area working face 2304 being superimposed in mining plan

為了清楚的顯示沉降情況，使用 ArcGIS 軟件生成沉降剖面圖，圖9所示為沿著直線(xiàn)A到B和C到D方向繪制的沉降剖面圖。

圖9 沉降剖面圖Fig.9 The subsidence profile

通過(guò)圖9可以看出：在2004年12月26日-2005年01月30日時(shí)間段的沉降圖中，葛亭區(qū)域的2304工作面上方已有明顯沉陷區(qū)域，形成了如圖8所示的沉降漏斗，而且沉降位置和開(kāi)采工作面的位置一致性較好；在整個(gè)監(jiān)測(cè)時(shí)段內(nèi)下沉盆地中心最大下沉值達(dá)到100 mm。由2304開(kāi)采工作面引起的沉降從邊緣到中心依次增大，形成如圖9所示的沉降漏斗，而實(shí)際的礦區(qū)開(kāi)采工作面及其附近區(qū)域的沉降也確實(shí)形似一個(gè)大的漏斗區(qū)。說(shuō)明利用雙軌D-InSAR技術(shù)可以得到具體工作面開(kāi)采引起的沉降分布和沉降量大小。

4 結(jié)論

本文闡述了利用雙軌D-InSAR技術(shù)獲取地表形變信息的基本原理與基于GAMMA軟件的雙軌D-InSAR數(shù)據(jù)處理流程。選取了2景2004年12月26日-2005年01月30日覆蓋濟(jì)寧礦區(qū)的C波段的ENVISAT ASAR作為干涉對(duì)，先獲取了影像覆蓋范圍內(nèi)的礦區(qū)分布情況，后對(duì)沉降量較大的葛亭礦區(qū)進(jìn)行了精細(xì)化的研究與分析，得到如下結(jié)論：

1)在2004年12月26日-2005年01月30日監(jiān)測(cè)周期內(nèi)，試驗(yàn)區(qū)內(nèi)出現(xiàn)了3個(gè)沉降區(qū)，地理坐標(biāo)顯示，它們分別是葛亭煤礦、運(yùn)河煤礦和岱莊煤礦沉降區(qū)。這充分說(shuō)明了利用雙軌D-InSAR技術(shù)可正確地得到礦區(qū)沉降的整體分布情況，對(duì)礦區(qū)沉陷的位置進(jìn)行探測(cè)與定位。其中，葛亭煤礦在監(jiān)測(cè)周期內(nèi)的沉降范圍最廣、沉降量最大。這也充分說(shuō)明了D-InSAR技術(shù)可以對(duì)礦區(qū)的沉降程度進(jìn)行監(jiān)測(cè)。

2)在監(jiān)測(cè)周期內(nèi)，A處葛亭煤礦的沉降最嚴(yán)重，形成了如圖7所示的沉降漏斗。D-InSAR技術(shù)監(jiān)測(cè)到的漏斗中心的最大沉降值為100 mm。葛亭礦區(qū)存在不同程度的沉降，一些輕微的沉降是地下水開(kāi)采和數(shù)據(jù)處理過(guò)程中不可避免的誤差引起的，但沉降量比較大的地方是由采礦引起的，與礦區(qū)的實(shí)際情況相符。

綜上所述，利用雙軌D-InSAR技術(shù)可以對(duì)礦區(qū)沉降分布進(jìn)行探測(cè)與定位，且對(duì)礦區(qū)的沉降程度進(jìn)行監(jiān)測(cè)。

[1] Papadaki E，Tripolitsiotis A，Steiakakis C，et al.Land movement monitoring at the Mavropigi lignite mine using spaceborne D-InSAR[J].Proc Spie，2013，8795(2)：26-28.

[2] Ji M W，Li X J，Wu S C，Gao Y T，Ge L L.Use of SAR interferometry for monitoring illegal mining activities：A case study at Xishimen Iron Ore Mine.Mining Science and Technology，2007，17(2)：262-266.

[3] 劉國(guó)林，張蓮蓬，成樞，等.合成孔徑雷達(dá)干涉測(cè)量與全球定位系統(tǒng)數(shù)據(jù)融合監(jiān)測(cè)礦區(qū)地表沉降的可行性分析[J].測(cè)繪通報(bào)，2005(11)：10-13.

[4] 董玉森，GE L L，CHANG H C，等.基于雷達(dá)干涉測(cè)量的礦區(qū)地面沉降監(jiān)測(cè)研究[J].武漢大學(xué)學(xué)報(bào)：信息科學(xué)版，2007，32(10)：888-891.

[5] 陶秋香，劉國(guó)林，劉偉科.L和C波段雷達(dá)干涉數(shù)據(jù)礦區(qū)地面沉降監(jiān)測(cè)能力分析[J].地球物理學(xué)報(bào)，2012，55(11)：3 681-3 689.

[6] 劉霖.礦區(qū)開(kāi)采沉陷大量級(jí)形變監(jiān)測(cè)與反演分析[D].西安：長(zhǎng)安大學(xué)，2013.

[7] 陳炳乾，鄧喀中，范洪冬.基于D-InSAR技術(shù)和SVR算法的開(kāi)采沉陷監(jiān)測(cè)與預(yù)計(jì)[J].中國(guó)礦業(yè)大學(xué)學(xué)報(bào)，2014，43(5)：880-886.

[8] 郭山川，候湖平，張紹良，等.D-InSAR的黃土高原礦區(qū)地表形變監(jiān)測(cè)[J].測(cè)繪科學(xué)，2016(6)：1-11.

[9] 陳炳乾.面向礦區(qū)沉降監(jiān)測(cè)的InSAR技術(shù)及應(yīng)用研究[D].徐州：中國(guó)礦業(yè)大學(xué)，2015.

[10] 陶秋香.PS InSAR關(guān)鍵技術(shù)及其在礦區(qū)地面沉降監(jiān)測(cè)中的應(yīng)用研究[D].青島：山東科技大學(xué)，2009.

[11] 黃寶偉.基于D-InSAR和GIS技術(shù)的煤礦區(qū)地面沉降監(jiān)測(cè)研究[D].青島：中國(guó)石油大學(xué)，2011.

ApplicationofDouble-TrackD-InSARinMonitoringGroundSubsidenceofMiningArea

LIU Tong-wen1，GAO Teng-fei2，YU Guang-ting1，XU Jin-da1

(1.ShandongInstituteofSurveyingandMappingofGeology，Ji’nanShandong250002，China； 2.ShandongUniversityofScienceandTechnology，QingdaoShandong266590，China)

2017-05-16

P 237

A

1007-9394(2017)03-0007-04

劉同文(1985～)，男，山東鄄城人，碩士，注冊(cè)測(cè)繪師，主要研究方向：3S技術(shù)研究及應(yīng)用。