基于SBAS-InSAR 的采煤區(qū)地表形變監(jiān)測(cè)與分析

宋旭斌，路 鑫

（潞安化工集團(tuán)有限公司 古城煤礦，山西 長治 046000）

0 引 言

煤炭作為我國的主體能源，長期擔(dān)負(fù)著國家能源安全和經(jīng)濟(jì)高速發(fā)展的重任。但大規(guī)模、高強(qiáng)度的地下開采往往會(huì)破壞地層的原始應(yīng)力平衡，而且隨著采空區(qū)范圍擴(kuò)大，直接頂在自身重力和上覆巖層荷載作用下破斷，最終引起地表的移動(dòng)變形[1]。因開采導(dǎo)致的地表沉陷易造成建筑物墻體開裂、交通道路損壞、農(nóng)田積水等地質(zhì)災(zāi)害問題，嚴(yán)重制約了礦區(qū)城市的可持續(xù)發(fā)展。因此，對(duì)采空區(qū)地表進(jìn)行實(shí)時(shí)精密形變監(jiān)測(cè)，獲取動(dòng)態(tài)沉降信息勢(shì)在必行。

常規(guī)開采監(jiān)測(cè)手段主要有精密水準(zhǔn)測(cè)量和全球?qū)Ш较到y(tǒng)（Global Navigation Satellite System,GNSS），但在實(shí)際工程中往往面臨觀測(cè)點(diǎn)易丟失、觀測(cè)精度易受干擾、成本投入巨大、時(shí)空分辨率低等問題，無法準(zhǔn)確獲取礦區(qū)整體的形變特征。合成孔徑雷達(dá)干涉測(cè)量(Interferometry Synthetic Aperture Radar,InSAR)是一種新型的對(duì)地觀測(cè)技術(shù)，它僅對(duì)同一區(qū)域的多景SAR 影像進(jìn)行復(fù)共軛相乘即可提取大范圍、毫米級(jí)的地面形變信息[2]，具有全天時(shí)、高精度、自動(dòng)化等優(yōu)勢(shì)。而小基線集（Small Baseline Subsets InSAR,SBAS-InSAR）技術(shù)在差分InSAR（Differential InSAR,D-InSAR）基礎(chǔ)上提出，有效克服了后者易受大氣效應(yīng)、時(shí)空失相干影響的缺陷[3]，它通過對(duì)SAR 影像上的高相干點(diǎn)進(jìn)行相位模型建立和相位變化分析，即使在永久散射體PS（Persistent Scatterer）密度較低的礦區(qū)也能獲取更為精確的長時(shí)序形變結(jié)果。

許多學(xué)者都曾將SBAS 技術(shù)運(yùn)用于礦區(qū)的形變監(jiān)測(cè)，并取得了一系列豐碩的成果。李達(dá)等[4]利用SBAS 技術(shù)對(duì)13 景TerraSAR-X 進(jìn)行處理，獲取了榆林某礦的地表沉降信息，發(fā)現(xiàn)地表下沉值與時(shí)間呈線性關(guān)系；欒元重等[5]借助19 景Sentinel-1A 影像和SBAS 技術(shù)對(duì)沉陷漏斗區(qū)域內(nèi)的工作面進(jìn)行時(shí)序監(jiān)測(cè)，并利用水準(zhǔn)數(shù)據(jù)驗(yàn)證了SBAS 技術(shù)的可靠性；柴華彬等[6]提出利用反距離加權(quán)插值方法（Inverse Distance Weighted,IDW）將實(shí)測(cè)水準(zhǔn)數(shù)據(jù)與SBAS 解譯結(jié)果相融合，以解決部分區(qū)域因失相干導(dǎo)致InSAR 形變值缺失的問題；張香凝等[7]以蒲河礦區(qū)為例，結(jié)合SBAS 技術(shù)和數(shù)值模擬共同探究了煤礦區(qū)地面沉降的時(shí)空變形規(guī)律和機(jī)制。

山西長治地區(qū)煤炭資源豐富，但因地下開采引起的地質(zhì)災(zāi)害頻發(fā)，亟需開展相關(guān)的地表沉陷監(jiān)測(cè)工作。本文在前人研究的基礎(chǔ)上，借助SBAS-InSAR 技術(shù)對(duì)29 景Sentinel-1A 影像進(jìn)行解譯處理，以期獲取古城煤礦S1306 工作面在2021年10 月12 日至2022 年10 月7 日期間的時(shí)序沉降結(jié)果和開采影響范圍。研究成果對(duì)礦區(qū)開采規(guī)劃和地表塌陷的預(yù)警防治均具有重要的現(xiàn)實(shí)意義。

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

