基于InSAR 的煤礦采空區(qū)地表形變監(jiān)測與預(yù)警

王鳳云，陶秋香，陳 洋，韓 宇，郭在潔

（山東科技大學 測繪與空間信息學院，山東 青島 266590）

近年來，能源需求不斷提升，煤炭開采規(guī)模逐漸擴大，經(jīng)過大幅度煤層開采后，地下會形成采空區(qū)，容易引發(fā)采空區(qū)的地表形變，即地表下沉、裂縫和地表抬升等現(xiàn)象，損壞農(nóng)用耕地、建筑群等基礎(chǔ)性設(shè)施；礦震是一種非天然地震，是由于地下持續(xù)開礦挖井，形成大面積空洞，導(dǎo)致采空區(qū)上方性質(zhì)不同的層狀巖層有規(guī)律性地彎曲、離層、斷裂下沉，從而引發(fā)巖層內(nèi)部儲存的彈性能釋放的現(xiàn)象，容易引起大面積地表塌陷，危害人民的生命財產(chǎn)安全[1-3]。2011 年11 月3 日，河南義馬煤業(yè)千秋煤礦因礦震引發(fā)重大沖擊壓事故；2016 年4 月19 日，北京房山發(fā)生2.7級礦震；2019 年6 月9 日，吉林長春龍家堡礦業(yè)發(fā)生2.3 級礦震，導(dǎo)致井下當班作業(yè)人員被困；2020 年12 月15 日，陜西榆林市榆陽區(qū)發(fā)生2.6 級礦震，導(dǎo)致煤礦停產(chǎn)。因此，在提高礦區(qū)生產(chǎn)效率的同時，保證煤礦安全開采，迫切需要對礦區(qū)地表形變進行有效監(jiān)測和穩(wěn)定性分析，建立礦震預(yù)警機制，預(yù)防潛在礦震等災(zāi)害的影響[4-5]。

合成孔徑雷達干涉測量InSAR（Interfero -metry Synthetic Aperture Radar）技術(shù)因其擁有全天候、全天時、高精度、高空間分辨率及覆蓋面廣等優(yōu)點逐步用于礦區(qū)的地表形變監(jiān)測[6]。隨著該技術(shù)的不斷完善，繁衍出DInSAR（Differential InSAR）、SBAS In－SAR（Small BAseline Subset InSAR）等多種地表形變監(jiān)測手段。許多研究人員將InSAR 技術(shù)應(yīng)用于采空區(qū)地表形變監(jiān)測，姚佳明等[7]選用升、降軌L 波段PALSAR-2 影像數(shù)據(jù)，利用InSAR 技術(shù)對煤礦采空區(qū)開展了短期動態(tài)地表沉降監(jiān)測，結(jié)合研究區(qū)開采信息對煤礦采空范圍及開采時間進行反演，驗證了InSAR技術(shù)對煤礦采區(qū)反演的合理性與可靠性；何榮興等[8]介紹了采空區(qū)災(zāi)害類型及不同類型的相應(yīng)案例，分析了采空區(qū)災(zāi)害發(fā)生的一般規(guī)律和特征，提出了盡量選擇不產(chǎn)生采空區(qū)的采礦方法；栗明明等[9]利用39景Sentinel-1 降軌數(shù)據(jù)，采用SBAS InSAR技術(shù)獲取隧道周邊采空區(qū)地表形變發(fā)展過程，證實SBAS In－SAR 技術(shù)能夠獲得采空區(qū)的毫米級沉降結(jié)果。

為深入探究InSAR 技術(shù)在采空區(qū)地表形變的監(jiān)測能力及對礦震的預(yù)警能力，驗證該技術(shù)的監(jiān)測精度，以近期發(fā)生礦震的山東某煤礦為研究區(qū)，依據(jù)礦震發(fā)生時間，選取2020-09-01—2020-12-31 覆蓋該煤礦的11 景C 波段Sentinel-1A SAR 影像，分別采用DInSAR、SBAS InSAR 技術(shù)獲取該時段內(nèi)的礦區(qū)形變監(jiān)測結(jié)果，并結(jié)合實際進一步分析時間演變生成的累計形變。

1 研究區(qū)與數(shù)據(jù)源

