基于D-InSAR與概率積分法的開采沉陷監測與預計

楊俊凱 陳炳乾 鄧喀中 何 強 趙偉穎

(1.中國礦業大學環境與測繪學院，江蘇 徐州 221116；2.礦山空間信息技術國家測繪地理信息局重點實驗室,河南 焦作 454000)

基于D-InSAR與概率積分法的開采沉陷監測與預計

楊俊凱1陳炳乾2鄧喀中1何 強1趙偉穎1

(1.中國礦業大學環境與測繪學院，江蘇 徐州 221116；2.礦山空間信息技術國家測繪地理信息局重點實驗室,河南 焦作 454000)

為掌握采動區附近高等級公路的變形趨勢，針對傳統開采沉陷監測方法存在諸多缺陷的問題，將合成孔徑雷達技術(D-InSAR)與概率積分法相結合，提出了一種礦山開采沉陷監測與預計方法。該方法首先利用D-InSAR技術對兗州某礦9308工作面進行了監測和分析，并將監測結果與路面水準測量數據進行了對比，然后基于D-InSAR的監測值，求取了9308工作面非充分采動下的概率積分法預計參數，再根據三下采煤規程對參數進行修正得到工作面全采時的預計參數，最后利用修正后的參數預計出工作面全采對附近高等級公路的影響程度。結果表明D-InSAR技術高精度的優勢可以為開采沉陷預計提供可靠的監測數據，該方法可以為實現礦山開采沉陷監測與預計一體化提供有力支持。

D-InSAR 概率積分法 開采沉陷 預計參數

煤炭是我國的主體能源，占據著我國能源消費結構中的70%左右，對于穿越高等級公路的采空區，煤層開采極易對公路造成嚴重的破壞，為了保護礦區附近重要建筑物，必須建立及時有效的開采沉陷監測與預計系統。D-InSAR技術作為近幾年新興的空間對地觀測技術之一，已經在礦區開采沉陷監測中取得了廣泛的應用，1996年Camec等人利用D-lnSAR技術監測了法國Gardanne附近煤礦的地表沉陷,證明D-lnSAR可以監測礦區地表沉陷[1]；其后國內外學者針對D-InSAR技術的不足做了許多改進,由此發展出了PS-InSAR、CPT、小基線集等技術[2-4],這些新技術往往在城市微小緩慢形變中能取得理想的效果，而在新興的礦區開采沉陷監測手段中，D-InSAR技術仍然有著不可替代的作用。而在開采沉陷預計方面，目前最常用的預計方法有概率積分法，典型曲線法和負指數函數法[5-6]，但這些預計方法都是基于傳統的測量手段獲取監測資料，傳統的測量手段雖然精度較高,卻存在工作量大、成本高、測點稀疏、人身安全不易保證等缺點[7]。針對傳統開采沉陷監測方法的諸多缺陷，為實現開采沉陷監測與預計一體化，本研究采用10景高分TerraSAR數據，利用D-InSAR技術獲取了山東兗州某礦9308工作面上方時序沉降圖，通過提取9308工作面走向和傾向2條觀測線的下沉曲線，分析了采空區上方地表移動規律，然后對2條下沉曲線進行了概率積分法擬合求參，結合三下采煤規程中預計參數與回采尺寸的關系，對預計參數進行了修正，最后利用修正后的參數對工作面全采時進行了穩態預計，預計出了工作面全采對附近高等級公路的影響程度，從而實現了礦區開采沉陷監測與預計一體化。

1 礦區開采沉陷監測與預計模型的建立

1.1 D-InSAR技術原理及監測礦區開采沉陷可行性分析

合成孔徑雷達差分干涉技術(D-InSAR)是以合成孔徑雷達復數圖像的相位信息獲取地表變化信息的技術[8]。D-InSAR技術常用的處理方法主要有二軌法和三軌法，在偏遠山區，外部DEM不易獲取，常常使用三軌D-InSAR技術，但是三軌差分法容易引入解纏相位噪聲[9]。二軌差分法是由Massonnet 等[10]提出的，它首先利用研究區域地表變化前后的2幅SAR影像生成干涉條紋圖，再用外部DEM數據模擬地形相位的條紋圖，從干涉條紋圖中去除地形信息，從而獲取地表形變信息。

