Hossfeld 模型在礦區(qū)地表動態(tài)沉降預測應用的可行性分析

趙 月，王志偉，張國建，王 翔，丁文壯

（山東建筑大學測繪地理信息學院，山東 濟南 250101）

0 引言

地下煤礦開采活動打破開采工作面原有應力平衡，導致巖體發(fā)生移動變形，極易誘發(fā)地質(zhì)災害問題，威脅人民生命財產(chǎn)安全[1-2]。因此，研究礦區(qū)地表沉降規(guī)律對評估該地區(qū)潛在地質(zhì)災害風險至關(guān)重要。為準確評估煤礦開采對地表沉降的影響程度，需要先在工作面上建立地表動態(tài)沉降預測模型，獲取模型未知參數(shù)，進而分析地表任意點、任意時刻的沉降變化[3]。時間模型是目前應用較為廣泛的礦區(qū)地表動態(tài)沉降規(guī)律研究方法之一[4-7]。

礦區(qū)動態(tài)沉降過程一般分為三個階段：初始期、主要期和殘余期，沉降曲線形狀呈“S 型”增長[8]。因此，為描述礦區(qū)沉降過程多采用“S 型”時間模型，如Knothe 函數(shù)[9]、Logistic 模型[10]和Usher 模型[11]等。然而，典型的時間模型中，由于自身特性的原因，導致其形狀曲線不經(jīng)過坐標原點（即：當t=0 時，沉降量、速度和加速度不全為0）[12]，這一點不能完全符合煤礦開采沉降特征。通常，針對這一問題，可以通過修正時間零點提高預測精度，但這種修正方式具有一定的經(jīng)驗性。目前，用來開展林木生長規(guī)律研究的Hossfeld 模型[12-13]，從模型形態(tài)上來看也屬于“S 型”曲線，且函數(shù)經(jīng)過坐標原點，符合礦區(qū)沉降特征，可以嘗試用來進行礦區(qū)地表動態(tài)沉降規(guī)律研究。喬思宇等[12]基于AIC 準則評價Hossfeld 模型，分析該模型在煤礦區(qū)地表動態(tài)沉降預測中的可靠度。除此以外，在以往的研究中，學者們利用水準數(shù)據(jù)對礦區(qū)進行單點地表動態(tài)沉降預測，這種少量水準數(shù)據(jù)具有偶然性，無法證明模型適用于礦區(qū)全盆地任意點，針對這一問題，楊澤發(fā)等[14]基于InSAR 時序沉降，利用Logistic 模型分析礦區(qū)全盆地沉降時空演化規(guī)律，充分發(fā)揮InSAR 技術(shù)高分辨覆蓋率特點，探索Logistic模型參數(shù)分布規(guī)律。

本文針對典型時間模型存在時間零點問題，采用Hossfeld 模型，分別利用菏澤某礦區(qū)水準監(jiān)測數(shù)據(jù)和門克慶某礦區(qū)D-InSAR 累積沉降量數(shù)據(jù)，通過與兩種典型時間模型（Usher 模型和Knothe 模型）進行對比，分析采用Hossfeld 模型進行礦區(qū)單點和任意點地表動態(tài)沉降預測的可行性分析。一方面，在基于水準監(jiān)測數(shù)據(jù)分析時間零點對典型時間模型影響的基礎(chǔ)上，采用均方根誤差（Root Mean Square Error，RMSE）和平均絕對誤差（Mean Absolute Error，MAE）對三種時間模型單點沉降預測精度進行評價，探討不同模型單點預測的適用性；另一方面，基于DInSAR 累積沉降量數(shù)據(jù)，根據(jù)D-InSAR 獲取高分辨沉降監(jiān)測結(jié)果的優(yōu)勢，研究Usher 模型、Knothe 模型和Hossfeld 模型參數(shù)相關(guān)性，探討利用少量點實現(xiàn)全盆地任意點沉降預測的可行性，并在此基礎(chǔ)上，采用RMSE、MAE、Bland-Altman 圖分別對三種時間模型全盆地任意點地表動態(tài)沉降預測精度進行評價，對比驗證了不同模型在礦區(qū)全盆地任意點動態(tài)沉降預計中的適用性和可行性。

