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

（山東科技大學(xué) 測(cè)繪與空間信息學(xué)院，山東 青島 266590）

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

合成孔徑雷達(dá)干涉測(cè)量InSAR（Interfero -metry Synthetic Aperture Radar）技術(shù)因其擁有全天候、全天時(shí)、高精度、高空間分辨率及覆蓋面廣等優(yōu)點(diǎn)逐步用于礦區(qū)的地表形變監(jiān)測(cè)[6]。隨著該技術(shù)的不斷完善，繁衍出DInSAR（Differential InSAR）、SBAS In－SAR（Small BAseline Subset InSAR）等多種地表形變監(jiān)測(cè)手段。許多研究人員將InSAR 技術(shù)應(yīng)用于采空區(qū)地表形變監(jiān)測(cè)，姚佳明等[7]選用升、降軌L 波段PALSAR-2 影像數(shù)據(jù)，利用InSAR 技術(shù)對(duì)煤礦采空區(qū)開展了短期動(dòng)態(tài)地表沉降監(jiān)測(cè)，結(jié)合研究區(qū)開采信息對(duì)煤礦采空范圍及開采時(shí)間進(jìn)行反演，驗(yàn)證了InSAR技術(shù)對(duì)煤礦采區(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ā)展過程，證實(shí)SBAS In－SAR 技術(shù)能夠獲得采空區(qū)的毫米級(jí)沉降結(jié)果。

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

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

1.1 研究區(qū)概況

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

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

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

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

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

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

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

2.1 DInSAR 原理

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

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

2.2 SBAS InSAR 原理

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

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

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

則相位為：

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

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

通過奇異值分解法和最小二乘法可求出平均形變速率相位值，得到地表累計(jì)線性形變量[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ù)處理、影像配準(zhǔn)及重采樣、干涉圖生成、基線估計(jì)、地形相位去除、差分干涉圖濾波、相位解纏、地理編碼等，通過相位-形變轉(zhuǎn)換，得到研究區(qū)的地表形變信息。

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

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

3 實(shí)驗(yàn)結(jié)果分析

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

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

3.1.1 特征點(diǎn)DInSAR影像疊加圖分析

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

由圖2、圖3 可以看出：

圖2 9 個(gè)特征點(diǎn)分布與DInSAR 各時(shí)段地表形變影像疊加圖Fig.2 Overlaying charts of 9 characteristic points and DInSAR subsidence images at each period

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

1）DInSAR 監(jiān)測(cè)到的地表形變中心（35.533°N，116.918°E）與官方發(fā)布的礦震中心“D”（35.54°N，116.92°E）相距約417 m。官方給出的礦震中心位置僅保留了小數(shù)點(diǎn)后2 位，是1 個(gè)大致位置，且“D”處于DInSAR 監(jiān)測(cè)的地表形變范圍之內(nèi)，2020-09—2020-12 期間，“D”處地表經(jīng)歷了“抬升-下沉-抬升”。因此，DInSAR 監(jiān)測(cè)到的地表形變中心在官方發(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 個(gè)時(shí)段內(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 個(gè)時(shí)段緩慢下沉。在第7 監(jiān)測(cè)時(shí)段（2020－11－12—2020-11-24）內(nèi)，地表出現(xiàn)-8～-10 mm 之間的不規(guī)則的大面積下沉，監(jiān)測(cè)區(qū)域內(nèi)最大形變量達(dá)到-11.5 mm。在第8 監(jiān)測(cè)時(shí)段（2020-11-24—2020-12-06），礦區(qū)中心下沉，周邊部分區(qū)域發(fā)生抬升，第9 監(jiān)測(cè)時(shí)段（2020-12-06— 2020-12-18）礦區(qū)中心仍持續(xù)下沉，周邊區(qū)域大面積抬升，地表高度差逐漸增大，地表結(jié)構(gòu)變形，導(dǎo)致在第10 個(gè)監(jiān)測(cè)時(shí)段（2020-12-18—2020-12-30）內(nèi)發(fā)生礦震。

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

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

表2 9 個(gè)DInSAR 特征點(diǎn)各時(shí)段形變量Table 2 Settlement data at each period of 9 characteristic points by DInSAR

表3 9 個(gè)DInSAR 特征點(diǎn)累計(jì)形變量Table 3 Cumulative settlement data of 9 characteristic points by DInSAR

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

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

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

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

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

9 個(gè)特征點(diǎn)分布與SBAS InSAR 累計(jì)形變影像疊加圖如圖5。9 個(gè)SBAS InSAR 特征點(diǎn)累計(jì)形變量見表4，9 個(gè)SBAS InSAR 特征點(diǎn)形變速率見表5，9個(gè)SBAS InSAR 特征點(diǎn)形變速率圖如圖6。

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

圖6 9 個(gè)SBAS InSAR 特征點(diǎn)形變速率圖Fig.6 Settlement rate diagrams of 9 characteristic points by SBAS InSAR

表4 9 個(gè)SBAS InSAR 特征點(diǎn)累計(jì)形變量Table 4 Cumulative settlement data of 9 characteristic points by SBAS InSAR

表5 9 個(gè)SBAS InSAR 特征點(diǎn)形變速率Table 5 Settlement rate of 9 characteristic points by SBAS InSAR

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

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

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

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

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

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

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

4 結(jié) 語(yǔ)

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

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

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