聯合InSAR和水準數據的礦區動態沉降規律分析

楊澤發，朱建軍，李志偉，汪長城，汪云甲，陳國良

?

聯合InSAR和水準數據的礦區動態沉降規律分析

楊澤發1，朱建軍1，李志偉1，汪長城1，汪云甲2，陳國良2

(1. 中南大學地球科學與信息物理學院，湖南長沙，410083；2. 中國礦業大學江蘇省資源環境信息工程重點實驗室，江蘇徐州，221116)

針對傳統礦區動態沉降測量中特征點缺失導致沉降規律不可靠問題，提出基于新型合成孔徑雷達干涉測量(interferometric synthetic aperture radar, InSAR)與傳統水準技術觀測的礦區地表沉降值，使用Logistic時間函數和Levenberg-Marquard (LM) 算法擬合分析礦區動態沉降規律以及其與地下開采的關系。研究結果表明：聯合InSAR和水準測量沉降值擬合的動態沉降函數、加速度函數的均方差比僅利用水準的均方差分別提高89.0%和97.9%；Logistic時間函數能夠較好地描述礦區地表單點的動態沉降過程，且各點的最大沉降速度和加速度極值發生時間與開采至該點的時間之間呈線性關系。

InSAR；水準；礦區動態沉降；Logistic時間函數

地下開采導致的變形預計和分析是開采沉陷研究的重要內容，也是評估其潛在地質災害的重要依據。然而，傳統的開采沉陷預計模型僅能估計地表穩定后的沉降和變形，忽略了開采過程中地表動態形變和破壞[1]，因此，研究礦區地表動態沉降規律對礦山開采損害評估及地表建構筑物防護有重要意義。目前，礦區動態沉降規律研究主要基于傳統的水準沉降數據[2]。由于該技術成本較高，其獲取的沉降數據時間分辨率通常不高，致使地表動態沉降過程拐點即“特征點”可能未被測量(即“特征點缺失”)，從而導致分析的動態沉降規律不可靠甚至錯誤[2?3]。近年來，合成孔徑雷達干涉測量(interferometric synthetic aperture radar，InSAR)技術不斷發展，其以低成本、全天候、全天時、高精度等優點在礦區地表沉降監測中得到廣泛應 用[4?6]。然而，由于當前可用SAR衛星均為重復軌道飛行，其固定的重訪周期也可能使InSAR數據未能監測動態沉降過程的特征點。考慮到單一InSAR或水準監測數據均可能出現特征點缺失，且兩者在時間上通常不同步，因此，若聯合傳統的水準與新型的InSAR監測的沉降，則將大大提高地表沉降數據的時間分辨率，增加多余觀測數，從而有效地降低單一數據源特征點缺失對動態沉降規律分析造成的不可靠程度。為此，本文作者提出聯合InSAR和水準監測沉降分析礦區地表動態規律。首先，利用模擬數據分析特征點缺失對擬合結果精度的影響，研究聯合InSAR和水準與僅利用水準監測沉降擬合的礦區地表沉降規律的差異。之后，基于錢營孜礦區水準與InSAR沉降測量值，使用Levenberg?Marquard(LM) 算法擬合Logistic時間函數的待估參數，并對比分析聯合InSAR和水準與僅利用水準數據擬合結果的差異，最后分析地表動態沉降過程特征點出現時間與地下開采之間的關系，以便為礦區合理制定地表監測方案及建筑物保護措施提供參考依據。

1 礦區地表單點的動態沉降過程與Logistic時間函數

1.1 礦區地表單點的動態沉降過程