D-InSAR采用重復軌道干涉測量模式，在其獲取的2幅影像干涉紋圖中，還包含了其他誤差，干涉相位主要由以下幾部分組成[11]：

(1)

式中，Φtop為地形相位，Φflt為平地相位，Φdef為視線方向(LOS)的形變相位，Φatm為大氣延遲相位，Φniose為噪聲相位，k為整周模糊度。其中Φtop和Φflt可以通過外部DEM來去除，Φatm可以用干涉圖疊加的方法來提高形變信息和大氣干擾信號之間的信噪比，降低其相位分量[12]。對于噪聲相位Φniose，一般采用多視處理和濾波抑制干涉圖的噪聲。而整周模糊度可通過相位解纏來還原真實相位，常用的解纏方法有枝切法[13]和最小費用流法[14]。最后得出了沿視線方向的地表形變相位。

在礦區開采沉陷監測應用中，D-InSAR主要受制于礦區地表形變梯度大、速度快，地表相干性低造成的失相干問題。對于地表植被引起的失相干，選取冬季、初春植被相對稀少的影像，可提高干涉對相干性。而由于形變梯度造成的空間失相干，仍然是D-InSAR技術的難點。不同的SAR衛星的臨界梯度主要受衛星的空間分辨率以及波長的影響。D-InSAR技術可監測到的臨界形變量可表示為[15-16]

(2)

式中，Wmax為地表可監測的臨界形變量；r為地表下沉主要影響半徑;g為雷達影像的地面分辨率,λ為傳感器波長。以TerraSAR衛星為例，其重訪周期為11d，若主要影響半徑r取100 m，根據式(2)可知，該衛星可監測最大沉降量為260 mm，可監測平均下沉速度為23.6 mm/d。本研究中煤層埋深620 m左右，屬于較大深部開采[17]，加之采空區上覆巖層堅硬，地表下沉量相對較小，D-InSAR技術可以在本研究實驗區內監測出完整的下沉盆地，可為概率積分法模型提供監測資料，從而實現礦區開采沉陷監測與預計集成化。

1.2 概率積分法原理

概率積分法是由我國學者劉寶琛、廖國華根據波蘭學者李特威尼申(J.Litwiniszyn)的隨機介質理論發展得來的[6]，對于半無限開采時，主斷面內任意點地表移動的變形預計公式為

(3)

2 實 驗

2.1 實驗區域及數據

為了驗證本研究所提出模型的可行性，以山東兗州某礦9308工作面為例，該工作面地表多被農田覆蓋，并且工作面正南方有2010年開采完的9306工作面和某高等級公路，9308工作面的開采必將對附近公路產生一定的影響。9308工作面設計切眼面長140m，推進長度800m，設計開采時間為2011年11月到2012年7月。煤層平均采厚5m，平均采深620m，煤層傾角為6°，開采方向自東北向西南方向，如圖1所示。

圖1 工作面位置

為了有效地實現開采沉陷監測，選取實驗區地表植被相對稀少時段的10景TerraSAR衛星影像，同時為了減小時間的失相干，將10景影像兩兩配對，共組成9組干涉對，具體影像參數見表1。

2.2 時序D-InSAR開采沉陷監測實驗

表1 干涉對組合與對應參數

為了實現實驗區域時序D-InSAR的開采沉陷監測，采用GAMMA軟件分別對9組干涉對進行二軌差分處理，此外D-InSAR易受噪聲的影響，采用自適應濾波方法多次濾波，來改善干涉條紋的質量。同時在去地形相位時選用SRTM90 m分辨率的高程數據，其垂直精度約為16 m，在二軌差分中，外部DEM10 m的垂直誤差只會導致對形變量估計約為0.9 mm的偏差，所以外部DEM的誤差可忽略[18]。然后將9組形變圖進行疊加，得到時序D-InSAR形變圖，最后經過地理編碼轉換成UTM坐標系的時序形變結果，如圖2所示。

圖2 D-InSAR時序形變圖

