基于SBAS-InSAR技術的采空區殘余變形規律分析

梁思語，胡海峰

(太原理工大學礦業工程學院，山西 太原 030024)

2020年，由于去產能政策調整和煤炭行業結構優化，全國關閉煤礦達428家，其中山西省關閉了32家[1]。而關閉礦山的環境恢復、改造和發展已成為礦業城市社會經濟發展的重要問題[2-3]。當礦井排水設備停止運行時，地下水位上升將改變采空區破碎巖石的應力和承載力。這導致地表發生二次變形，威脅到在封閉礦井頂部建造的建筑物的安全[3-4]。礦山關閉前后發生的形變相互關聯，隨著時間的推移，地表變形呈現出連續性。因此，獲取礦區關閉或者停止開采后發生的地表變形，對于確定地表變形規律、完善開采沉陷理論、土地利用規劃和建筑物穩定性評價具有重要意義。傳統的監測手段(如水準測量、全球衛導航衛星系統等)成本高、耗時長、地表變形的空間覆蓋度低，只能獲得點位的變形，很難獲得大范圍的地表變形的時空發展規律，而且監測點很難長期保存，受外界環境影響較大。與傳統方法相比，合成孔徑雷達干涉測量(InSAR)具有效率高、成本低、范圍大、所有天氣條件下可全天候使用等優點[5]。因此，SAR圖像對于探明礦區關閉前后地表變形的位置和范圍具有很好的應用價值。基于DInSAR技術的小基線集，InSAR技術實現了長時間、大范圍區域、連續微小形變的監測，且很大程度提高了地表形變監測的精度，得到了廣泛的應用[6]。張艷梅等[7]將多期哨兵數據采用SBAS-InSAR技術對西安市城區地表進行沉降監測，并利用已有研究資料，證明了方法的有效性；李達等[8]通過采用SBAS-InSAR技術對陜西一礦區進行了殘余變形的監測，并用DInSAR的結果驗證了SBAS-InSAR監測結果的精度；余昊等[9]用SBAS-InSAR技術對徐州西部地區關閉礦井的地表變形規律進行了研究，發現礦井關閉的地區地面有抬升的趨勢；何榮等[10]對大采深條帶開采的礦區進行地表沉降的監測研究，發現條帶工作面開采會誘發相鄰老采空區地表再次發生變形。

現有研究多針對礦區開采過程，對工作面開采結束后地表變形的研究較少，且對殘余變形的研究大多只針對某一個工作面，時間跨度小，無法反映礦井關閉后長期的變化情況。本文以山西省某礦區為研究對象，工作面的停采時間跨度為2005—2018年，選取24景Sentinel-1A影像，時間跨度為2018年1月—2018年6月、2018年10月—2019年4月，應用SBAS-InSAR技術對該礦區進行了地表形變監測，獲取了自2018年1月—2018年6月和2018年10月—2019年4月該地區地表形變信息。通過對2018年1月—6月的InSAR監測結果與實測數據進行對比分析并計算差值，發現SBAS-InSAR反演出的下沉變化情況與實際情況相符，結果表明采用SBAS-InSAR技術反演礦區殘余變形具有較好的可靠性。對2018年10月—2019年4月的反演結果進行分析，得到了其形變規律，為預測和評價采空區殘余形變提供了基礎。

1 研究區概況

研究區位于山西省長治市某礦區內，該研究區域的地表由第四系覆蓋，煤層上依次是二疊系山西組、下石盒子組上段、下石盒子組下段、上石盒子組下段。區內地形以低山-丘陵為主，地勢西高東低，為典型的北方山地丘陵地貌，屬東亞季風區半干旱大陸性氣候，四季分明，夏季多雨，春秋季多風少雨，冬季寒冷。

研究區域選取的是該礦區內已開采的9101工作面、90105工作面、90207工作面～90209工作面，上述工作面均為已經停采工作面，開采時間為2005—2018年，位置如圖1所示。本次實驗選取90208工作面、90209工作面、90105工作面以及9101工作面進行研究分析，開采煤層為9號煤層和10號煤層，開采深度約為150 m，開采平均厚度為2.4 m，煤層傾角為2°～5°。本次研究工作面的范圍跨度大，停采時間涵蓋了10 a以上，可以獲取比較完整的殘余形變規律。