1.1 研究區(qū)概況

山東濟寧市境內(nèi)的魯西煤炭生產(chǎn)基地是全國煤炭生產(chǎn)基地之一，礦產(chǎn)資源豐富，經(jīng)過長期持續(xù)高強度開采，導(dǎo)致大面積土地沉陷，礦震多發(fā)頻發(fā)，生態(tài)環(huán)境遭到嚴重破壞。因此，為保護耕地及礦區(qū)生態(tài)環(huán)境，減少由于煤礦采空區(qū)巖層的移動變形導(dǎo)致的地表塌陷以及地表建筑損傷倒塌現(xiàn)象，避免或減輕礦震災(zāi)害，迫切需要有效監(jiān)測礦區(qū)采空區(qū)地表形變情況，查清采空區(qū)沉陷現(xiàn)狀，加大對采煤塌陷地區(qū)的監(jiān)測監(jiān)管力度[10]。

2020 年12 月23 日，山東省濟寧市曲阜（35.54°N，116.92°E）發(fā)生M2.4 級礦震，震中位于已停采的濟寧某煤礦采空區(qū)。此次礦震，無開采工作面，井下安全，地表無塌陷，無地表建筑物和人員財產(chǎn)損失，各項情況正常。該煤礦位于兗州市以東約15 km，曲阜市西南約10 km，陵城鎮(zhèn)附近。

1.2 SAR 衛(wèi)星數(shù)據(jù)

Sentinel-1 衛(wèi)星，載有C 波段合成孔徑雷達，具備多種工作模式，是歐洲航天局哥白尼計劃（Global Monitoring for Environment and Security，GMES）中的地球觀測衛(wèi)星。該衛(wèi)星實現(xiàn)單、雙極化等若干種不同的極化方式，能提供連續(xù)、全天候的雷達影像，擁有高重訪頻率、高覆蓋能力以及極好的時效性和可靠性[11]。選取了覆蓋礦區(qū)的11 景C 波段、VV 極化的Sentinel-1A 升軌數(shù)據(jù)，時間跨度為2020-09-01—2020-12-30。SAR 影像的具體參數(shù)見表1。

表1 SAR 影像的具體參數(shù)Table 1 Specific parameters of SAR images

為滿足本次研究，去除地形相位，選用了由美國太空總署（National Aeronautics and Space Administration, NASA）和國防部國家測繪局（National Imagery and Mapping Agency, NIMA）聯(lián)合測量的地面分辨率為90 m、平均精度16 m 的SRTM3-DEM[12-13]。

2 InSAR 技術(shù)基本原理

2.1 DInSAR 原理

DInSAR 技術(shù)是通過對同一地區(qū)不同時間的2 幅SAR 影像進行差分干涉處理獲取干涉相位，利用多圖像重復(fù)干涉或者引入外部DEM 模擬地形信息去除地形相位，從而獲取地表微量形變的測量技術(shù)[14-15]。DInSAR 獲取的相位可以表示為[13]：

對相位φdef進行相位解纏處理，得到真實形變相位φreal，然后提取雷達視線向上的地表形變量△Rtow，表示為[16]：

2.2 SBAS InSAR 原理

假設(shè)雷達傳感器在同一研究區(qū)不同時刻獲取n+1 幅SAR 影像，通過給定時間和空間閾值，生成m 幅差分干涉圖。用ta、tb（ta>tb）時刻獲取的SAR 影像生成去除地形相位的第k 幅差分干涉圖，則第k幅影像在方位-距離像元坐標系（x，r）中的差分干涉相位δφk（x，r）為：

式中：φ（tb，x，r）、φ（ta，x，r）為tb、ta時刻相位；d（ta，x，r）、d（tb，x，r）為視線向累計形變量。

去除地形殘余相位、大氣延遲相位以及各種噪聲相位后，地表形變的平均地表形變速率νT為：

則相位為：

式中：Ej、Sj為主、從影像獲取時間；νk為k 時刻像元形變速率。

由此，定義A（j，k）=tk-tk-1，且A 是1 個m×n 的秩虧矩陣，得到矩陣方程：

通過奇異值分解法和最小二乘法可求出平均形變速率相位值，得到地表累計線性形變量[17-18]。

