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

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

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

1 理論與方法

1.1 D-InSAR基本原理

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

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

由式（1），經(jīng)相位解纏提取地表LOS向形變相位φdeformation'，并將其轉(zhuǎn)化為SAR衛(wèi)星沿視線方向的形變值dLOS，即：

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

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

1.2 經(jīng)典GM（1，1）模型理論

式中，a,b為待估參數(shù)。

求解式（7）可得：

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

本研究提出了2種改進的GM（1，1）模型，分別為基于補償最小二乘法估計半?yún)?shù)GM（1，1）模型（BGM（1，1））和賦予相對權(quán)重的補償最小二乘法估計半?yún)?shù)GM（1，1）模型（WGM（1，1））。

1.3.1 BGM（1，1）模型

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

式（9）中加入的S也是未知量，使得BGM（1，1）變?yōu)橹忍澞Ｐ?，若按最小二乘方法求解，無法得到唯一解，所以必須添加約束條件即補償最小二乘估計準(zhǔn)則，從而得到唯一解，即求解方法由VTPV=min轉(zhuǎn)變?yōu)?/p>

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

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

1.3.2 WGM（1，1）模型

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

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

2 D-InSAR監(jiān)測分析

2.1 研究區(qū)概況及數(shù)據(jù)源

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

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

2.2 D-InSAR處理方法與監(jiān)測結(jié)果

差分干涉常有“兩過差分”、“三過差分”和“四過差分”［21］3種方法。本研究采用較為成熟的“兩過差分”法，利用監(jiān)測時間段內(nèi)相鄰時間獲取的2幅SAR影像形成干涉對和對應(yīng)的外部DEM數(shù)據(jù)，處理得到僅包含地表形變信息的干涉圖，并對得到的5期相繼時段內(nèi)地表沉降量進行時序疊加處理，最終得到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工作面下方區(qū)域，7221工作面上方村莊范圍內(nèi)未發(fā)生明顯沉降，說明所采用的注漿充填開采工藝（注漿孔位于工作面東北部）有效發(fā)揮了局域減沉減變作用。

2.3 D-InSAR監(jiān)測成果驗證與分析

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

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

3 改進的模型監(jiān)測與驗證分析

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

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

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

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

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

3.2 WGM（1，1）模型預(yù)計試驗

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

4 結(jié)語

基于D-InSAR技術(shù)獲取的沉降結(jié)果，融合GM（1，1）模型理論，實現(xiàn)了地表沉降監(jiān)測和預(yù)計的一體化。將優(yōu)化后的BGM（1，1）模型與WGM（1，1）模型應(yīng)用于礦區(qū)沉降預(yù)測可有效減小經(jīng)典GM模型誤差，提高了擬合精度和預(yù)測精度，有助于提升礦區(qū)地表沉降監(jiān)測與預(yù)計效果。