圖1 研究區位置Fig.1 Location of the study area

2 數據與方法

2.1 數據源

實驗數據為24景分辨率為5 m×20 m的Sentinel-1A衛星影像數據，影像時間跨度為2018年1月—2018年6月，2018年10月—2019年4月，數據為C波段，波長為5.6 cm，具體的影像信息見表1。同時，為了減少地形相位的影響，實驗選取了分辨率為30 m×30 m的DEM數據。

表1 Sentinel-1A衛星影像數據Table 1 Satellite image data of Sentinel-1A

2.2 方法

SBAS-InSAR技術能夠有效地削弱時空失相干對時序沉降結果的影響，從而使得到的形變圖在時間和空間上更為連續[11]。 由于哨兵數據的波長是C波段(5.6 cm)， 礦區夏季植被茂密， 為避免出現嚴重的失相干現象，本次研究將數據處理分為兩次重復性實驗，避開失相干嚴重的時期。其中處于采空區下沉時期的影像做一次數據處理，獲得與水準測量時間重合的形變情況；將沉降衰退期結束后的影像做第二次數據處理，獲取到工作面的殘余形變加以分析。

首先，對影像進行預處理，將ESA提供的原始數據轉換成SLC格式的數據，提取覆蓋研究區的公共burst，并對研究區域進行裁剪。利用軌道參數、SAR影像的強度信息來配準圖像，使配準的精度達到千分之一個像元。將主影像、從影像分別去斜，再利用精密軌道數據對軌道信息進行修正。為了提高監測的精度，將時間基線閾值設置為90 d，空間基線閾值設置為150 m，通過對干涉圖濾波、相位解纏、基線精細化后，剔除掉相干性差的干涉對，分別組成30對和33對干涉對進行差分干涉處理，干涉對情況如圖2所示。為了避免低相干點引起的誤差，在差分干涉相位圖中，利用影像的相位穩定性、振幅離散指數和空間相干性等作為影像參考因子，尋找高相干的目標點。在此基礎上，對相干點進行差分、濾波、解纏[12]，通過逐次迭代，去除高程相位、大氣誤差、噪聲誤差，提取真實的形變相位，最后利用奇異值分解獲得高相干點的時間序列的沉降速率。將解得的各時段相位速率在時間域上積分，即可得到整個觀測時段的形變時間序列。

圖2 時空基線組合Fig.2 Spatiotemporal baseline combination

3 可靠性分析

為了驗證SBAS-INSAR反演結果的可靠性，本次實驗共選取工作面開采過程中布設的觀測線上的九個點用于驗證Sentinel-1A的監測結果，水準點的位置如圖3所示。由于從InSAR獲得的形變是LOS向的形變，精度分析時通過將LOS形變轉換為垂直方向的形變而不考慮水平運動來進行的。

圖3 沉降衰退期SBAS累計沉降值Fig.3 Cumulative settlement value of SBAS during settlement decline period

提取InSAR結果中對應水準點位置的沉降值[13]，InSAR處理結果與水準測量結果對比情況如圖4所示。由圖4可知，1號點、2號點、4號點～8號點的Sentinel-1A圖像產生的時間序列累計沉降趨勢和大小與水準數據均一致。雖然3號點和8號點的沉降趨勢與水準數據一致，但在大小上存在顯著差異。3號點的變形最大值達到62 mm，9號點的變形最大值達到81 mm，而InSAR監測得到的變形值只達到35.7 mm和46.9 mm。

圖4 SBAS-InSAR結果與水準測量結果對比Fig.4 Comparison between results of SBAS-InSAR and leveling

為了定量評價監測結果的準確性，將礦區的InSAR時間序列形變值線性擬合后進行插值，得到與水準數據時間重合的沉降值。采用最大偏差(MaxD)、最小偏差(MinD)、均方根誤差(RMSE)、標準差(STD)和相關系數(R2)對InSAR結果的準確性進行評價，評價結果見表2。