2.3 數(shù)據(jù)處理主要步驟

數(shù)據(jù)處理流程圖如圖1。

圖1 數(shù)據(jù)處理流程圖Fig.1 The primary data processing procedure

1）DInSAR 。DInSAR 基于干涉相位獲取地表形變信息，關(guān)鍵技術(shù)主要包括：主、輔影像預(yù)處理、影像配準及重采樣、干涉圖生成、基線估計、地形相位去除、差分干涉圖濾波、相位解纏、地理編碼等，通過相位-形變轉(zhuǎn)換，得到研究區(qū)的地表形變信息。

2）SBAS InSAR。選擇2020 年9 月13 日的影像為主影像，其余為輔影像，建立連接圖，設(shè)置合適的時間基線和空間基線閾值，將滿足時空基線閾值條件的2 幅影像進行差分干涉處理，得到時序差分干涉圖，影像濾波后選取高相干像元，并進行相位解纏，選擇無殘余地形條紋且遠離形變區(qū)的地面控制點進行相位修正去除相位偏移，估算形變速率和殘余地形，利用二次解纏優(yōu)化干涉圖，進行大氣濾波估算，去除大氣相位，得到時間序列上最終位移結(jié)果。

3）對比DInSAR 和SBAS InSAR 得到的各成像時刻的形變量，研究2 種InSAR 技術(shù)對采空區(qū)地表形變的監(jiān)測和礦震預(yù)警能力。

3 實驗結(jié)果分析

3.1 DInSAR 形變監(jiān)測結(jié)果

結(jié)合Sentinel-1A 影像的獲取時間，以12 d 為1個監(jiān)測時段，利用該技術(shù)獲取各時段礦區(qū)地表形變信息及累計形變信息。

3.1.1 特征點DInSAR影像疊加圖分析

根據(jù)DInSAR 監(jiān)測到的地表形變分布，選擇形變特征明顯并且能夠監(jiān)測到形變數(shù)據(jù)的9 個特征點進行數(shù)據(jù)提取與研究分析。9 個特征點分布與DInSAR各時段地表形變影像疊加圖如圖2，9 個特征點分布與DInSAR 累計地表形變影像疊加圖如圖3。

由圖2、圖3 可以看出：

圖2 9 個特征點分布與DInSAR 各時段地表形變影像疊加圖Fig.2 Overlaying charts of 9 characteristic points and DInSAR subsidence images at each period

圖3 9 個特征點分布與DInSAR 累計地表形變影像疊加圖Fig.3 Overlaying charts of 9 characteristic points and cumulative subsidence images by DInSAR at each period

1）DInSAR 監(jiān)測到的地表形變中心（35.533°N，116.918°E）與官方發(fā)布的礦震中心“D”（35.54°N，116.92°E）相距約417 m。官方給出的礦震中心位置僅保留了小數(shù)點后2 位，是1 個大致位置，且“D”處于DInSAR 監(jiān)測的地表形變范圍之內(nèi)，2020-09—2020-12 期間，“D”處地表經(jīng)歷了“抬升-下沉-抬升”。因此，DInSAR 監(jiān)測到的地表形變中心在官方發(fā)布的礦震中心位置的誤差范圍內(nèi)。

2）在第1（2020－09－01—2020－09－13）、第4（2020－10－07—2020－10－19）、第6（2020－10－31—2020－11－12）3 個時段內(nèi)礦區(qū)上覆地表存在大面積抬升，在第2（2020－09－13—2020－09－25）、第3（2020－09－25—2020－10－07）、第5（2020－10－19—2020－10－31）3 個時段緩慢下沉。在第7 監(jiān)測時段（2020－11－12—2020-11-24）內(nèi)，地表出現(xiàn)-8～-10 mm 之間的不規(guī)則的大面積下沉，監(jiān)測區(qū)域內(nèi)最大形變量達到-11.5 mm。在第8 監(jiān)測時段（2020-11-24—2020-12-06），礦區(qū)中心下沉，周邊部分區(qū)域發(fā)生抬升，第9 監(jiān)測時段（2020-12-06— 2020-12-18）礦區(qū)中心仍持續(xù)下沉，周邊區(qū)域大面積抬升，地表高度差逐漸增大，地表結(jié)構(gòu)變形，導(dǎo)致在第10 個監(jiān)測時段（2020-12-18—2020-12-30）內(nèi)發(fā)生礦震。

