基于D-InSAR技術和改進GM（1，1）模型的礦區沉降監測與預計

石曉宇 魏祥平 楊可明 王 劍 姚樹一

（1.中國礦業大學（北京）地球科學與測繪工程學院，北京100083；2.淮北礦業股份有限公司，安徽淮北235000）

礦產資源大規模開采所引起的地面沉降會危及礦區及周邊范圍內建構（筑）物安全，甚至會影響礦區生態環境與居民生活安全，因此對礦區地表沉降進行監測預警很有必要。相較于傳統監測手段，合成孔徑雷達差分干涉測量（Differential Interferometry Synthetic Aperture Radar，D-InSAR）技術具有全天候全天時、高精度、大范圍連續空間覆蓋、成本低、零風險等優勢。該技術通過雷達干涉圖差分獲得視線向微小形變［1-2］，被廣泛應用于火山運動［3］、山體滑坡［4］、地震形變［5-6］、采礦區塌陷［7-8］以及城市地面沉降［9］等環境與災害監測中。就開采沉陷預計而言，由開采引起的地表沉陷量會隨著時空的變化而變化［10］，是一個關于多元變化的因變值。常規預計方法主要為概率積分法，但該方法無法較好地反映開采沉陷隨時間變化的動態過程，在非充分采動或極不充分采動時，預計結果往往適用性不強［11-13］。灰色模型（Grey Model，GM（1，1））以部分信息已知、部分信息未知的小數據和貧信息不確定系統為研究對象［14］，通過對已有的時間序列數據進行累加、累減等變換處理，生成具有變化規律的新時間序列，可以弱化原始時間變量的隨機性。因而可以通過建立GM（1，1）方程預測未來某一時刻的特征量，即可以有效演化地表沉降量隨時間的變化規律。已有研究表明：礦區變形監測數據具有一定的灰色特性［15-17］，然而該類研究往往是基于觀測站的水準或GPS測量結果，且與后續的演變預計環節彼此獨立，不能形成一個完整且系統化的監測與預測模型。本研究聯合D-In-SAR技術與GM（1，1）模型，建立描述沉降量與時間關系的模型方程，實現開采沉陷監測與預計的一體化，為實現礦區地面沉陷的動態監測與實時預測預警以及后續治理提供參考。

1 理論與方法

1.1 D-InSAR基本原理

D-InSAR是利用同一地區不同時相SAR影像間的差分干涉，獲取該地區的地表形變信息的一種監測技術［18］。其主要思路是借助外部DEM數據以及利用形變前后兩景SAR影像間雷達成像具有幾何微小差異的干涉信號，獲得干涉相位。通過去除其余相位的干擾影響，獲取視線（LOS）方向的地表形變信息。兩次雷達衛星成像記錄的相位差值φ可近似表示為

式中，φflat為無高程變化時平坦地形引起的相位變化；φtopography為應用干涉或外部DEM數據可消除的地形相位；φdeformation為在LOS方向上地面形變引起的相位變化；φatmosphere為大氣延遲引起的相位偏移；φnoise為噪聲相位。

由式（1），經相位解纏提取地表LOS向形變相位φdeformation'，并將其轉化為SAR衛星沿視線方向的形變值dLOS，即：

式中，k為相位整周模糊度；λ為雷達波波長。

由于垂直向形變是井下開采引起地表形變的主要貢獻分量，故忽略水平向形變對雷達視線方向形變的影響，經過解三角形可得到形變前后SAR影像間的地表垂直形變值w為

1.2 經典GM（1，1）模型理論

式中，a,b為待估參數。

求解式（7）可得：

1.3 改進的GM（1，1）模型

本研究提出了2種改進的GM（1，1）模型，分別為基于補償最小二乘法估計半參數GM（1，1）模型（BGM（1，1））和賦予相對權重的補償最小二乘法估計半參數GM（1，1）模型（WGM（1，1））。

1.3.1 BGM（1，1）模型

BGM（1，1）相較于經典GM（1，1）模型，在方程中加入了非參數分量S，較好地顧及了模型中可能存在的系統誤差，則模型方程為