1.1 研究區(qū)概況

本文研究區(qū)位于長治市古城煤礦東南方向一采區(qū)，整體地勢(shì)呈西高東低，地貌類型以河谷平原為主。研究區(qū)分布著大面積農(nóng)田、村莊房屋以及電力線等重要的基礎(chǔ)設(shè)施，且覆蓋有較厚的第四系黃土沉積；地下山西組3 號(hào)煤層目前有多個(gè)工作面正在回采，其中S1306 工作面走向長約1 982 m，傾向?qū)捈s363 m，平均采厚6.4 m，煤層傾角約為0～6°，底板標(biāo)高-352—-397 m。開采方法以傾斜長壁式綜采法為主，全部垮落法管理頂板。目前，該采區(qū)周邊已出現(xiàn)大量明顯的地裂縫及路面損壞情況，亟需對(duì)該區(qū)域進(jìn)行地表沉陷監(jiān)測(cè)工作。研究區(qū)的地理位置如圖1 所示。

圖1 研究區(qū)地理位置Fig.1 Geographical location of the study area

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

本文從歐空局（European Space Agency,ESA）哥白尼開放數(shù)據(jù)訪問中心（https://scihub.copernicus.eu/dhus/#/home）收集了29 景Sentinel-1A 升軌單視復(fù)數(shù)影像，時(shí)間跨度為2021 年10 月12 日至2022 年10 月7 日，工作模式均為干涉寬幅（Interferometric Wide-swath，IW）模式，該模式基于遞進(jìn)式地形掃描(Terrain Observation with Progressive Scans SAR，TOPSAR)技術(shù)，能夠?qū)? 個(gè)子條帶合并成1 景高質(zhì)量的SAR 雷達(dá)數(shù)據(jù)[11]，其地面分辨率可達(dá)5 m×20 m。實(shí)驗(yàn)影像的其他參數(shù)見表1。另外，在數(shù)據(jù)處理過程中，引入了POD精密軌道文件用來去除衛(wèi)星軌道的系統(tǒng)性誤差，并利用30 m 分辨率的STRM1 DEM 數(shù)據(jù)消除了干涉圖中地形相位的影響。

表1 Sentinel-1A 基本參數(shù)Table 1 Basic parameters of Sentinel-1A

2 技術(shù)原理與數(shù)據(jù)處理流程

2.1 SBAS-InSAR 原理

小基線集技術(shù)最早由國外學(xué)者Berardino 等人在2002 年提出，其基本原理如下[12-15]：假設(shè)選?。╰0,t1,t2,…,tN）時(shí)段內(nèi)覆蓋同一區(qū)域的N+1 幅SAR影像，基于設(shè)定的時(shí)空基線閾值，組合成M 幅多視差分干涉對(duì)，其中M 的取值范圍滿足：

若在tA和tB（tB>tA）時(shí)刻對(duì)2 幅SAR 影像進(jìn)行干涉處理，并生成了第j 幅已解纏的差分干涉圖，在去除地形相位、大氣延遲相位和噪聲誤差之后，任意像元位置（x,r）（x 和r 分別為方位向和距離向坐標(biāo)）處的差分干涉相位可表示為：

式中：λ 為雷達(dá)波長；φ(tB,x,r)和φ(tA,x,r)分別表示2 幅影像上的相位值；d(tB,x,r)和d(tA,x,r)分別為tA和tB時(shí)刻相對(duì)于初始時(shí)刻t0在雷達(dá)視線向(Line of Sight,LOS)上的形變，一般假設(shè)t0時(shí)刻的形變相位值為0；vj為tA至tB時(shí)段內(nèi)的平均形變速率。

將用于干涉處理的M 對(duì)影像數(shù)據(jù)按照時(shí)序進(jìn)行排列，則所有干涉圖相位可表示為如下矩陣形式：

式中：B 為1 個(gè)M×N 的系數(shù)矩陣；V 為速度矢量。

由于矩陣B 的秩小于N，根據(jù)最小二乘法無法得到唯一的結(jié)果，因此需要采用奇異值分解法（Singular Value Decomposition,SVD）聯(lián)合多個(gè)基線集，來解決法方程的秩虧問題，并得到速度V 的廣義逆解。之后將各個(gè)時(shí)間段內(nèi)的形變速率在時(shí)間域上進(jìn)行積分，即可獲取整個(gè)監(jiān)測(cè)時(shí)段內(nèi)地表在LOS 向上的累積時(shí)序形變結(jié)果。

2.2 SBAS-InSAR 數(shù)據(jù)處理流程