礦區地表單點沉降是一個復雜的時空過程，大致可分為3個時期：初始沉降期、主要沉降期和殘余沉降期[7]，如圖1所示。目前，研究礦區地表動態移動規律的主要方法為時間函數法，具體有Knothe時間函 數[8?9]、雙參數時間函數[10]、廣義時間函數[11]、“S”型增長時間函數[12?13]等。其中：Knothe時間函數的一階導數(沉降速度)、二階導數(沉降加速度)與礦區實際沉降過程的對應值不符[13]，因此，其結果精度不高；雙參數時間函數和廣義時間函數待估參數較多，其實際應用受到制約；“S”型增長時間函數由于更加符合礦區地表沉降的3個時期[13]，因而被廣泛研究和使用。在礦區地表單點動態沉降過程中存在3個重要特征點，即最大沉降速度點(見圖1中點2)、最大和最小加速度點(見圖1中點1和3)，其不僅是約束沉降曲線形狀的關鍵點，而且是評估地表建構筑物破壞持續時間及最大破壞強度出現時間的主要依據。因此，研究礦區地表動態沉降規律、分析其特征點與地下開采的關系有重要意義。

圖1 礦區地表沉降過程的3個階段

1.2 Logistic時間函數

Logistic時間函數為典型的“S”型增長曲線[13]，其基本形式為

式中：()為地表在時刻的沉降值；0為地表最大沉降值；和為Logistic時間函數的形狀參數。0，和為函數的待估參數。

由式(1)分別對時間求一階、二階導數，得到沉降速度()和加速度()表達式：

由于Logistic時間函數為典型的“S”型增長曲線，其也存在與礦區單點動態沉降過程相似的3個特征點，即最大沉降速度點和2個加速度極值點。3個特征點中任意1個點缺失均會導致擬合的Logistic函數存在 誤差。

2 特征點缺失對Logistic時間函數擬合的影響

為了驗證沉降監測中特征點缺失對Logistic時間函數擬合的影響，首先假定式(1)中0=?1 m，=1 200，=0.08，=0~300 d，并利用式(1)計算模擬的動態沉降曲線，如圖2中曲線1所示。之后，假設該點有6組地表沉降測量數據，其監測周期分別為15，35，46，60，80和95 d(如圖2(a)~2(f)中圓圈點所示)，且首次觀測時間均為第1天。為了使模擬的監測數據符合實際情況，在每組監測的沉降值中加入3 cm正態隨機誤差。隨后，基于含誤差的6組監測數據，利用LM算法分別估計6組Logistic時間函數的待估參數(即0，和)，其擬合結果如圖2中曲線2所示。

監測周期/d：(a) 15；(b) 35；(c) 46；(d) 60；(e) 80；(f) 95

為了更清晰地顯示6組數據中監測時間與沉降過程特征點的分布關系，利用上述擬合的Logistic時間函數待估參數并按照式(3)分別計算6組對應的動態沉降加速度。圖3所示為模擬沉降加速度曲線(由0=?1 m，=1 200，=0.08，代入式(3)計算結果，如圖3中曲線1所示)與6組沉降加速度曲線(如圖3中曲線2所示)對比結果。

監測周期/d：(a) 15；(b) 35；(c) 46；(d) 60；(e) 80；(f) 95

為了定量比較不同周期的監測數據對Logistic函數擬合結果的影響，計算模擬和擬合的動態沉降曲線之間的均方差(mean square error，MSE)，其結果如表1所示。

表1 模擬與擬合沉降曲線均方差

從表1可以看出：由于監測周期為15 d(如圖2(a)所示)和35 d(如圖2(b)所示)的監測時間基本覆蓋了動態沉降曲線的特征點，其擬合曲線與模擬曲線吻合較好且精度相當(均方差分別為4.79 mm和7.65 mm)。然而，因為多余觀測量減少，監測周期為46 d(如圖2(c)所示)和60 d(如圖2(d)所示)的擬合曲線均方差(分別為21.72 mm和15.23 mm)均大于15 d和35 d的結果。此外，雖然監測周期分別為80 d(如圖2(e)所示)和92 d(如圖2(f)所示)的觀測次數均為4次，但后者更接近沉降過程的特征點，因此，其均方差明顯低于前者(分別為104.19 mm和39.09 mm)。另外，從表1 還可發現：監測周期為15和35 d的觀測次數(分別為21和9次)相差1倍多，但其擬合結果的精度卻非常接近(其均方差之差僅為2.86 mm)，其原因主要為后者的監測時間幾乎覆蓋了動態沉降曲線的所有特征點。

