融合D-InSAR與GIS技術的礦區開采沉陷形變監測及預測方法

陳銀翠，徐良驥，余禮仁

(1. 安徽礦業職業技術學院，安徽 淮北 235000； 2. 安徽理工大學，安徽 淮南 232001)

1 數據處理方法和流程

1.1 D-InSAR二軌法數據處理流程

二軌法利用形變前后的兩景SAR影像差分干涉處理再導入外部DEM可得到地表形變相位[7-9]；干涉數據處理是基于ENVI SARscape平臺，其數據處理流程如圖1所示。

1.2 GIS軟件中數據處理方法及流程

D-InSAR獲得的柵格結果和GIS中矢量數據進行疊加分析，在ArcMAP中加載出D-InSAR柵格格式形變圖與礦區開采工作面、礦區界限矢量圖并顯示在一個界面范圍。在完成空間坐標系的轉換和D-InSAR柵格形變圖與GIS矢量數據融合后，對融合數據可進行沉降區的提取、剖面分析和疊加分析[11]。

2 基于D-InSAR監測值的三次樣條插值與GM(1,1)組合預測模型的建立

2.1 灰色理論預測模型

經典的GM(1,1)模型為等間距序列預測模型[12-13]。設X0為大于零的原始數據序列，將X0進行一次累加得到AGO序列X1，Z1是X1的鄰均值等權序列。即

X(0)={X(0)(1),X(0)(2),X(0)(3),…,X(0)(k)}

(1)

X(1)={X(1)(1),X(1)(2),X(1)(3),…,X(1)(k)}

(2)

其中

X(1)(n)=X(0)(1)+X(0)(2)+X(0)(3)+…+X(0)(n)

(3)

則GM(1,1)模型相應的白化微分方程為

酒店認知實踐是酒店管理專業大一學生的入學教育。課程通過學生參觀酒店、訪談酒店管理層、酒店經理人講座等具體形式來開展實踐課程。學校安排學生參觀優質校企合作單位，使學生設身處地地接觸實踐教學基地的真實對客服務場景，讓學生主動調研高端酒店的基本情況，增加同學對酒店行業的認同感。

(4)

式中，a為發展系數；b為內生控制灰數。

若令c=[a,b]T,那么c可以表示為：c=(BTB)-1BTYn

其中

t=1,2,…,n

(5)

Yn=[x(0)(2),x(0)(3),…,x(0)(n)]T

(6)

白化微分方程的解為

(7)

式(7)求出的是一次累加后的生成值，對其進行累減還原，就能夠得到預測值，即

(8)

2.2 結合三次樣條插值的非等間距GM(1,1)模型

在開采沉陷監測中，相鄰兩次監測數據時間間隔一般不相等，因而不能使用經典GM(1,1)模型對礦區沉降結果進行有效的預測分析。為了克服單一模型算法所帶來的局限性，提出將三次樣條插值與經典GM(1,1)模型相結合的模型算法[14]。模型算法的具體步驟為：

(1) 將原始影像形變圖中獲取的非等時間間隔沉降數據進行三次樣條插值。

(2) 由三次樣條插值算法獲得每個雷達重訪周期內的下沉值，然后處理生成等時間間距序列，計算級比，確定所有樣本值級比都落在區間[exp(-2/(n+1))，exp(2/(n+1))]內；如果不在，則需對數據做平滑處理。

(3) 使用經典GM(1,1)模型對處理后的數據進行計算得到白化微分方程。

(4) 求解白化微分方程得到預測序列并進行累減，對預測數據進行殘差檢驗。

2.3 精度評定

表1 后驗方差模型精度分級標準

3 D-InSAR融合GIS技術數據處理與結果分析

3.1 研究區概況及數據選擇

袁店二礦位于安徽省亳州市內，該地區地勢平坦，地面平均標高為30 m。7225工作面長630～637 m，寬144～145 m，標高為-411～-446 m。煤層平均厚度為3.81 m，煤層傾角為2～17.5°。由于工作面北部有村莊，采用注漿減沉方式進行回采。本文采用2景C波段Radarsat-2和6景Sentinel-1衛星影像數據對研究區進行分析，采用90 m分辨率的SRTM DEM數據。6個干涉對參數見表2。

表2 Radarsat-2和Sentinel-1數據干涉對參數

3.2 D-InSAR融合GIS技術干涉差分處理結果分析

采用二軌法對沉陷區范圍影像數據進行差分干涉處理獲取工作面形變圖。將6組形變圖和地表移動觀測站圖導入ArcGIS進行疊加分析，得到6個工作面的時序形變圖如圖2所示。

從圖2(a)中可知，隨著7225工作面開采工作的推進，在工作面的中間部分出現了下沉盆地。由圖2(b)—(f)可知，從2015年7月至2016年1月，7225工作面上方每一個干涉對時間段內都形成了下沉區域，而且下沉較大的區域隨著地下煤炭的開采逐漸上前推進。

4 基于三次樣條插值與GM(1,1)組合模型的礦區沉降預測分析

選取2015年4月19日至12月21日6個監測點的沉降數據為研究對象，分別為傾向線和走向線工作面上方的A14、A15、A16和B10、B11、B12，見表3。將2015年4月1日間隔24 d的下沉值作為首次監測數據，通過三次樣條插值獲得2015年4月1日至12月21日每個運行周期內下沉值。將每個周期內下沉值逐次累加得到相對于首幅影像日期的累計下沉值見表4。

將此等間距數據作為GM(1,1)模型的原始數據前9期數據(2015-04-25—2015-11-03)建模預測后兩期(2015-11-27和2015-12-21)累計下沉量，并與原始等間序列的后兩期數據進行對比，驗證預測模型的精度，結果見表5。由表5可知：到2015年11月27日的累計下沉值最大相對誤差出現在B10點處，其值為11.5%，其余都在10%之內；到2015年12月21日的累計下沉值最大相對誤差出現在B10點和B12點其值為11.2%。最大方差比為0.097，預測模型的精度能夠滿足D-InSAR定量分析的精度要求，因此可以進行沉陷預計分析。

m

表4 等間距下沉序列 m

表5 6個監測點預測值與實測值對比分析

由于7225工作面采用注漿減沉開采，因此預計地表的累計下沉量對安全開采有著重要的作用。由概率積分經驗公式可知[15]，注漿開采的下沉系數一般在0.06～0.2之間。結合7225工作面地下采煤地質條件及工作面回采情況可以計算出地表充分采動時的最大下沉量。根據下沉系數公式得到地表最大下沉值在227～756 mm之間。因此根據位于中心盆地的B12點下沉值430 mm可知，7225工作面注漿減沉效果明顯，在2015年5月到下沉穩定期間一般不會發生重大的塌陷事故。

5 結 論

(1) 本文采用融合D-InSAR和GIS技術對淮北礦業集團袁店二礦7225工作面進行分析，獲得6個工作面的時序形變圖。結果表明工作面下沉區域隨著工作面開采工作的推進逐漸上前推進，秋冬季節與夏秋季節相比相干性較好。

(2) 利用三次樣條插值與GM(1,1)組合模型，以7225工作面上6個監測點的沉降數據為研究對象進行預測。D-InSAR預測值與實測值最大相對誤差和方差比分別為11.5%和0.097，該精度能夠滿足D-InSAR定量分析的精度要求；通過下沉系數計算最大下沉值430 mm處在227～756 mm之間，一般不會發生重大的塌陷事故。