SBAS-InSAR 數(shù)據(jù)處理關(guān)鍵步驟和參數(shù)設(shè)置如下：①選取2021 年11 月29 日的影像作為主影像，其余從影像與之配準(zhǔn)，并設(shè)置空間基線百分比閾值為45%，時(shí)間基線閾值為60 d，多視比為4∶1（距離向∶方位向），最后共生成104 個(gè)干涉對(duì)；②利用最小費(fèi)用流(Minimum Cost Flow,MCF)方法進(jìn)行相位解纏，Goldstein 濾波技術(shù)去除干涉圖噪聲，并逐個(gè)檢查調(diào)整不理想的像對(duì)；③參考干涉處理步驟生成的相干圖和解纏相位圖，選取20～30 個(gè)穩(wěn)定且分布均勻的地面控制點(diǎn)（Ground Control Point,GCP）完成軌道精煉，去除殘余軌道誤差相位；④選擇線性模型（Linear Model）進(jìn)行二次解纏，并估算平均位移速率和殘余地形相位；⑤根據(jù)大氣延遲相位、形變相位和噪聲相位表現(xiàn)出的不同的時(shí)空特性，利用時(shí)間域的高通濾波（365 d）和空間域的低通濾波（1 200 m）分離并剔除大氣延遲誤差，得到覆蓋整個(gè)觀測(cè)時(shí)間的地表形變時(shí)間序列；⑥將最終解譯結(jié)果從斜距坐標(biāo)系轉(zhuǎn)換至WGS-84 坐標(biāo)系，即完成地理編碼。本文的數(shù)據(jù)處理流程如圖2所示。

圖2 SBAS 數(shù)據(jù)處理流程Fig.2 Data dealing process of SBAS

3 解譯結(jié)果與分析

3.1 空間演化分析

利用SBAS-InSAR 技術(shù)對(duì)覆蓋S1306 工作面的29 景影像進(jìn)行解譯處理，獲取了監(jiān)測(cè)時(shí)段內(nèi)研究區(qū)LOS 向的時(shí)序累積形變結(jié)果。之后忽略水平位移的影響，利用公式dv=dLOS/cosθ 將解譯結(jié)果從視線向投影至垂直向，其中θ 為衛(wèi)星入射角。另外，解譯結(jié)果中的正值表示地面抬升，負(fù)值表示地面下沉。

S1306 工作面自2022 年1 月開始回采，目前僅開采到距離始采線約850 m 的位置。但圖3 顯示，在監(jiān)測(cè)時(shí)段內(nèi)，該工作面周邊地表已經(jīng)發(fā)生了明顯沉陷。

圖3 S1306 工作面時(shí)序形變Fig.3 Time series deformation of No.S1306 Face

在形變監(jiān)測(cè)的初期階段，即2021 年10 月12日至2021 年12 月23 日（圖3a～圖3c），研究區(qū)地表總體較為穩(wěn)定，零星區(qū)域發(fā)生小幅沉降，量值約-13 mm，產(chǎn)生該現(xiàn)象的原因可歸結(jié)為工作面所處的農(nóng)田區(qū)域相干性較差，導(dǎo)致SBAS 在反演形變的過程中發(fā)生解纏誤差，影響了結(jié)果的準(zhǔn)確性；從2022 年1 月16 日起（圖3d），工作面的始采線位置開始出現(xiàn)明顯下沉，最大沉降量約-32.6 mm，符合工作面開始回采的時(shí)間；隨著工作面持續(xù)推進(jìn)，地表沉降范圍逐漸朝西南方向（即工作面開采方向）擴(kuò)張；到了2022 年5 月16 日（圖3h），始采線附近地表已經(jīng)形成了一個(gè)較為明顯的沉降盆地，盆地中心可探測(cè)到的最大下沉值約-178.6 mm；隨著采動(dòng)影響逐漸傳遞至工作面中部位置（圖3i～圖3k），始采線區(qū)域的下沉速度逐漸趨于平緩，但仍處在下沉狀態(tài)；直至2022 年10 月7 日（圖3m），此時(shí)S1306 工作面的開采影響邊界稍超過實(shí)際開采進(jìn)度位置，但總體與實(shí)際情況相符，說明SBAS 技術(shù)在礦區(qū)形變監(jiān)測(cè)中具有一定的可靠性和準(zhǔn)確性。

3.2 時(shí)間演化分析

為了進(jìn)一步探究研究區(qū)地表在采動(dòng)影響下的時(shí)序形變過程，在S1306 工作面的走向和傾向上共選取8 個(gè)特征點(diǎn)（編號(hào)點(diǎn)A～點(diǎn)H），并以2021 年10 月12 日影像為參考基準(zhǔn)，生成各個(gè)點(diǎn)的時(shí)序累積沉降曲線。具體點(diǎn)位分布如圖4 所示。