以上實驗表明：礦區地表沉降監測時間是否覆蓋該過程的特征點對擬合函數精度有較大影響，監測時間間隔太大，多余觀測不足，且易缺失特征點，使得擬合的函數精度較低。因此，聯合多源監測數據不僅能夠增加函數擬合的多余觀測量，同時增大了覆蓋動態沉降過程特征點的可能性，從而有效提高礦區地表動態沉降規律的可靠性。

3 聯合InSAR和水準數據的礦區沉降規律分析方法

鑒于傳統InSAR技術監測的地表形變為雷達視線(line of sight, LOS)方向，且礦區變形主要以沉降為 主[1]，因此，忽略水平移動分量對LOS向形變的貢獻，直接將InSAR形變LOS轉換為沉降值，即=LOS/cos(式中，為雷達傳感器的入射角)。聯合InSAR和水準數據分析礦區地表沉降規律的方法流程如下。

2) 合并InSAR與水準的監測結果，得到合并后的監測時間和沉降監測數據：

(5)

3) 從式(1)可以看出Logistic時間函數為含有指數的非線性函數，直接求解較復雜，因此，本文基于聯合的InSAR和水準沉降數據，使用LM非線性算法擬合Logistic時間函數的待估參數0，和。求出各觀測點的待估參數后，即可利用式(2)和(3)計算該點的沉降速度與加速度。

4) 估計所有監測點最大沉降速度、極大和極小加速度(3個特征點)的出現時間與地下開采經過該點的時間，分析地下開采與地表動態沉降特征點之間的關系函數，以便為合理地制定礦區變形監測方案以及地表建構筑物防護提供參考依據。

4 模擬數據驗證

為了驗證聯合InSAR和水準沉降監測數據的優勢和可靠性，假定第2節中監測周期為46 d的沉降值為InSAR監測值(圖4中三角形)，同時將之前擬合結果精度最低的一組數據(監測周期為80 d)設定為水準監測值(圖4中圓圈點)。之后，按照第3節中描述的方法擬合Logistic時間函數的待估參數，并利用擬合的參數計算沉降加速度曲線，其結果如圖4所示。

(a) 動態沉降曲線；(b) 加速度曲線

從圖4可以看出：聯合InSAR和水準監測數據擬合的動態沉降曲線(圖4(a)中曲線2)、加速度曲線(圖4(b)中曲線2)與對應的模擬曲線(圖4(a)、4(b)中曲線1)較吻合，其均方差分別為和，而僅利用水準測量擬合(見圖2(e)和3(e))的均方差分別為和。與后者相比，聯合InSAR和水準測量數據的均方差分別提高了92.68 mm和2.431 mm/d2，其提高比例(即)分別為89.0%和97.9%。另外，相對于僅利用InSAR監測數據擬合的沉降曲線(圖2(c))和加速度曲線(圖3(c))的均方差21.72 mm和0.147 mm/d2，其精度也提高了47.8%和64.6%。實驗結果表明：聯合InSAR和水準監測數據后，多余觀測量大大增加，其擬合結果的精度比僅利用InSAR或水準有大幅度提升，從而說明聯合InSAR和水準分析礦區動態沉降規律比兩者單獨使用更可靠。

5 錢營孜礦區地表動態沉降規律

5.1 錢營孜礦區InSAR監測

選取4景覆蓋安徽錢營孜礦區的ALOS PALSAR數據(幅號為660，軌道號為449)，并將時間相鄰的兩景影像組成干涉對，其干涉參數如表2所示。為了將所有影像統一到同一空間坐標系，將所有SAR影像與2010?02?28獲取的影像配準，然后對每個干涉對分別利用傳統的“二軌法”(two-pass)獲取礦區地表LOS向形變。其中，利用3弧秒shuttle radar topography mission (SRTM) digital elevation model (DEM)去除干涉對中的地形相位，使用改進的Goldstein濾波[14]對干涉圖濾波。由于礦區范圍較小，且大氣擾動在1 km左右的距離內是相關的[15]，因此，大氣擾動引起的誤差可以忽略不計。此外，由于軌道誤差和大氣相位長波部分在干涉圖中通常以線性趨勢呈現[16?17]，因此，可利用多項式擬合削弱其對形變結果的影響。最后，利用最小費用流法解纏差分干涉圖并將相位轉換為LOS向形變值，得到3個時間段內的LOS向差分形變。為了獲得礦區累計LOS向形變，將干涉對中相干性均高于0.3的點的形變值累加，然后忽略水平移動貢獻而直接將LOS向形變轉換為沉降值。