由表2可知，除Z2點外，所有點的InSAR結果與水準測量結果的相關系數都達到了0.9，Z2點相關性較低是由于Z2點所處位置坡度變化較大，沉降值受到了地形影響，而SBAS-InSAR監測結果由于分辨率較低，并不能準確地反映Z2點這一點位的形變值。 除Z9點外，所有點的精度均小于20 mm。所有點的平均均方根誤差和標準差分別為11.05 mm和20.22 mm，整體來看，影響InSAR技術影響監測結果精度的因素有以下幾個方面。①監測期間只采集了少量的Sentinel-1A圖像，且由于處于工作面下沉的衰退期，進行水準測量的次數較少。雖然C波段具有較大的可檢測變形梯度，且對植被和時間去相關相對不敏感，但Sentinel-1A圖像之間較長的時間間隔導致了較大的形變梯度。②進行對比分析的數據處于衰退期，沉降速率在逐漸減小，使用線性插值來比較SBAS-InSAR的結果與水準數據將不可避免地導致較低的精度水平。③每個水準點的InSAR形變值因為分辨率的原因，并不能準確的反映水準點也就是觀測點處的形變值，而是通過對水準點附近的相干點的變形進行平均得到的，這影響了對比結果。盡管仍存在偏差，但Sentinel-1A圖像監測結果的精度可以達到或優于厘米級，足以監測開采引起的地表殘余變形。

表2 SBAS-InSAR結果精度分析Table 2 Accuracy analysis of SBAS-InSAR results

4 試驗結果及分析

4.1 沉降速率

通過對影像處理，獲得了該礦區殘余形變時期2018年10月—2019年4月期間沿雷達視線方向的年平均沉降速率，將雷達視線方向的沉降速率換算到垂直方向上，得到形變速率圖如圖5所示。礦區的右側緊鄰處于開采過程中的工作面，根據開采沉陷中沉陷盆地影響半徑計算公式見式(1)。

(1)

式中：H為開采深度，m；tanβ為主要影響角正切值。

圖5 研究區域垂直沉降速率Fig.5 Vertical settlement rate in the study area

90207工作面和90208工作面的停采線方向受到其影響，無法準確獲得其沉降規律，90209工作面與其距離大于其開采沉陷半徑，未受到其影響，可以對整個工作面進行規律分析。結果顯示，2018年2月停采的90209工作面沉降速率較大，下沉速度為50～88 mm/a，2017年和2016年停采的90208工作面、90207工作面沉降速率較小，為30～59 mm/a。90105工作面開采時間為2013年，沉降速率約為10 mm/a。早期開采的9101工作面沉降速率接近0，說明地面已經基本穩定，附近沒有正在進行開采的工作面，在9號煤層和10號煤層下方的煤層也未進行開采，整體處于穩定狀態。

4.2 時序累計沉降值分析

圖6是以2018年10月28日為起始時間，其他時間相對于起始時間的形變結果。右側沉降區域的范圍在向南擴大，90209工作面未受正在開采的工作面的影響。90209工作面于2018年1月開采結束，下沉值連續6個月小于50 mm，但衰退期結束后仍在不斷下沉。由圖6可以看出，自2018年10月28日起，90209工作面出現明顯沉降，沉降區域在逐漸變大，逐漸出現一個下沉盆地，最終盆地范圍與采空區的范圍相近，沉降值最大為32 mm，從下沉量來看，90209工作面的殘余下沉仍會持續一段時間。90208工作面和90207工作面忽略鄰采的影響，靠近開切眼方向的沉降值為10 mm左右，最大的沉降值為13 mm，90105工作面沉降值較小，僅8 mm，9101工作面基本沒有沉降。通過分析可知該地質采礦條件下工作面在開采結束1 a內沉降較為明顯，在開采沉降結束1 a以上，5 a以內，仍然會有較小的沉降，開采結束5 a以上，沉降值接近于0。

圖6 研究區域時序累計沉降Fig.6 Time series cumulative settlement in the study area

圖7展示了觀測線位置。如圖7所示，選取工作面在開采過程中布設的觀測站作為時序累計沉降分析的觀測點，并延伸至下沉盆地邊緣生成東西向和南北向的剖面線，提取剖面線的時序沉降值。忽略受植被和噪聲影響的觀測點，根據90209工作面的走向剖面線和傾向剖面線上的點，繪制下沉曲線如圖8所示。