式（9）中加入的S也是未知量，使得BGM（1，1）變為秩虧模型，若按最小二乘方法求解，無法得到唯一解，所以必須添加約束條件即補償最小二乘估計準則，從而得到唯一解，即求解方法由VTPV=min轉變為

根據螺旋CT成像分析和頸動脈內膜狹窄測量標準來確定患者的狹窄或閉塞程度:狹窄在30%之內為輕度狹窄,范圍在31~69%之間為中度狹窄,70%~99%為重度狹窄,100%為閉塞。

式中，α為平滑參數［20］；R為正則矩陣。

1.3.2 WGM（1，1）模型

WGM（1，1）模型是在BGM（1，1）模型的基礎上，基于平滑參數α求取方法對補償最小二乘法進行相應改進，求解方法為

WGM（1，1）模型相較于BGM（1，1）模型，將平滑參數α的范圍由[0 ,+∞]區間限定在[0 ,1]區間，在很大程度上，使平滑參數α的搜索范圍減小，從而節省了搜索時間。

2 D-InSAR監測分析

2.1 研究區概況及數據源

淮北礦區袁二煤礦位于安徽省渦陽縣境內，采區地面縱橫分布有較多人工溝渠和村莊等建筑。研究區7221工作面為俯斜開采，于2017年12月6日開始回采，2018年6月3日停采。當工作面進行頂板完全垮落式回采時，工作面上方地表將形成塌陷區并對地面公共設施和村莊等建筑造成一定程度破壞。為了精確監測地表形變，分別在地表公路或建筑體上設置了一些水準測量與變形監測點，監測地表及其重要建筑的變形情況。由于該類監測點的數量與監測范圍有限，且常常被人為破壞，無法對重要建筑物形變進行精確監測，進而無法對該類建筑物的損害情況進行有效評估。

研究區影像采用Sentinel-1A數據，外部數字高程模型選用分辨率為90 m的SRTM（Shuttle Radar Topography Mission，航天飛機雷達地形測繪使命）數據。從歐空局網站（https://scihub.copernicus.eu/）下載了2017年11月16日—2018年1月27日期間共6幅SLC數據，形成5個干涉對，干涉對影像的詳細信息如表1所示，其中VV為同向垂直極化。

2.2 D-InSAR處理方法與監測結果

差分干涉常有“兩過差分”、“三過差分”和“四過差分”［21］3種方法。本研究采用較為成熟的“兩過差分”法，利用監測時間段內相鄰時間獲取的2幅SAR影像形成干涉對和對應的外部DEM數據，處理得到僅包含地表形變信息的干涉圖，并對得到的5期相繼時段內地表沉降量進行時序疊加處理，最終得到7221工作面2017年11月16日—2018年1月27日地表沉降量的時間序列形變圖，如圖1所示。圖1（a）至圖1（e）是以2017年11月16日為參考時間得出的地表形變，可以看出，在回采前期以及回采初期（圖1（a）和圖1（b）），地表下沉量及影響范圍均較小，隨著開采時間推移，開始形成下沉盆地（圖1（c）至圖1（e）），自西南向東北方向擴張（與7221工作面推進方向一致），平面范圍逐漸增大。至2018年1月27日，最大下沉量為71 mm，位于7221工作面下方區域，7221工作面上方村莊范圍內未發生明顯沉降，說明所采用的注漿充填開采工藝（注漿孔位于工作面東北部）有效發揮了局域減沉減變作用。

2.3 D-InSAR監測成果驗證與分析

為了驗證D-InSAR技術監測成果的可靠性，從D-InSAR面狀監測結果中提取了各水準點的下沉值，與已有水準監測數據進行對比。由于水準數據監測時段為2017年10月8日—2018年1月30日，而D-InSAR研究時段為2017年11月16日—2018年1月27日，兩者時間基準不一致。為了盡可能減少因時間間隔不統一而帶來的誤差，以地表下沉速率為依據將水準數據起止日期內插至與研究時段相同的日期，即為2017年11月16日—2018年1月27日。采區7221工作面水準點布設位置如圖2所示，水準監測成果與D-InSAR監測結果對比如圖3所示。