表2 干涉對參數

5.2 錢營孜礦區地表動態沉降估計

錢營孜礦區地表布設有一條走向觀測線，在地下工作面開采期間共進行了10次觀測(觀測時間如圖6圓圈點所示)。為了減少內插引入的誤差，選取同時有水準和InSAR監測數據的6個觀測點(如圖5所示)，分析該礦區地表動態沉降規律。之后，按照第3節中描述的方法分別獲得僅利用水準與聯合InSAR和水準監測數據擬合的Logistic函數待估參數，其結果分別如圖6中曲線1和2所示。

2009?01?01和2009?03?01為觀測日期

從圖6可以看出：無論對于水準數據或聯合InSAR和水準數據，Logistic時間函數均能較好地擬合監測數據。此外，圖6中點I~III(圖6(a)~(c))由于水準測量幾乎沒有缺失沉降特征點，所以，聯合InSAR和水準與僅利用水準數據擬合的曲線差異不大。而對于點Ⅳ~Ⅵ(圖6(d)~(f))，僅基于水準測量數據擬合的動態沉降曲線相對于聯合InSAR和水準測量值的結果均出現了主要沉降期時間增長、起始和終止沉降時刻滯后、最大沉降速度減小等問題。其主要原因是水準數據在主要沉降期(見圖1)缺少觀測值，從而無法約束該過程。若僅利用水準數據擬合的規律評估地表動態破壞和潛在地質災害，容易使結果偏小，從而增加地表建構筑物的潛在危險。

觀測點(見圖5)：(a) I；(b) II；(c) III；(d) IV；(e) V；(f) VI

1—僅利用水準擬合的Logistic曲線；2—聯合InSAR與水準數據擬合的Logistic曲線；

○—水準監測數據；△—InSAR監測數據

圖6 2009?11?19水準測量數據與聯合InSAR和水準測量數據擬合6個觀測點的沉降曲線

Fig. 6 Fitted Logistic curves of six observation points based on leveling measurements and integration of InSAR and leveling measurements

對于聯合InSAR和水準監測值后的點Ⅳ~Ⅵ(圖6(d)~(f)曲線2)，由于增加了2010?05?31的InSAR監測數據，其對Logistic函數的整體形狀起到了較好的約束作用(特別是對主要沉降期)，因此，聯合InSAR和水準數據比僅利用水準得出的規律更可靠。實驗結果表明：Logistic模型能夠較好地描述礦區地表動態沉降過程，而聯合InSAR和水準數據分析礦區動態沉降規律由于增加了額外的監測時間和觀測值，其結果比僅利用水準擬合結果可靠。

5.3 錢營孜礦區地下開采與地表動態沉降的關系

在礦區地表動態沉降過程中，最大沉降速度出現的時間是地表變形最迅速的時間，此時，地表建構筑物破壞最劇烈；而沉降加速度的2個極值點不僅是沉降從初始沉降到加速沉降、從加速沉降進入殘余沉降的2個重要時間點，也是地表建構筑物經歷潛在風險的起始和結束時間點。因此，研究地下開采與最大沉降速度、加速度極值出現時間的關系，對于提前制定地表建筑物防護措施和地表監測方案具有重要意義。

為了分析該礦區不同監測點的動態沉降過程特征點與開采進度的關系，首先分別計算該工作面(如圖5)所示)推進至6個觀測點時的時間1。隨后，基于聯合InSAR和水準數據擬合的Logistic模型待估參數，利用式(2)和(3)估計6個監測點的最大沉降速度發生時間2、加速度極大值、加速度極小值發生時間3和4，其結果如表3所示。最后，分別擬合1與2以及3和4之間的關系，其結果如圖7所示。