圖2(a)～圖2(i)所示為實驗區域2011年12月25日—2012年4月23日時序D-InSAR結果，在圖2(a)中，工作面上方未有明顯沉陷區域，表明雖然工作面起始開采時間為2011年11月，距起始監測時間工作面已經推進約100 m，但是由于工作面推進距離小于起動距(0.25～0.5)H0，所以圖2(a)中未形成明顯的下沉區域。隨著工作面的推進，有明顯的下沉盆地，并且下沉盆地影響范圍慢慢變大，直至不再增大。在整個監測時段內下沉盆地中心最大下沉值達到了682 mm，如圖2(i)，下沉盆地中心由于受工作面傾角和附近9306工作面重復采動的影響，明顯的向下山方向和9306工作面偏移。

為了定量地分析礦山開采沉陷的地表形變規律，沿9308工作面走向、傾向各布置1條觀測線，如圖3所示，然后在時序D-InSAR形變圖中分別提取出2條觀測線上的下沉值，如圖4。從圖4(a)可以看出，走向始終處于非充分采動狀態，從2011年12月25日至2012年4月2日走向累計最大下沉達到663 mm，在監測初期，下沉值普遍較小，說明未受提前開采的影響，可以不考慮其提前采動影響，下沉速度在2011-12-25—2012-01-16監測時段處于增加的狀態，在2012-01-16—2012-01-27下沉速度明顯減小，后期下沉速度出現了減小后增加，分析可能是由于受到工作面推進速度的影響。在圖4(b)中看出，在監測時段內傾向線處于非充分采動狀態，傾向整個監測時段累計最大的下沉達到682 mm，其下沉速度趨勢與走向一致，隨著工作面的開采，傾向下沉曲線的最大下沉值沿著工作面的推進方向向前移動。在走向方向下沉影響范圍已經成形，不再有明顯的變化，后期只是下沉量值的增加，如圖4(a)所示。

圖3 觀測線A、B以及路面水準點位置圖

圖4 走向A和傾向B下沉曲線

2.3 D-InSAR監測精度評定

為了驗證D-InSAR監測精度，取布設于高等級公路一側的21個保存比較好的水準點進行驗證，具體點位如圖3所示中點L01-L21。這21個水準點于2012年1月5日至2012年4月3日進行了2次觀測，其觀測時間與D-InSAR監測時段2012年1月5日至2012年4月2日較為一致，以穩定點L01為參考點，求取了其余20個點相對于L01的沉降量，具體對比結果如圖5所示。

圖5 水準與D-InSAR監測值對比

由圖5可看出，D-InSAR與水準測量值整體上保持一致，D-InSAR監測值相對于水準值的最大誤差出現在點L21，為20.6 mm，最小誤差在點L17，為0.1 mm，其他點誤差分布較均勻，中誤差為7.7 mm。經過以上驗證，證明D-InSAR可以滿足礦區開采沉陷監測的精度要求。

2.4 融合概率積分法的開采沉陷預計

為了有效地實現開采沉陷預計,選取我國應用最廣泛的概率積分法進行開采沉陷預計，每隔20 m等距提取9308工作面整個監測時段內走向和傾向曲線上的下沉信息，共計69個監測點，如圖6所示，利用開采沉陷軟件(MSPS)進行預計參數，求得9308工作面開采到2012年4月2日時的下沉系數q=0.33，主要影響角正切tanβ=2.28，水平移動系數b=0.3，最大下沉角θ=78.25°，拐點偏移距S左=-37.18 m，S右=33.30 m，S上=-0.99 m，S下=-2.53 m。

圖6 D-InSAR提取值與MSPS擬合值對比

由圖4可知，9308工作面在整個監測時段走向和傾向均處于非充分采動，經過概率積分法求參，只得到了非充分采動情況下的預計參數。在《建筑物、水體、鐵路及主要井巷煤柱留設與壓煤開采規程》中指出,當回采尺寸較小時，要得到充分采動情況的預計參數，應對非充分采動預計參數進行修正，其中水平移動系數和開采影響傳播角與回采尺寸關系不明顯，可以不用修正。采用規程中下沉系數、主要影響角正切、拐點偏移距與回采尺寸的關系并結合該礦區移動覆巖巖性進行預計參數修正，修正后的參數q=0.47， tanβ=2.3，b=0.3，θ=78.25°，S左=-41.3 m，S右=37 m，S上=-1.1 m，S下=-2.5 m。最后采用修正后的參數進行9308工作面全采時概率積分法預計，如圖7所示。