結合圖2和圖3可知：D-InSAR監測結果與常規水準測量成果相比，下沉趨勢基本一致，D-InSAR監測的工作面各觀測線上分別在H3和I2點處（離工作面較近）下沉量最大；雖二者下沉趨勢一致，但就下沉量而言，存在一定差異，在I觀測線上表現得最為明顯，H和C觀測線上沉降趨勢則大致相同；在I觀測線上，I2、I3、I4點的水準測量下沉值與D-InSAR監測下沉值相差較大，下沉值之差最大約40 mm，推測可能是由于水準測量中人為誤差傳遞或特別時間段的失相干導致；在F觀測線上，F1點的水準數據與D-In-SAR監測數據相差最大，差值為25 mm。

3 改進的模型監測與驗證分析

3.1 改進的GM（1，1）模型擬合試驗

由D-InSAR技術獲取的7221工作面開采階段性地表沉降量，對該工作面上監測精度較高的H觀測線上的4個觀測點（H1、H2、H3、H4點）數據進行后續的建模與預測。從各形變圖中提取2017年12月10日—2018年1月27日期間各觀測點對應的像元值作為原始時間序列值，分別采用GM（1，1）、BGM（1，1）和WGM（1，1）模型進行擬合，擬合結果采用殘差絕對值均值、相對精度、后驗差等指標進行精度分析，3個模型與D-InSAR反演值的擬合精度分析見表2。

表2中，經典GM（1，1）模型加入系統誤差之后，各點的擬合精度均有一定程度提高，就殘差絕對值均值而言，減少了0.074～1.588 mm，平均相對誤差以及相對精度最大分別減少和提高了2.51%和1.69%。對3個模型擬合結果進行綜合分析，發現BGM（1，1）與WGM（1，1）模型均優于經典GM（1，1）模型；BGM（1，1）模型與WGM（1，1）模型雖然在擬合精度上相差很小，但在平滑因子α的選取方法上，WGM（1，1）模型確定了搜索上界k（k=0.5、0.8、1、5、10，步長均為0.01），因而節省了搜索時間。該模型在求解平滑參數α的過程中，最大搜索上界為1，因此可直接設置k等于1。兩者模型平滑參數求取結果對比如表3所示。

注：殘差絕對值=|實測監測值-模型擬合值|；相對誤差=殘差絕對值與實測值的百分比；關聯度表示擬合曲線與原始曲線的相似程度；相對精度=1-平均相對誤差；后驗差為殘差方差與原始數據方差之比。

根據上述分析，本研究采用WGM（1，1）模型進行沉降預測。

3.2 WGM（1，1）模型預計試驗

采用WGM（1，1）模型預測4個像元第7期（2018年2月20日）沉降數據，并采用同期的水準監測成果進行對比驗證，第7期預測值與水準數據對比見表4。由表4可知：H3和H4點第7期預測值與水準監測值基本吻合，相對誤差在4%以下，精度較高，而H1和H2點預測相對誤差最大達26.64%。原因為：利用D-In-SAR技術獲取沉降預測樣本的同時，也將沉降誤差帶入模型；時間序列周期較短，僅為12 d，監測點下沉規律復雜，6期數據中包含了沉降3個階段（初始期、活躍期、衰退期），使得原序列隨機性增強（如監測點H1），所以原序列經過模型處理后，仍然存在較大的隨機性；H2點的原序列值雖然隨機性較弱，但前4期沉降速度較快，形變至第5，6、7期階段，沉降反而變為由緩慢至穩定階段，原始序列值趨勢呈“S”型曲線狀態，而灰色模型適用于呈現指數增加的下凸曲線，致使模型方程不能準確地進行沉降預測。

4 結語

基于D-InSAR技術獲取的沉降結果，融合GM（1，1）模型理論，實現了地表沉降監測和預計的一體化。將優化后的BGM（1，1）模型與WGM（1，1）模型應用于礦區沉降預測可有效減小經典GM模型誤差，提高了擬合精度和預測精度，有助于提升礦區地表沉降監測與預計效果。