3.1.2 特征點DInSAR 形變數(shù)據(jù)分析

9 個DInSAR 特征點形變數(shù)據(jù)見表2、表3，9 個DInSAR 特征點形變數(shù)據(jù)圖如圖4。

表2 9 個DInSAR 特征點各時段形變量Table 2 Settlement data at each period of 9 characteristic points by DInSAR

表3 9 個DInSAR 特征點累計形變量Table 3 Cumulative settlement data of 9 characteristic points by DInSAR

由表2、表3、圖4 分析可得：

1）該采空區(qū)工作面已停止開采，且附近無正在開采的工作面，因此監(jiān)測初期（2020-09-01—2020 -11-12）研究區(qū)域的上覆地表形變趨于穩(wěn)定。至礦震發(fā)生前夕（2020 年11 月下旬到2020 年12 月中下旬），監(jiān)測點中最大累計形變量達-27.1 mm。此次礦震發(fā)生時，周邊環(huán)境感受到異常，但未造成地表塌陷，因此監(jiān)測到形變值較小。

2）圖4 中，9 個監(jiān)測點在前6 個監(jiān)測時段（2020-09-01—2020-11-12）內(nèi)，地表緩慢下沉和緩慢抬升交替性出現(xiàn)，且抬升范圍在0.5～7.5 mm 之間，下沉數(shù)值小于-7 mm；在第7 個監(jiān)測時段（2020-11-12—2020-11-24）時，9 個特征點均監(jiān)測到形變的急劇變化，無抬升，最小下沉數(shù)據(jù)為-7.8 mm，最大下沉數(shù)據(jù)達-11.2 mm，9 個特征點在該時段的下沉數(shù)值均高于前6 個時段的下沉數(shù)值；從第7 個監(jiān)測時間段到第10 個監(jiān)測時間段（2020-11-12—2020-12-30），采空區(qū)持續(xù)下沉。綜合來看，2020-11-12—2020-12-30，采空區(qū)經(jīng)過形變量急劇變化后依然持續(xù)下沉，經(jīng)過長時間的下沉累計，采空區(qū)內(nèi)部承受不住地表重力，于2020 年12 月23 日發(fā)生礦震。

圖4 9 個DInSAR 特征點形變數(shù)據(jù)圖Fig.4 Settlement data diagrams of 9 characteristic points by DInSAR

3.2 SBAS InSAR 形變監(jiān)測結(jié)果

9 個特征點分布與SBAS InSAR 累計形變影像疊加圖如圖5。9 個SBAS InSAR 特征點累計形變量見表4，9 個SBAS InSAR 特征點形變速率見表5，9個SBAS InSAR 特征點形變速率圖如圖6。

圖5 9 個特征點分布與SBAS InSAR 累計形變影像疊加圖Fig.5 Overlaying charts of 9 characteristic points and the cumulative subsidence images by SBAS InSAR at each period

圖6 9 個SBAS InSAR 特征點形變速率圖Fig.6 Settlement rate diagrams of 9 characteristic points by SBAS InSAR

表4 9 個SBAS InSAR 特征點累計形變量Table 4 Cumulative settlement data of 9 characteristic points by SBAS InSAR

表5 9 個SBAS InSAR 特征點形變速率Table 5 Settlement rate of 9 characteristic points by SBAS InSAR

由圖5、表4、表5、圖6 分析可得，

1）根據(jù)WGS84 坐標系的參數(shù)進行計算，SBAS InSAR 監(jiān)測到發(fā)生地表形變的中心位置（35.532°N，116.918°E）與官方發(fā)布的礦震中心“D”的距離在460 m 左右。同樣，“D”處于SBA InSAR 監(jiān)測的地表形變范圍之內(nèi)，在2020 年9 月至2020 年12 期間，“D”附近地表經(jīng)歷了“抬升-下沉-抬升”。因此，SBAS InSAR 監(jiān)測到的地表形變中心在官方發(fā)布的礦震中心位置的誤差范圍內(nèi)。