圖4 特征點(diǎn)位及其時(shí)序累積沉降值Fig.4 The location of feature points and their time-series cumulative subsidence values

由圖4 可知，在2021 年10 月12 日至2022 年1 月16 日期間，走向線上4 個(gè)特征點(diǎn)沒有發(fā)生明顯形變，由于點(diǎn)A 和點(diǎn)B 分別位于始采線和沉陷盆地中心，易受開采擾動(dòng)，這一時(shí)段內(nèi)兩點(diǎn)均出現(xiàn)小幅度下沉，最大沉降值約-10 mm；而點(diǎn)C 和點(diǎn)D 分布在工作面中部和停采線的位置，點(diǎn)位表現(xiàn)出反復(fù)微弱抬升的跡象，可將這種情況其視為SAR圖像噪聲帶來的解算誤差。

2022 年1 月16 日之后，點(diǎn)A 和點(diǎn)B 開始持續(xù)加速下沉，且分別于2022 年3 月17 日和2022 年5 月16 日，其形變速率才逐漸放緩，但形變量仍在持續(xù)增加，這是由于采動(dòng)作用對(duì)始采線區(qū)域地表的影響即將達(dá)到“臨界點(diǎn)”，該區(qū)域的形變經(jīng)歷了“開始下沉—?jiǎng)×蚁鲁痢徛Ｖ瓜鲁痢钡倪^程，一定程度上符合常規(guī)開采沉陷規(guī)律。點(diǎn)C 和點(diǎn)D在2022 年8 月8 日之前，下沉曲線吻合度較高，但之后點(diǎn)C 的沉降量陡增至-87.6 mm，說明工作面中部區(qū)域此時(shí)開始受采動(dòng)的強(qiáng)烈影響，與實(shí)際工程進(jìn)度一致。同時(shí)，發(fā)生迅速沉降的原因還可能是由于點(diǎn)C 處在土質(zhì)松軟的農(nóng)田，在夏季強(qiáng)降水的滲透下，土層的抗剪強(qiáng)度會(huì)大大降低，從而引起地表的濕化崩解和垮落[12]。

傾向線上的4 個(gè)特征點(diǎn)均勻分布在工作面的兩側(cè)，在2022 年1 月16 日之前，由于點(diǎn)E 和點(diǎn)H遠(yuǎn)離開采中心，點(diǎn)F 和點(diǎn)G 位于沉陷盆地邊緣，因此后者呈現(xiàn)小幅持續(xù)下沉的態(tài)勢(shì)，前者沉降與起伏交替，受采動(dòng)影響較??；之后點(diǎn)F 和點(diǎn)G 開始加速下沉，但點(diǎn)F 的沉降幅度與速率要大于點(diǎn)G。同樣，點(diǎn)E 和點(diǎn)H 于2022 年4 月10 日也開始經(jīng)歷一個(gè)緩慢下沉的過程，但總體上點(diǎn)E 的沉降量級(jí)要高于點(diǎn)H。

因此結(jié)合圖3 的時(shí)序形變結(jié)果，可以發(fā)現(xiàn)沿開采方向右側(cè)的地表形變值要大于左側(cè)區(qū)域，即始采線附近的沉陷盆地實(shí)際上有朝西北方向發(fā)育的趨勢(shì)，后期需要對(duì)東李高村附近地表進(jìn)行重點(diǎn)監(jiān)測(cè)。

4 結(jié) 論

（1）在采動(dòng)作用下，S1306 工作面周邊地表自2022 年1 月16 日開始出現(xiàn)下沉，之后在始采線附近迅速產(chǎn)生了一個(gè)沉陷盆地，且沉陷盆地有朝西北方向發(fā)育的趨勢(shì)。

（2）目前采動(dòng)影響僅波及至工作面的中部區(qū)域，與實(shí)際開采到達(dá)位置基本吻合，說明SBAS-InSAR 技術(shù)在礦區(qū)形變監(jiān)測(cè)中具有良好的適用性。

（3）除了利用時(shí)序InSAR 技術(shù)獲取不同時(shí)段的地表沉陷情況，仍可以考慮利用UAV 激光雷達(dá)掃描技術(shù)或建立地面移動(dòng)觀測(cè)站來對(duì)局部重點(diǎn)區(qū)域進(jìn)行詳細(xì)監(jiān)測(cè)，真正發(fā)揮“空-天-地”協(xié)同監(jiān)測(cè)技術(shù)在開采沉陷領(lǐng)域的重要作用。