1 原理與方法

1.1 礦區(qū)開采沉陷規(guī)律

由煤礦開采引起的礦區(qū)地表動態(tài)沉降是一個復雜的過程，大致可以分為三個階段：初始沉降期、主要沉降期和殘余沉降期，如圖1 所示。①初始沉降期：隨著地下煤炭資源開采地表上覆地某一點會發(fā)生位置變化，此時t=0、v=0、ɑ=0、w=0，當?shù)叵麻_采影響該點時，t＞0，該點逐漸產(chǎn)生沉降，沉降量逐漸增大，速度和加速度也逐漸增加，此階段沉降量較?。虎谥饕两灯冢撼两盗侩S著時間的推移迅速增大，速度隨時間增大至峰值然后減小，加速度也隨時間增大到極值然后減小至相對應的負值，此階段地表沉降達到充分采動狀態(tài)或超充分采動狀態(tài)，沉降范圍達到最大；③殘余沉降期：速度逐漸減小至0，加速度也逐漸由負值減小至0，沉降量緩慢增加直至達到極限值趨于穩(wěn)定，此階段持續(xù)時間最長。由圖1（a）可知，沉降量滿足“S”型曲線，但曲線為非對稱性線形，曲線前端短、后端長，且單調(diào)遞增。除此之外，由于礦區(qū)不同地質(zhì)采礦條件，使煤礦開采過程存在差異，曲線表現(xiàn)為近“S”型，可通過調(diào)整參數(shù)改變曲線的陡峭程度[3]。

圖1 礦區(qū)某點沉降過程Fig.1 Subsidence process at a point in the mine

1.2 典型時間模型

作為典型“S”時間模型，Usher 模型和Knothe 模型均可以通過參數(shù)調(diào)節(jié)曲線形狀，具有較強的適應性[15-16]。然而，它們都存在t=0 時，沉降量、速度和加速度不全為0 的情況，不符合礦區(qū)沉降初始規(guī)律。為此，本文根據(jù)這兩個模型的沉降量表達式推導其對應的速度表達式和加速度表達式，并計算t=0 時，沉降量、速度和加速度特征，結(jié)果見表1。Usher 模型中，wm為最大下沉量；ɑ、b、c為沉降參數(shù)。當t=0 時，w≠0，v≠0，a≠0，不符合沉降初期規(guī)律。Knothe模型中，當t=0 時，w（t）=0，v（t）≠0，a（t）≠0，不滿足沉降初期，沉降量、速度和加速度都為0 的規(guī)律。

表1 Usher 模型和Knothe 模型沉降特征Table 1 Subsidence characteristics of Usher model and Knothe model

1.3 Hossfeld 模型

Hossfeld 模型是一種常用的理論生長模型[12-13]，可表示為式（1）。

式中：w（t）為t時刻的下沉量；wm為最大下沉量；ɑ和b為沉降參數(shù)。對式（1）進行一階求導可得速度公式，二階求導可得加速度公式，其表達式分別為式（2）和式（3）。

由式（1）～式（3）可知，當t=0 時，w（t）=0，v（t）=0，a（t）=0 ；當t=+∞ 時，w（t）=wm，v（t）=0，a（t）=0。進一步分析發(fā)現(xiàn)，參數(shù)ɑ代表時間的冪次方對沉降量的影響程度，當ɑ＞1 時，時間的影響更加顯著，w對t的響應會更快；當0＜ɑ＜1 時，時間的影響較小，w對t的響應會相對緩慢；參數(shù)b作為一個常數(shù)項，代表了時間對沉降量的影響程度，其會對公式中的分母部分進行調(diào)節(jié)，從而影響w的數(shù)值，但當參數(shù)b=0 時，w與wm相等，這表示wm是唯一影響w的變量；參數(shù)wm為最大沉降量。