由圖7可知，當工作面全采時，采用概率積分法預計出高等級公路受采動影響長度約為982 m，最大下沉出現在L8點，最大下沉值1 017 mm，最大傾斜在L3點，達到4.2 mm/m,最大曲率0.04 mm/m2,最大拉伸變形1.8 mm/m,最大壓縮變形3.5 mm/m,傾斜和水平移動略微超過了參考文獻[19]中規定的一級公路地基移動變形容許值，已經對公路造成了輕度影響，所以建議協調好路、礦有關部門及時對路面進行定期監測和修復，以保障行車安全。

圖7 概率積分法預計的等值線圖

3 結 論

本研究利用D-InSAR技術對9308工作面進行了開采沉陷監測，并且分析了采空區上方地表移動規律，同時將D-InSAR監測值與水準實測值，對比發現，D-InSAR監測值相對于水準值的最大誤差為20.6 mm，最小誤差僅為0.1 mm。在保證D-InSAR監測精度的前提下，獲取了9308工作面非充分采動下的概率積分法預計參數，根據“三下”采煤規程對預計參數進行修正，最后預計出了工作面全采時對附近高等級公路的影響程度預計得出高等級公路最大下沉值達到1 017 mm，最大傾斜達到4.2 mm/m,最大曲率0.04 mm/m2,最大水平變形3.5 mm/m。結果證明融合D-InSAR技術和概率積分法可以有效地實現礦山開采沉陷監測與預計的集成化，對礦山安全生產具有重要意義。

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

[2] Ferretti A,Claudio P,Rocca A.Nonliner subsidence rate estimation using permanent scatters in differential SAR interferometry[J].IEEE Transactions on Geoscience and Remote Sensing，2000,38(5):202-212.

[3] Mora O，Mallorqui J J，Broquetas A.Linear and nonlinear terrain deformation maps from a reduced set of interferometric SAR images[J].IEEE Transactions on Geoscience and Remote Sensing，2003，41(10)：2243-2253.

[4] Berardino P，Fornaro Q,Lanari R，et al.A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms[J].IEEE Transactions on Geoscience and Remote Sensing,2002,40(11)：2375-2383.

[5] 鄒友峰,鄧喀中,馬偉民.礦山開采沉陷工程[M].徐州：中國礦業大學出版社，2003:69-103. Zou Youfeng，Deng Kazhong，Ma Weimin.Mining Subsidence Engineering in Mine[M].Xuzhou：China University of Mining and Technology Press，2003：69-153.

[6] 何國清,楊 倫,凌賡娣,等.礦山開采沉陷學[M].徐州:中國礦業大學出版社,1995:72-90. He Guoqing，Yang Lun，Ling Gengdi，et al.Mining Subsidence Theory[M].Xuzhou：China University of Mining and Technology Press,1995:72-90.

[7] 范洪冬,鄧喀中,祝傳廣,等.基于時序SAR技術的采空區上方高速公路變形監測及預測方法[J].煤炭學報，2012，37(11)：1841-1846. Fan Hongdong，Deng Kazhong，Zhu Chuanguang，et al.Deformation monitoring and prediction methods for expressway above goaf based on time series SAR technique[J].Journal of China Coal Society,2012，37(11)：1841-1846.

[8] 焦明連,蔣廷臣.合成孔徑雷達干涉測量理論與應用[M].北京:測繪出版社,2009. Jiao Minglian，Jiang Tingchen.Synthetic Aperture Radar Interferometry Measurement Theory and Application[M].Beijing:Surveying and Mapping Press,2009.

[9] 周建民,李 震,李新武.基于ALOS /PALSAR雷達干涉數據的中國西部山谷冰川冰流運動規律研究[J].測繪學報,2009,38(4):341-347. Zhou Jianmin，Li Zhen，Li Xinwu.Research on rules of the valley glacier motion in western China based on ALOS/PALSAR interferometry[J].Acta Geodaetica et Cartographica Sinica.2009,38(4):341-347.

[10] Massonnet D，Rossi M，Carmona C，et al.The displacement field of the Landers earthquake mapped by radar interferometry[J].Nature，1993，364：138-142.

[11] 朱建軍,邢學敏,胡 俊,等.利用InSAR技術監測礦區地表形變[J].中國有色金屬學報，2011,21(10)：2564-2576. Zhu Jianjun，Xing Xuemin，Hu Jun，et al.Monitoring of ground surface deformation in mining area with InSAR technique[J].The Chinese Journal of Nonferrous Metals,2011,21(10):2564-2576.