表3 開采至觀測點、最大沉降速度及加速度極值出現時間

(a) 最大沉降速度出現時間(T2)；(b) 加速度曲線極大值出現時間(T3)；(c) 加速度曲線極小值出現時間(T4)

從圖7可以看出：開采至觀測點的時間(1)與最大沉降速度出現時間(2)(如圖7(a)所示)、與加速度2個極值點出現時間(3和4)(分別如圖7(b)和7(c)所示)均呈線性函數關系，其相關系數為0.92~0.99。根據上述關系和地下工作面開采進度即可提前估計出地表建構筑物破壞持續時間(3到4時間段)、最劇烈的時間(2)，合理地對建構筑物采取變形防護措施和有效地制定地表監測方案。

6 結論

1) 鑒于單一數據源易導致動態沉降過程特征點缺失，從而降低沉降規律的可靠性問題，提出了聯合InSAR和水準監測數據分析礦區地表動態沉降規律的方法。首先通過模擬數據研究了特征點缺失對Logistic時間函數擬合精度的影響，分析得出聯合InSAR和水準比僅用水準數據擬合的沉降規律更可靠。

2) 監測數據是否覆蓋動態沉降過程對擬合結果有較大影響，而聯合InSAR和水準監測數據擬合的Logistic時間函數、加速度函數的均方差相對于僅利用水準數據的均方差分別提高了89.0%和97.9%。

3) Logistic時間函數能較好地描述了礦區地表動態沉降過程，地下開采至觀測點的時間與最大速度、加速度極值發生時間之間滿足線性函數關系。該關系將地下開采與地表動態沉降的特征點聯系起來，可為更好地制定地表監測方案、確定合理的地表建構筑物防護措施提供參考依據。

[1] Kratzsch H. Mining subsidence engineering[M]. Berlin: Springer, 1983: 1?20.

[2] Karmis M, Agioutantis Z. Enhancing mine subsidence prediction and control methodologies for long-term landscape stability[M]. United States, OSMRE: Pittsburgh Center, 2009: 12?31.

[3] 閻躍觀, 戴華陽, GE Linlin, 等. DInSAR動態沉降監測錯失點問題研究[J]. 煤炭學報, 2012, 37(12): 2038?2042. YAN Yueguan, DAI Huayang, GE Linlin, et al. Problem of dynamic of monitoring subsidence feature points missed by DInSAR technology[J]. Journal of China Coal Society, 2012, 37(12): 2038?2042.

[4] Jiang L M, Lin H, Ma J W, et al. Potential of small-baseline SAR interferometry for monitoring land subsidence related to underground coal fires: Wuda (Northern China) case study[J]. Remote Sensing of Environment, 2011, 115: 257?268.

[5] 尹宏杰, 朱建軍, 李志偉, 等. 基于SBAS的礦區形變監測研究[J]. 測繪學報, 2011, 41(1): 52?58.YIN Hongjie, ZHU Jianjun, LI Zhiwei, et al. Ground subsidence monitoring in mining area using DInSAR SBAS algorithm[J]. Acta Geodaetica et Cartographica Sinica, 2011, 41(1): 52?58.

[6] Miguel C, Andrew J, Hooper R. Surface deformation induced by water influx in the abandoned coal mines in Limburg, the Netherlands observed by satellite radar interferometry[J]. Journal Applied Geophysics, 2013, 88: 1?11.

[7] Knothe S. Time influence on a formation of a subsidence surface[J]. Archiwum Gornicwai HuticwaKrakow, 1952, 1(1): 1?3.

[8] Lbrahim D, Yasuhiro M, Hiro L. GIS-Based computational method simulating the components of 3D dynamic ground subsidence during the process of undermining[J]. International Journal of Geomechanics, 2012, 12: 43?53.

[9] Kwinta A, Hejmanowski R, Sroka A. A time function analysis used for the prediction of rock mass subsidence[C]//Proceeding of the International Symposium on Mining Science and Technology. Xuzhou, China, 1996: 419?421.

[10] Komalski A. Surface subsidence and rate of its increments based on measurements and theory[J]. Archives of Mining Science, 2001, 46(4): 391?406.