為進一步探究Hossfeld 模型中不同參數(shù)對沉降量曲線的影響，模擬不同參數(shù)取值下沉降結(jié)果曲線，如圖2 所示，其中，在圖2（a）中，ɑ分別為2.6、2.8、3.0、3.2 和3.4，wm=1 000 mm，b=1 000；在圖2（b）中，b分別為50、500、1 000、2 500 和5 000，wm=1 000 mm，ɑ=3。從圖2 中可以看出，參數(shù)ɑ越大，到達最大沉降量的時間越短，曲線彎曲程度越大，也更加陡峭；參數(shù)b越大，到達最大沉降量時間越長，曲線彎曲程度越小，也更加平緩。綜上所述，從模擬的曲線來看，Hossfeld 模型符合礦區(qū)沉降曲線前端短后端長且單調(diào)遞增的特征；且Hossfeld 模型未知參數(shù)的大小對曲線擬合程度有較大影響，需要選擇合適的參數(shù)才能準確擬合礦區(qū)開采沉降過程。

圖2 Hossfeld 模型沉降參數(shù)對曲線的影響Fig.2 Influence of parameters on the subsidence curves for the Hossfeld model

2 聯(lián)合水準數(shù)據(jù)的礦區(qū)開采單點動態(tài)沉降預測

2.1 研究區(qū)概況

選取菏澤某礦為研究區(qū)（圖3）。圖3 中部圓圈標記部分為自南向北的走向線水準點，共86 個。水準數(shù)據(jù)時間跨度為2016 年12 月2 日—2018 年2 月27 日，時間周期最短間隔7 d，最長間隔43 d，總共觀測18 次。

圖3 研究區(qū)概況Fig.3 Overview of the study area

2.2 時間零點對典型時間模型沉降預測精度影響

為了說明時間零點對典型時間模型沉降預測精度影響，本文以實測水準數(shù)據(jù)為數(shù)據(jù)源，以Usher 模型和Knothe 模型為例，分別開展修正時間零點和未修正時間零點實驗，其中，修正時間零點閾值為2 cm。選取研究區(qū)范圍內(nèi)任意一點的水準觀測沉降值，采用遺傳算法（Genetic Algorithm，GA）[17]對上述兩種模型未知參數(shù)進行反演，根據(jù)反演出的參數(shù)構(gòu)建Knothe 模型和Usher 模型并進行沉降值預測，預測結(jié)果如圖4 所示。由圖4 可知，由于零點的設置不同，導致樣本點的擬合度存在差異，較多水準數(shù)據(jù)偏離擬合曲線，修正時間零點的模型明顯優(yōu)于未修正時間零點的模型。

圖4 Knothe 模型和Usher 模型和預測結(jié)果擬合度對比Fig.4 Comparison of the fitting for the predictions of the Knothe model and the Usher model

本文以RMSE 和MAE 為評價指標，分別統(tǒng)計分析全部樣本點所有預測值與實測值差值?？紤]到選用的水準數(shù)據(jù)累計沉降量部分已達2 000 mm，而對于沉降量級較大區(qū)域，預測值與實測值之間的差異通常較大，因此，選擇100 mm 為評價指標閾值，結(jié)果見表2。由表2 可知，在Usher 模型沉降預測結(jié)果中，在＜100 mm 范圍的RMSE 和MAE，未修正時間零點的Usher 模型精度低于修正時間零點的Usher 模型，經(jīng)過修正時間零點后的模型精度分別提高2.33%和3.49%。在Knothe 模型沉降預測結(jié)果中，經(jīng)過修正時間零點的模型精度分別提高4.63%和3.49%。綜上所述，對于Usher 模型和Knothe 模型而言，設置合適的時間零點會提高預測精度。然而，在面對不同地質(zhì)條件的礦區(qū)時，單純從模型公式出發(fā)，很難給出統(tǒng)一的標準來修正時間零點。因此，對于時間零點的修正問題，經(jīng)驗性因素起著至關(guān)重要的作用。

表2 修正和未修正時間零點模型預測值與水準數(shù)據(jù)實測結(jié)果對比Table 2 Comparison of predictions of corrected and uncorrected time-zero model and results of leveling data 單位：%

2.3 不同時間模型沉降預測精度對比