2）SBAS InSAR 監(jiān)測結(jié)果顯示采空區(qū)工作面早已停止開采，上覆地表仍在持續(xù)下沉，地表形變活動一直在發(fā)生。礦震發(fā)生前夕（2020 年11 月下旬到2020 年12 月中下旬），SBAS InSAR 監(jiān)測到礦區(qū)最大累計形變達到-41.3 mm。因此，為防止由于持續(xù)形變影響地表結(jié)構(gòu)而產(chǎn)生礦震等強地表活動，防止危害附近人民安全和破壞生態(tài)環(huán)境，對采空區(qū)的監(jiān)測是必要的。

3）圖6 中，9 個監(jiān)測點在前7 個監(jiān)測時段（2020-09-01—2020-11-24）內(nèi)，形變速率在0～-0.4 mm/d之間平緩浮動；到第8 個監(jiān)測時段（2020-11-24—2020-12-06），監(jiān)測點的地表形變速率增大，其中3、5、6 3 個監(jiān)測點的形變速率有較明顯的增大，形變速率在數(shù)值上分別增大了0.14、0.13、0.21 mm/d；從第8 個監(jiān)測時間段到第10 個監(jiān)測時間段（2020-11-24—2020-12-30），監(jiān)測點的形變速率持續(xù)增大，形變量也持續(xù)增加?？梢姡摰V區(qū)采空區(qū)經(jīng)過長時間的形變累計，內(nèi)部承受不住地表重力而發(fā)生礦震。

為更進一步分析該采空區(qū)形變變化，對該區(qū)域進行剖面分析，考察其在時間序列上的形變變化。提取時序累計形變量繪制的采空區(qū)地表形變剖面圖如圖7。

圖7 采空區(qū)地表形變剖面圖Fig.7 Settlement profile of surface

圖7 清楚地反映出該采空區(qū)在時間域的形變量變化情況。由圖7 可見，隨著時間的推移，該采空區(qū)形變量在逐漸增加，2020 年11 月24 日，最大累計形變量達到-30.4 mm，自此之后，形變量開始增大。由此可以推斷，該采空區(qū)由于長期持續(xù)下沉，導(dǎo)致2020 年11 月24 日后，形變速率加快，地表結(jié)構(gòu)加速破壞，導(dǎo)致礦震的發(fā)生。

4 結(jié) 語

1）受礦震影響前，DInSAR、SBAS InSAR 技術(shù)均監(jiān)測到礦區(qū)上覆地表緩慢變化，監(jiān)測到的地表形變分別在-0.2 ～-6.9、-1.0 ～-5.2 mm 之間；礦震發(fā)生前夕（2020 年11 月下旬到2020 年12 月中下旬），2種InSAR 技術(shù)都監(jiān)測到采空區(qū)上覆地表形變發(fā)生了較明顯的變化，形變持續(xù)加快，形變量持續(xù)增大，地表持續(xù)下沉，DInSAR 監(jiān)測到最大累計形變量達到-27.0 mm，SBAS InSAR 監(jiān)測到礦區(qū)最大累計形變達到-41.3 mm；礦震發(fā)生期間，礦區(qū)地表仍持續(xù)下沉。

2）DInSAR 技術(shù)采用相鄰成像時刻的2 幅SAR影像兩兩差分干涉處理獲取相鄰時刻之間的地表形變信息，處理耗時，低相干點的監(jiān)測精度不高；SBAS InSAR 技術(shù)采用的是時間序列的差分干涉處理得到各成像時刻的累計地表形變信息，處理過程復(fù)雜，無法得到低相干點的地表形變信息；二者各有優(yōu)缺點。

3）2 種InSAR 技術(shù)的監(jiān)測結(jié)果均體現(xiàn)出煤礦采空區(qū)經(jīng)過持續(xù)形變，會引發(fā)地表形變加劇，地表結(jié)構(gòu)遭到破壞，內(nèi)部承受不住地表重力，導(dǎo)致礦震發(fā)生。研究結(jié)果對DInSAR、SBAS InSAR 技術(shù)應(yīng)用于礦震預(yù)測提供一定的參考，但尚需結(jié)合其它更多的礦震實例做進一步的深入研究。