改進的殘差修正DGM(1,1)模型在煤礦地表沉降預測中的應用

袁德寶 閆若鵬, 張 玲 張力彪, 耿程興,

(1. 中國礦業(yè)大學(北京) 地球科學與測繪工程學院, 北京 100083;2. 中國自然資源航空物探遙感中心, 北京 100083)

0 引言

我國是以煤炭為主要能源的國家,煤炭儲量和采礦產(chǎn)量居世界第一[1]。煤炭的大規(guī)模開發(fā)和使用對地表土地產(chǎn)生了嚴重影響[2]。采空區(qū)引起的土地沉陷會導致下沉盆地、裂縫、臺階和塌陷坑以及露采跡地和固體廢棄物堆積壓占,給人民生命財產(chǎn)安全帶來嚴重威脅[3]。因此,礦區(qū)地表形變的長期監(jiān)測和準確預測對于安全生產(chǎn)有著重要的意義。

傳統(tǒng)的形變監(jiān)測手段如精密水準測量等都是基于點觀測的,適用于小范圍監(jiān)測,不利于開展大規(guī)模應用[4]。合成孔徑雷達干涉測量(interferometric synthetic aperture radar,InSAR)技術(shù)能夠提取地表微小形變信息。且具有受天氣影響小、覆蓋范圍大、成本低、測量精度高等特點[5-6]。

在地下采煤的過程中,會破壞周圍覆蓋巖層的應力平衡,從而導致巖層滑移引起地表形變[7]。概率積分法是目前應用最為廣泛的一種沉陷預計方法,但由于該方法需要大量、翔實的現(xiàn)場資料[8],在實際應用中工程技術(shù)人員很難對礦山進行詳細監(jiān)測,這給概率積分法參數(shù)的確定帶來了較大困難。

基于上述問題,本文選用基于小基線集(small baseline subset,SBAS)和相干點目標分析(interferometric point target analysis,IPTA)技術(shù)對內(nèi)蒙古鄂爾多斯市某礦進行地表沉降監(jiān)測。并利用殘差修正的離散灰色和基于Fisher最優(yōu)分割法改進的馬爾科夫組合模型對InSAR監(jiān)測結(jié)果進行預測。

1 研究區(qū)概況與實驗數(shù)據(jù)

1.1 研究區(qū)概況

實驗礦區(qū)位于內(nèi)蒙古自治區(qū)鄂爾多斯市境內(nèi),地理位置如圖1所示,區(qū)內(nèi)地表大部被第四系風積沙所覆蓋,植被稀疏,本文主要研究2-2上煤層06A工作面開采引起的地表沉降,工作面范圍如圖1黑框所示。

1.2 實驗數(shù)據(jù)

本文實驗數(shù)據(jù)選用28景Sentinel-1A單視復數(shù)(single look complex,SLC)影像。數(shù)字高程模型(digital elevation model,DEM)選用30 m SRTM DEM。影像獲取時間如表1所示。

圖1 研究區(qū)地理位置

表1 Sentinel-1基本信息

2 研究方法

2.1 基于SBAS的IPTA方法

2.1.1 SBAS-InSAR原理

SBAS是由Berardino[9]和Lanari[10]等人提出的,其原理為：首先對N+1景影像按時間順序排序并配準;然后選擇恰當?shù)臅r空基線依次選取N幅影像為主影像,生成M個干涉對,對干涉對進行差分干涉處理、濾波、解纏;之后對所有干涉圖組成相位方程并采用最小二乘法或者奇異值分解法進行形變參數(shù)的估計[11]。

2.1.2 IPTA技術(shù)

IPTA是瑞士GAMMA遙感公司開發(fā)的用于時序InSAR數(shù)據(jù)處理的模塊,該模塊利用雷達影像的相位和幅度信息識別大量相干點目標,并計算形變速率以及時序地表形變信息[12]。

2.2 改進的離散灰色馬爾科夫模型

2.2.1 離散灰色模型

灰色系統(tǒng)理論主要用于研究“小樣本、貧信息”的不確定性系統(tǒng)[13]。謝乃明等人進一步提出了離散灰色模型(discrete grey model,DGM),并用純指數(shù)序列驗證了其無偏性[14-15]。DGM建模過程如下：

(2)建立DGM(1,1)模型

(1)

則DGM(1,1)的最小二乘估計參數(shù)列滿足

(2)

(4)取x(1)(1)=x(0)(1),將上式所求參數(shù)a、b代入式

(3)

(5)求累減還原值

(4)

(5)