為了彌補典型時間模型需要修正時間零點這一缺陷，采用無需修正時間零點的Hossfeld 模型。同樣以2.2 部分水準數(shù)據(jù)為數(shù)據(jù)源，采用GA 進行未知參數(shù)反演，根據(jù)反演參數(shù)構(gòu)建Hossfeld 模型，并進行沉降值預測。通過RMSE 和MAE 這兩個評價指標探究基于Hossfeld 模型對礦區(qū)開采沉降預測結(jié)果精度，結(jié)果見表3。由表3 可知，采用Hossfeld 模型進行沉降預測結(jié)果的絕大部分樣本點誤差較小。對比表2來看，在＜100 mm 的范圍內(nèi)的RMSE，Hossfeld 模型沉降預測精度比經(jīng)過修正時間零點的Usher 模型沉降預測精度降低1.16%，比未經(jīng)過修正時間零點的沉降預測精度提高1.17%，比經(jīng)過修正時間零點和未經(jīng)過修正時間零點的Knothe 模型沉降預測精度分別提高25.28%和30.24%，說明Hossfeld 模型預測精度略低于修正時間零點的Usher 模型，高于未修正時間零點的Usher 模型，遠遠高于修正時間零點和未修正時間零點的Knothe 模型；在＜100 mm 的范圍內(nèi)的MAE，Hossfeld 模型沉降預測精度比經(jīng)過修正時間零點和未經(jīng)過修正時間零點的Usher 模型的沉降預測精度分別提高3.49%和6.98%，比經(jīng)過修正時間零點和未經(jīng)過修正時間零點的Knothe 模型的沉降預測精度分別提高33.72%和37.21%，說明Hossfeld 模型沉降預測精度均高于Usher 模型和Knothe 模型。

表3 Hossfeld 模型預測值與水準數(shù)據(jù)實測結(jié)果對比Table 3 Comparison of predictions of Hossfeld model and results of leveling data 單位：%

此外，從上述的原理與方法中，可知Usher 模型需要反演四個未知參數(shù)，而Hossfeld 模型只需要反演三個未知參數(shù)，降低了模型的復雜度；Knothe 模型雖然只有兩個未知參數(shù)，模型復雜度降低，但是模型精度有所降低。因此，相較于Usher 模型和Knothe 模型，Hossfeld 模型無論是在精度還是模型復雜度上，都具有一定優(yōu)勢。

3 聯(lián)合D-InSAR 的礦區(qū)全盆地開采任意點動態(tài)沉降預測

3.1 D-InSAR 累積沉降監(jiān)測

選取內(nèi)蒙古自治區(qū)門克慶某煤礦為研究區(qū)（圖5），圖5 圖中長方形框為Sentinel-1A 數(shù)據(jù)覆蓋范圍。選取的SAR 影像時間跨度為2017 年10 月2 日—2018 年5 月30 日，具體參數(shù)見表4。該地區(qū)地處干旱與半干旱過渡地帶，土地沙漠化和水土流失較為嚴重，植被覆蓋率較低，因此，Sentinel-1A 數(shù)據(jù)受空間失相干影響較小[18]。利用D-InSAR 技術(shù)對11 景覆蓋研究區(qū)的影像進行兩兩差分干涉處理，其中，兩兩干涉處理的干涉對，主影像為前一個時間的影像?；谲壍佬畔τ跋裣冗M行粗配準，再基于頻譜差異法精配準，通過不斷迭代，直到達到0.001 像素[19]。通過多視處理消除由單個像元散射的雷達回波信號相干疊加導致強度信息的大量噪聲[20]，并采用自適應濾波方法[21]進一步消除噪聲影響。采用最小費用流法[22]對消除噪聲后的纏繞相位進行相位解纏，并將解纏后的相位轉(zhuǎn)換為雷達視線方向的沉降量，最后通過地理編碼得到地圖坐標系下的沉降量。數(shù)據(jù)處理過程中，采用在美國國家航空航天局（National Aeronautics and Space Administration，NASA）獲得的30 m 分辨率的SRTM DEM 數(shù)據(jù)來消除地形相位的影響[23]，通過歐州航天局（European Space Agency，ESA）提供的POD 精密軌道數(shù)據(jù)（Precise Orbit Ephemerides）來消除軌道誤差帶來的誤差[24]。