[11] 王正帥, 鄧喀中. 采動區地表動態沉降預測的Richards模型[J]. 巖土力學, 2011, 32(6): 1664?1668.WANG Zhengshuai, DENG Kazhong. Richards model of surface dynamic subsidence prediction in mining area[J]. Rock and Soil Mechanics, 2011, 32(6): 1664?1668.

[12] 張文志, 鄒友峰, 任筱芳. Logistic模型在開采沉陷單點預測中的研究[J]. 采礦與安全工程學報, 2009, 26(4): 486?489. ZHANG Wenzhi, ZOU Youfeng, REN Xiaofang. Research on logistic model in forecasting subsidence single-point during mining[J]. Journal of Mining & Safety Engineering, 2009, 26(4): 486?489.

[13] Li Z W, Ding X L, Zheng D W, et al. Least squares-based filter for remote sensing image noise reduction[J]. IEEE Transactions on Geoscience and Remote Sensing, 2008, 46(7): 2044?2049.

[14] LI Zhiwei, DING Xiaoli, CHEN Dawei, et al. Assessment of atmospheric effects on repeat-pass insar measurements in southern China region[J]. Journal of Geospatial Engineering, 2005, 7(1): 1?10.

[15] Tarayre H, Massonnet D. Atmosphereic propagation heterogenetities revealed by ERS-1 interferometry[J]. Geophysical Research Letters, 1996, 23(9): 989?992.

[16] LI Zhiwei, DING Xiaoli, HUANG Chen, et al. Modeling of atmospheric effects on InSAR measurements by incorporating terrain elevation information[J]. Journal of Atmospheric and Solar-Terrestrial Physics, 2006, 68(1): 1189?1194.

Analysis of law of kinematic mining subsidence by integrating InSAR and leveling measurements

YANG Zefa1, ZHU Jianjun1, LI Zhiwei1, WANG Changcheng1, WANG Yunjia2, CHEN Guoliang2

(1. School of Geosciences and Info-Physics, Central South University, Changsha410083, China; 2. Key Laboratory of Resources and Information Engineering of Jiangsu Province, China University of Mining and Technology, Xuzhou 221116, China)

To reduce the unreliability of the law of kinematic mining subsidence derived from leveling measurements with temporal feature points loss, the integrated new interferometric synthetic aperture radar (InSAR) and traditional leveling measurements with Logistic time function and Levenberg-Marquard (LM) algorithm were proposed, and the law of kinematic mining subsidence and the relationship between the law and underground extraction were analyzed. The results show that the mean square errors (MSE) of fitted subsidence and acceleration using InSAR and leveling measurement are 89.0% and 97.9% higher, respectively, than that only using leveling observations with temporal feature points loss. Logistic model can fit the kinematic subsidence of a ground point greatly well, and there is linear relationship between the occuring time of the maximum velocity and acceleration extremum of fitted Logistic time function and that of advancing working panel under the observation point, respectively.

InSAR; leveling; kinematic mining subsidence; Logistic time function

10.11817/j.issn.1672-7207.2015.10.026

TP237

A

1672?7207(2015)10?3743?09

2014?10?10；

2014?12?22

國家自然科學基金資助項目(41222027，40974006)；國家高技術研究發展計劃(863計劃)項目(2012AA121301)；湖南省杰出青年科學基金資助項目(13JJ1006)；國土資源部公益基金資助項目(201211011)；中央高校基本科研業務費專項資金資助項目(2013zzts249)(Projects (41222027, 40974006)supported by the National Natural Science Foundation of China; Project(2012AA121301)supported by the National High Technology Research and Development Program(863 Program) of China; Project(13JJ1006)supported by Outstanding Young Science Foundation of Hunan Province; Project(201211011) supported by the Public Welfare Fund of Ministry of Land and Resources;Project(2013zzts249) supported bythe Fundamental Research Funds of Central Universities)

朱建軍，教授，從事InSAR、測量平差及數據處理研究；E-mail：zjj@mail.csu.edu.cn

(編輯 陳燦華)