基于SBAS-InSAR技術(shù)的采空區(qū)殘余變形規(guī)律分析

梁思語，胡海峰

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

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

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

1 研究區(qū)概況

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

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

圖1 研究區(qū)位置Fig.1 Location of the study area

2 數(shù)據(jù)與方法

2.1 數(shù)據(jù)源

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

表1 Sentinel-1A衛(wèi)星影像數(shù)據(jù)Table 1 Satellite image data of Sentinel-1A

2.2 方法

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

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

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

3 可靠性分析

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

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

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

圖4 SBAS-InSAR結(jié)果與水準(zhǔn)測量結(jié)果對比Fig.4 Comparison between results of SBAS-InSAR and leveling

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

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

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

4 試驗結(jié)果及分析

4.1 沉降速率

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

(1)

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

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

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

4.2 時序累計沉降值分析

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

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

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

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

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

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

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

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

4.3 研究區(qū)單點沉降特征分析

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

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

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

5 結(jié) 論

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

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

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

3) 工作面的殘余變形經(jīng)過線性擬合，符合二維高斯曲線擬合的結(jié)果，說明工作面的殘余形變?nèi)匀环细怕史e分法。

綜上所述，SBAS-InSAR技術(shù)在研究采空區(qū)殘余變形時具有很好的應(yīng)用前景，可以有效監(jiān)測到礦區(qū)在開采結(jié)束后的地形的變化規(guī)律，為今后進行重復(fù)開采等工作提供參考。