表4 影像參數(shù)Table 4 Parameters of the SAR imaging

圖5 研究區(qū)概況Fig.5 Overview of the study areas

圖6 為礦區(qū)2017 年10 月2 日—2018 年5 月30日部分累積沉降結(jié)果。隨著地下煤炭不斷開采，地面沉降的量級和影響范圍都不斷擴大，造成如圖6所示的沉降盆地中心出現(xiàn)空白區(qū)域，一方面是因為礦區(qū)大沉降梯度的特點，導致礦區(qū)中心地表沉降超過Sentinel-1A 影像可監(jiān)測的沉降范圍[25]；另一方面，SAR 影像的波長與分辨率決定了可監(jiān)測最大沉降梯度，而本文選擇的Sentinel-1A 影像波長本身較短，且為去除斑點噪聲，對影像進行多視處理，導致像元分辨率降低，進一步限制可監(jiān)測最大沉降梯度的范圍[26]。為了驗證D-InSAR 技術(shù)監(jiān)測結(jié)果的可靠性，收集了圖6（c）中黑點所示的水準監(jiān)測數(shù)據(jù)，其中，走向線自東向西為MA1 號水準點～MA38 號水準點，共38 個，傾向線自北向南為MD1 號水準點～MD18 號水準點，共18 個。除空白區(qū)域外，共有14 個水準點與DInSAR 監(jiān)測結(jié)果位置相重合。圖7 為水準數(shù)據(jù)沉降量與D-InSAR 累積沉降量的差值對比圖。由圖7 可知，在MA 線上水準點與D-InSAR 結(jié)果在邊緣處有八個點重合，最大差值不超過23 mm；在MD 線上有六個點重合，其中，有三個點在邊緣處，最大差值不超過38 mm，還有三個點在中心位置重合，最大差值不超過280 mm。由此可見，礦區(qū)邊緣處差值較小，D-InSAR 精度較為可靠，而礦區(qū)中心位置差值較大，已超出D-InSAR 沉降監(jiān)測范圍，結(jié)果不可靠。后續(xù)時間模型參數(shù)求解過程中以沉降邊緣區(qū)域數(shù)據(jù)為主。

圖6 部分D-InSAR 累積沉降結(jié)果Fig.6 Partial cumulative subsidence results from D-InSAR

圖7 水準數(shù)據(jù)沉降量與D-InSAR 累積沉降量的差值對比Fig.7 Difference between the subsidence from leveling data and D-InSAR

3.2 礦區(qū)全盆地任意點預測過程分析

為利用少量地表點的監(jiān)測數(shù)據(jù)實現(xiàn)礦區(qū)全盆地任意點沉降預測，需要探究時間模型未知參數(shù)內(nèi)在規(guī)律，本文選用相關(guān)系數(shù)作為評判標準，取值范圍為-1 到1，當接近1 時，表示存在強正相關(guān)關(guān)系；當接近-1 時，表示存在強負相關(guān)關(guān)系；當接近0 時，則表示兩個變量之間幾乎沒有線性關(guān)系。通過對上述三種時間模型的全部模型參數(shù)進行統(tǒng)計，在Usher 模型中，將未知參數(shù)ɑ作為自變量，參數(shù)b作為因變量，獲取的相關(guān)系數(shù)為4.21×10-3；將未知參數(shù)ɑ作為自變，參數(shù)c作為因變量，獲取的相關(guān)系數(shù)為7.64×10-3；將未知參數(shù)b作為自變量，參數(shù)c作為因變量，獲取的相關(guān)系數(shù)為0.425，相關(guān)性低。在Hossfeld 模型中將參數(shù)ɑ作為自變量，b參數(shù)作為因變量，ɑ的取值范圍從0.730 到3.650，b的取值范圍從0.210 到1 000，獲取的相關(guān)系數(shù)為0.796。在Knothe 模型中，將參數(shù)ɑ作為自變量，b參數(shù)作為因變量，ɑ的取值范圍從106.650 到652.268，b的取值范圍從4.99×10-3到0.005，獲取的相關(guān)系數(shù)為-0.158。三種模型中，Hossfeld 模型未知參數(shù)相關(guān)系數(shù)最高，更有可能用少量地表點的監(jiān)測數(shù)據(jù)實現(xiàn)任意點的動態(tài)預計。