為了補償取絕對值帶來的正負號影響,在后續(xù)殘差修正過程中仍取原始殘差的正負號。

2.2.2 馬爾科夫模型

馬爾科夫建模過程為[16]：

(1)計算狀態(tài)判定的分界值,劃分狀態(tài)[17]。

(2)記由狀態(tài)i轉(zhuǎn)變?yōu)闋顟B(tài)j的概率為Pij,建立狀態(tài)轉(zhuǎn)移矩陣P(*)。

(3)根據(jù)馬爾科夫的無后向性構(gòu)建當前狀態(tài)矩陣P(0),并計算m步轉(zhuǎn)移概率

(6)

(4)記當前狀態(tài)預測值為X(0),用馬爾科夫模型對預測值進行修正,計算公式為

(7)

其中T為評價權(quán)值矩陣,通常[18]

(8)

建立馬爾科夫模型的關(guān)鍵在于狀態(tài)劃分的合理性,本文引入Fisher最優(yōu)分割法對馬爾科夫狀態(tài)進行劃分[19]。

3 實驗結(jié)果分析

3.1 IPTA-InSAR結(jié)果分析

本次實驗提取了研究區(qū)域內(nèi)2017年9月至2018年8月底沉降速率和累計沉降值,結(jié)果如圖2和圖3所示。06A工作面最大下沉量達到188.75 mm,最大下沉速率達到166.25 mm/a,且該工作面受到北部工作面影響,其下沉盆地有逐漸向西北方向移動的趨勢。

圖2 沉降速率圖

圖3 累計沉降圖

為驗證InSAR結(jié)果的準確性,本文分別選取沿傾向線和走向線均勻分布的全球定位系統(tǒng)(global positioning system,GPS)觀測點以及和GPS觀測點地理位置最接近的IPTA-InSAR相干目標點共30個進行比對,選取兩者監(jiān)測的累計沉降量作為比對值。比對結(jié)果為：走向線平均絕對誤差為18.5 mm,傾向線平均絕對誤差為17.7 mm,總體下沉趨勢基本相同。但由于GPS點通常反映的是以點為范圍的時序變化情況,而InSAR所反映的是以像元為單位的地表變化趨勢,所以二者獲取的沉降范圍、沉降速率等結(jié)果會存在一定差異。

3.2 離散灰色馬爾科夫模型預測結(jié)果分析

本文在06A工作面開采所引起的下沉盆地范圍內(nèi),均勻選取4個監(jiān)測點作為樣本點進行預測模型的驗證,樣本點位置如圖4所示。將整個研究過程按兩個衛(wèi)星重訪周期等時序分為14期,并用前10期數(shù)據(jù)建模擬合后4期的結(jié)果。

圖4 樣本點位置示意圖

3.2.1 建立離散灰色模型DGM(1,1)

根據(jù)2.2.1節(jié)提到的建模方式建立離散灰色模型,4個樣本點11到14期的原始序列、擬合結(jié)果和殘差如表2、表3所示。

表2 樣本點1、2離散灰色模型擬合結(jié)果 單位：mm

表3 樣本點3、4離散灰色模型擬合結(jié)果 單位：mm

3.2.2 建立馬爾科夫模型修正預測值

依據(jù)Fisher最優(yōu)分割法進行的分類結(jié)果和2.2.2節(jié)中所述的原理,對殘差修正后的離散灰色模型預測結(jié)果進行馬爾科夫修正,最終修正結(jié)果如表4所示。

表4 馬爾科夫模型修正預測結(jié)果 單位：mm

3.2.3 結(jié)果分析

傳統(tǒng)的離散灰色模型的預測結(jié)果并不理想,而進行馬爾科夫修正后,預測結(jié)果顯著提升。但針對第13期地面突然抬升的情況,無論哪種方法的預測結(jié)果都不夠準確,在進行精度分析時,將這個值作為不可預測的粗差剔除,具體精度對比分析結(jié)果見表5。

表5 精度對比結(jié)果 單位：mm

如表5所示,在加入馬爾科夫鏈修正后,均方根誤差分別為：1.38 mm、0.85 mm、1.37 mm、1.73 mm,相較于傳統(tǒng)模型精度有很大提高。

4 結(jié)束語

(1)該礦區(qū)在2017年9月20日至2018年8月22日,地表形變明顯,最大沉降速率達到了166.25 mm/a,最大累積沉降量為188.75 mm。

(2)采用殘差修正后的離散灰色馬爾科夫模型對礦區(qū)沉降值的預測結(jié)果要明顯優(yōu)于傳統(tǒng)離散灰色模型預測結(jié)果。