圖7 觀測線位置示意圖Fig.7 Position diagram of observation line

由圖8可以看出，90209工作面最大下沉值是32 mm，逐漸趨于穩定。從下沉曲線整體上看，在工作面內有下沉，在工作面邊界以外基本無下沉或者下沉值較小。工作面的邊界位于27號點和28號點以及58號點和59號點之間，通過走向剖面線提取的線上可以看出，在31號點和55號點，即工作面邊界附近的點出現了較大的下沉值，這是由于在開采結束后采用全部垮落法處理頂板，在工作面的邊緣會由于上方巖石的抗拉性等出現空洞，可能會由于上方覆巖體應力作用下，下方的空洞處的巖石失穩，空洞被巖石壓填，覆巖移動，造成該區域下沉較大。

圖8 下沉盆地剖面時序變化Fig.8 Time series change of subsidence basin profile

由圖8可以看出，下沉曲線保持了較好的連續性。采用二維高斯曲線擬合沉降值，并根據曲線形態分析下沉盆地的沉降特征，擬合結果如圖9所示。根據擬合結果可知，沉降值與高斯擬合曲線的相關系數為0.911 8，高斯擬合曲線與SBAS結果達到較高擬合度。該礦地表變形特征在空間上類似于高斯曲線，表明礦井地表運動特征在一定范圍內符合概率積分法模型的特征。

圖9 下沉盆地高斯擬合結果Fig.9 Gauss fitting results of subsidence basin

4.3 研究區單點沉降特征分析

為了分析研究區域的沉降性質，選取下沉盆地中心的6個點作為研究對象，其中1號點～3號點位于90209工作面內，4號點～6號點位于90208工作面內，考慮相鄰工作面開采的影響，4號點～6號點選取的位置靠近開切眼。提取出這6個點的時間序列沉降值，并用線性擬合模型擬合下沉值與下沉時間的關系。擬合后發現，90209工作面的下沉值與下沉時間呈明顯的二次相關關系，但是二次項的系數非常小，接近于線性關系。90208工作面的下沉值與下沉時間為一次線性關系，得到的線性關系式如圖10所示。

圖10 研究區域單點沉降特征Fig.10 Characteristics of single point settlement in the study area

根據多個工作面下沉值的提取發現，停采時間在1 a以內的工作面的下沉值呈二次相關關系減小，停采時間在1 a以上的工作面下沉值趨于線性形變，停采時間越長，短期內的下沉值越符合線性關系，最終下沉值都趨近于0。

5 結 論

本文采用SBAS-InSAR方法對采空區的殘余變形進行研究，通過采用SBAS-InSAR方法對24景Sentinel-1A影像分別進行兩次數據處理，分別得到了工作面開采結束后衰退期的累計沉降值，并與水準觀測得到的沉降值進行了精度對比；同時得到了地表衰退期結束以后的年平均沉降速率、時序累計沉降值等高空間分辨率的時間序列形變信息。通過分析，得到結論如下所述。

1) 將SBAS-InSAR得到的監測結果與水準測量得到的結果進行對比分析，發現二者得到的結果大小、趨勢一致，所有點的平均均方根誤差和標準差分別為11.05 mm和20.22 mm，并且與開采沉陷規律一致，說明SBAS-InSAR在進行礦區殘余變形監測方面具有較高的可靠性。

2) 研究表明，煤層采深和采厚均較小的工作面，在無外界因素干擾(如重復開采、地震等)條件下，開采結束后一年內沉降值與沉降時間符合二次相關關系，并逐漸趨于線性關系；開采結束1 a以上，沉降值與時間趨于線性關系；開采結束時間達到2 a以上，5 a以內的工作面仍然會有微小的沉降；開采結束時間在5 a以上的工作面，年沉降值接近于0。

3) 工作面的殘余變形經過線性擬合，符合二維高斯曲線擬合的結果，說明工作面的殘余形變仍然符合概率積分法。

綜上所述，SBAS-InSAR技術在研究采空區殘余變形時具有很好的應用前景，可以有效監測到礦區在開采結束后的地形的變化規律，為今后進行重復開采等工作提供參考。