Hossfeld 模型在礦區地表動態沉降預測應用的可行性分析

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

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

0 引言

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

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

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

1 原理與方法

1.1 礦區開采沉陷規律

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

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

1.2 典型時間模型

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

表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為沉降參數。對式（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。進一步分析發現，參數ɑ代表時間的冪次方對沉降量的影響程度，當ɑ＞1 時，時間的影響更加顯著，w對t的響應會更快；當0＜ɑ＜1 時，時間的影響較小，w對t的響應會相對緩慢；參數b作為一個常數項，代表了時間對沉降量的影響程度，其會對公式中的分母部分進行調節，從而影響w的數值，但當參數b=0 時，w與wm相等，這表示wm是唯一影響w的變量；參數wm為最大沉降量。

為進一步探究Hossfeld 模型中不同參數對沉降量曲線的影響，模擬不同參數取值下沉降結果曲線，如圖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 中可以看出，參數ɑ越大，到達最大沉降量的時間越短，曲線彎曲程度越大，也更加陡峭；參數b越大，到達最大沉降量時間越長，曲線彎曲程度越小，也更加平緩。綜上所述，從模擬的曲線來看，Hossfeld 模型符合礦區沉降曲線前端短后端長且單調遞增的特征；且Hossfeld 模型未知參數的大小對曲線擬合程度有較大影響，需要選擇合適的參數才能準確擬合礦區開采沉降過程。

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

2 聯合水準數據的礦區開采單點動態沉降預測

2.1 研究區概況

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

圖3 研究區概況Fig.3 Overview of the study area

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

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

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

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

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

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

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

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

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

3 聯合D-InSAR 的礦區全盆地開采任意點動態沉降預測

3.1 D-InSAR 累積沉降監測

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

表4 影像參數Table 4 Parameters of the SAR imaging

圖5 研究區概況Fig.5 Overview of the study areas

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

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

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

3.2 礦區全盆地任意點預測過程分析

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

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

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

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

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

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

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

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

4 結論

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

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

2）通過D-InSAR 技術獲取大量地表累積沉降量數據，來統計 礦區Hossfeld 模型、Usher 模型和Knothe 模型任意點未知參數并分析未知參數相關性，Hossfeld 模型相干性最高；通過RMSE、RMAE 和Bland-Altman 圖三種評價指標發現，Hossfeld 模型預測礦區任意點動態沉降在精度方面具有優勢。

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