將D-InSAR 技術(shù)獲取的礦區(qū)累積沉降量作為數(shù)據(jù)源，利用GA 分別建立修正時間零點的Usher 模型和Knothe 模型，其中，時間零點閾值為3 cm，以及無需修正時間的Hossfeld 模型。選擇2018 年5 月30 日預測的所有樣本點統(tǒng)計Bland-Altman 圖，結(jié)果如圖8所示。由圖8 可知，絕大部分樣本點都位于上下限內(nèi)，部分點出現(xiàn)偏離（圖8 中橢圓圈出部分），主要原因為GA 算法反演參數(shù)時過早收斂，陷入局部最優(yōu)解，導致預測值存在偏差?？傮w而言，Hossfeld 模型上下限范圍明顯小于Usher 模型和Knothe 模型，說明預測得到結(jié)果差別較小，表示預測結(jié)果一致性越高。

圖8 Bland-Altman 圖Fig.8 Results of the Bland-Altman

表5～表7 分別給出了Usher 模型、Knothe 模型和Hossfeld 模型預測沉降量與D-InSAR 累積沉降量的RMSE 和MAE 分布情況，Usher 模型和Hossfeld 模型的預測樣本點精度較高，而Knothe 模型的預測樣本點精度相對較差。綜合分析發(fā)現(xiàn)，在＜20 mm 范圍內(nèi)，Hossfeld 模型較Usher 模型和Knothe 模型在RMSE 和MAE 比例中分別提升2.55%、2.55%、88.99%和75.12%。

表5 Usher 模型沉降預測結(jié)果與D-InSAR 沉降結(jié)果對比Table 5 Comparison of Usher model subsidence prediction results and D-InSAR subsidence results

表6 Knothe 模型沉降預測結(jié)果與D-InSAR 沉降結(jié)果對比Table 6 Comparison of Knothe model subsidence prediction results and D-InSAR subsidence results

表7 Hossfeld 模型沉降預測結(jié)果與D-InSAR 沉降結(jié)果對比Table 7 Comparison of Hossfeld model subsidence prediction results and D-InSAR subsidence results

綜上所述，Knothe 模型在描述礦區(qū)全盆地沉降方面的精度較低，而Usher 模型和Hossfeld 模型的預測精度相對較高。進一步對比發(fā)現(xiàn)，相比于Usher 模型，Hossfeld 模型預測礦區(qū)沉降結(jié)果精度略高。

4 結(jié)論

本文針對Usher 模型和Knothe 模型存在時間零點問題，采用Hosfeld 模型，基于水準數(shù)據(jù)和D-InSAR數(shù)據(jù)，構(gòu)建時間模型，探討不同模型在礦區(qū)開采沉降過程中預測精度，得出結(jié)論如下所述。

1）通過水準數(shù)據(jù)驗證，相較于未修正時間零點的Usher 模型和Knothe 模型，在＜100 mm 范圍內(nèi)Hossfeld 模型的RMSE 和MAE 有所提高；對于修正時間零點的Usher 模型和Knothe 模型，在＜100 mm 范圍內(nèi)Hossfeld 模型的RMSE 和MAE 高于Knothe 模型，與Usher 模型相當。

2）通過D-InSAR 技術(shù)獲取大量地表累積沉降量數(shù)據(jù)，來統(tǒng)計 礦區(qū)Hossfeld 模型、Usher 模型和Knothe 模型任意點未知參數(shù)并分析未知參數(shù)相關(guān)性，Hossfeld 模型相干性最高；通過RMSE、RMAE 和Bland-Altman 圖三種評價指標發(fā)現(xiàn)，Hossfeld 模型預測礦區(qū)任意點動態(tài)沉降在精度方面具有優(yōu)勢。

3）相較于Usher 模型，Hossfeld 模型少了一個未知參數(shù)，降低了模型復雜性的同時可以獲取相當?shù)念A測精度；相較于Knothe 模型，Hossfeld 模型多了一個未知參數(shù)，但可以提升預測精度。