[12] 閆建偉,汪云甲,陳國良,等.錢營孜煤礦地表沉陷的D-InSAR監測[J].金屬礦山，2011(3):105-107. Yan Jianwei，Wang Yunjia，Chen Guoliang，et al.Ground subsidence monitoring with D-InSAR in Qianyingzi Coal Mine[J].Metal Mine，2011(3):105-107.

[13] Goldstein R M，Zebker H A，Werner C L.Satellite radar interferometry：two-dimensional phase unwrapping[J].Radio Science，1988，23(4)：713-720.

[14] Buckland J R，Huntley J M，Turner S R E.Unwrapping noisy phase maps by use of a minimum-cost-matching algorithm[J].Applied Optics，1995，34(23)：5100-5108.

[15] 武繼峰,楊志強.InSAR技術應用于礦區開采沉陷監測的數據優選研究[J].測繪通報，2013(S):92-94. Wu Jifeng，Yang Zhiqiang.Data optimization of mining subsidence monitoring by InSAR technology[J].Bulletin of Surveying and Mapping，2013(S):92-94.

[16] 姚丹丹,吳 侃,何 強.基于D-InSAR技術的煤礦區沉陷監測[J].金屬礦山，2014(11):151-155. Yao Dandan，Wu Kan，He Qiang.Coal mine subsidence monitoring based on D-InSAR technology[J].Metal Mine,2014(11):151-155.

[17] 李培現.深部開采地表沉陷規律及預測方法研究[D].徐州：中國礦業大學,2012. Li Peixian.Study on Regularity and Prediction Method of Surface Subsidence Due to Deep Mining[D].Xuzhou:China University of Mining and Technology,2012.

[18] 陳炳乾,鄧喀中,范洪冬.基于D-InSAR技術和SVR算法的開采沉陷監測與預計[J].中國礦業大學學報,2014(5):880-886. Chen Bingqian，Deng Kazhong，Fan Hongdong.Combining D-InSAR and SVR for monitoring and prediction of mining subsidence[J].Journal of China University of Mining and Technology,2014(5):880-886.

[19] GB/T7714—2011JTG/TD31 采空區公路設計與施工技術明細[S].北京：人民交通出版社，2011. GB/T7714—2011JTG/TD31 Guidelines for Design and Construction of Highway Engineering in the Mined-out Area[S].Beijing:China Communications Press,2011.

(責任編輯 石海林)

Monitoring and Prediction of Mining Subsidence based on D-InSAR and Probability Integral Method

Yang Junkai1Chen Bingqian2Deng Kazhong1He Qiang1Zhao Weiying1

(1.SchoolofEnvironmentalScienceandSpatialInformatics，ChinaUniversityofMining&Technology，Xuzhou221116，China;2.KeyLaboratoryofMineSpatialInformationTechnologies，NationalAdministrationofSurveying，MappingandGeoinformation，Jiaozuo454000，China)

In order to obtain the deformation law of high grade highway near mining area and overcome the limitations of traditional mining subsidence monitoring methods，a mining subsidence monitoring and predicting method based on the combination of the Differential SAR Interferometry (D-InSAR) technique and probability integral method was presented.The D-InSAR technology was used to monitor and analyze the subsidence trend of the 9308 working face in Yanzhou，and the monitor results were compared with leveling data.Before modifying the prediction parameters into sufficiency mining condition，the probability integral method parameters of the 9308 working face under insufficiency mining were calculated using D-InSAR monitoring results.Then the influence of high grade highway due to coal mining was predicted using the modified prediction parameters.Experimental results verified that D-InSAR technology had high monitoring accuracy which can provide reliable data for the prediction of mining subsidence，and this method could provide strong support to realize the integration of monitoring and prediction of mining subsidence.

D-InSAR，Probability integral method，Mining subsidence，Prediction parameters

2015-02-12

測繪地理信息公益性行業科研專項(編號：201412016),江蘇高校優勢學科建設工程項目(編號：SZBF2011-6-B35),礦山空間信息技術國家測繪地理信息局重點實驗室開放基金項目(編號：KLM201311)。

楊俊凱(1990—)，男，碩士研究生。

TD325，TD17

A

1001-1250(2015)-04-195-06