基于SBAS-InSAR 技術(shù)的楊伙盤礦區(qū)地表形變監(jiān)測

張 帆，常 樂，荀張媛

（1.西安煤航遙感信息有限公司，陜西 西安 710100；2.沈陽城市建設(shè)學(xué)院，遼寧 沈陽 110167）

0 引言

我國是礦產(chǎn)資源大國，煤炭等資源一直在能源開發(fā)利用中占據(jù)重要地位。煤炭資源的長期開采不僅引起地面塌陷等一系列問題，還會造成建筑物受損、農(nóng)田被破壞[1]，以及誘發(fā)崩塌或滑坡等次生災(zāi)害[2]，對人民的生產(chǎn)生活造成嚴(yán)重的經(jīng)濟損失，甚至?xí)＜吧踩玔3-4]。因此，重視礦區(qū)地表形變，實現(xiàn)礦區(qū)生態(tài)環(huán)境可持續(xù)發(fā)展是科學(xué)開采煤炭等礦產(chǎn)資源的重要任務(wù)。

礦區(qū)地表形變是一個基于面區(qū)域的大范圍的地表運動，但是在實際生產(chǎn)過程中，更多的還是依賴于傳統(tǒng)的水準(zhǔn)測量、GPS 測量等點監(jiān)測技術(shù)，在生產(chǎn)工作面上每隔一段距離布設(shè)一個水準(zhǔn)點或GPS 點。一方面，布設(shè)較多的點，人工消耗成本高、周期長；另一方面，水準(zhǔn)測量和GPS 測量都是基于點觀測，對于礦區(qū)大范圍大區(qū)域監(jiān)測，需要大量加密點，通過擬合獲取最終結(jié)果，無法保證精度[5-7]。

合成孔徑雷達(dá)干涉測量（InSAR）技術(shù)是一種高效對地觀測技術(shù)，具有全天時、全天候?qū)Φ爻上?，監(jiān)測結(jié)果精確度高，測量范圍廣的優(yōu)勢，大大降低人力、物力、財力的消耗[8-10]。特別是時間序列InSAR 技術(shù)能夠更快速獲取地表大范圍高精度時間序列形變信息，為礦區(qū)開采沉陷監(jiān)測提供了新的方向[11-13]。

越來越多的學(xué)者已經(jīng)將InSAR 技術(shù)應(yīng)用于煤礦地表監(jiān)測研究。尹宏杰等[14]利用時間序列InSAR 技術(shù)對湖南冷水江錫煤礦地面形變進行監(jiān)測，證實了煤礦地表時間序列形變漏斗的發(fā)展過程；趙偉穎等[15]利用SBAS-InSAR 技術(shù)獲取了徐州礦區(qū)地表的時間序列形變及累積形變量，并與水準(zhǔn)測量結(jié)果進行了外部精度符合驗證，且兩者的相關(guān)性較高；張香凝等[16]采用SBAS-InSAR 技術(shù)獲取了煤礦地表形變的時空分布特征，并利用數(shù)值模擬技術(shù)對觀測形變結(jié)果進行模擬分析，獲取了地面沉降在時間和空間上的變形規(guī)律和機制；姚佳明等[17]選用升軌、降軌L波段PALSAR-2 影像數(shù)據(jù)，利用InSAR 技術(shù)對煤礦采空區(qū)開展了短期動態(tài)地表沉降監(jiān)測，結(jié)合研究區(qū)開采信息對煤礦采空范圍及開采時間進行反演，驗證了InSAR 技術(shù)對煤礦采區(qū)反演的合理性與可靠性；何榮興等[18]介紹了采空區(qū)災(zāi)害類型及不同類型的相應(yīng)案例，分析了采空區(qū)災(zāi)害發(fā)生的一般規(guī)律和特征，提出了盡量選擇不產(chǎn)生采空區(qū)的采礦方法。

本文以楊伙盤煤礦為例，采用SBAS-InSAR 技術(shù)，利用中分辨率Sentinel-1A 數(shù)據(jù)對煤礦地表進行時間序列形變監(jiān)測分析，提取礦區(qū)開采沉陷形變場，并結(jié)合光學(xué)影像數(shù)據(jù)，研究其沉陷變化規(guī)律。

1 研究區(qū)概況

陜北煤炭生產(chǎn)基地是全國煤炭生產(chǎn)基地之一，尤其是神木、府谷等區(qū)域，煤炭資源豐富。但是長期高強度的資源開采，加之陜北地區(qū)大多為黃土地貌，最終導(dǎo)致礦區(qū)的土地大面積沉陷，并引發(fā)了滑坡、崩塌等多種次生災(zāi)害。因此，為保護礦區(qū)生態(tài)環(huán)境，減少人民的財產(chǎn)損失，需要對礦區(qū)的采空區(qū)形變進行長時間的有效監(jiān)測；結(jié)合多種技術(shù)，加大對采煤塌陷區(qū)域的監(jiān)測監(jiān)管力度，形成一套行之有效的陜北榆神府礦區(qū)采空區(qū)沉陷監(jiān)測體系。

楊伙盤煤礦位于陜西省神木市城北約30 km 的府店一級公路北側(cè)，行政區(qū)劃屬神木市店塔鎮(zhèn)，靠近府谷縣（圖1），工作區(qū)面積約26.92 km2。礦區(qū)地表廣泛覆蓋著現(xiàn)代風(fēng)積沙及第四系黃土層、新近系紅土層等，僅部分地區(qū)有基巖出露，基巖以砂巖、泥巖為主。根據(jù)某一鉆孔資料，3-1煤層主要位于侏羅系中統(tǒng)延安組第三段，煤層距地表約160 m。根據(jù)現(xiàn)有巷道布置情況，結(jié)合各煤層開采范圍及煤層賦存特點，一水平3-1煤層劃分4 個盤區(qū)。圖2 展示了3-1煤層開采工作面分布情況。

圖1 楊伙盤煤礦地理位置及高程Fig.1 Geographical location and elevation of Yanghuopan Coal Mine

圖2 工作面分布圖Fig.2 Distribution map of working face

2 InSAR 數(shù)據(jù)及技術(shù)方法

2.1 InSAR 技術(shù)方法

BERARDINO 等[19]提出的SBAS-InSAR 技術(shù)，是目前較為常用的一種時間序列InSAR 技術(shù)，技術(shù)流程如圖3 所示。首先，根據(jù)SAR 數(shù)據(jù)量，通過設(shè)置時間基線閾值和空間基線閾值組合干涉對，一般根據(jù)相干系數(shù)圖和干涉圖挑選質(zhì)量較高的干涉對；其次，通過多視和濾波等操作提高所選干涉圖的信噪比，基于具有穩(wěn)定后向散射特性的高相干點去除地形相位和平地相位，并進行解纏；之后，根據(jù)干涉相位及時空基線去除DEM 誤差，并分離線性形變相位，解纏殘余相位，通過時間域高通濾波和空間域低通濾波分離出殘余相位中的大氣相位，剩下即為非線性相位，補償線性形變相位即可得到完整的形變相位；最后，采用奇異值分解，連接多個子集，求得影像序列地表形變速率的最小范數(shù)、最小二乘解[20]。該方法能有效去除或減弱時空失相干等誤差的影響，并在一定程度消除大氣相位和限制地形誤差的影響。

圖3 SBAS 技術(shù)流程圖Fig.3 Flow chart of SBAS technology

獲取研究 區(qū)t0······tn時間段內(nèi)N+1幅SAR 影像，將影像配準(zhǔn)到同一個坐標(biāo)系，設(shè)定時間基線閾值及空間基線閾值，得到M個差分干涉對。假設(shè)其中第i幅干涉圖由tA、tB（tA＜tB）兩時期影像生成，此干涉圖上雷達(dá)坐標(biāo)系下（a,r）處的干涉相位可表示為式（1）。

式中：δφi（disp）為形變相位；δφi（topo）為殘余地形相位；δφi（aps）為大氣延遲相位；δφi（noise）為噪聲相位。忽略殘余地形相位、大氣延遲相位及噪聲相位，將形變相位以LOS向形變表示，式（1）可簡化為式（2）。

式中：λ為雷達(dá)波長；dtA（a,r）、dtB（a,r）為（a,r）處在tA時期、tB時期相對于初始t0時期的形變。針對（a,r）這一像元處，假設(shè)IE為主影像，IS為從影像，對應(yīng)全部的M個干涉對，見式（3）。

式（3）滿足IEi＞ISi，其中，i=1,···,M。則對于任意干涉圖，對應(yīng)的有式（4）。

式（4）為n個未知數(shù)M個方程所組成的方程組，以矩陣形式表示 δφ=Aφ，其中，A為M*n矩陣。對于小基線子集內(nèi)部，M≥n，方程可由最小二乘法求解。各子集之間M接近或者小于n時方程會出現(xiàn)多解，采用奇異值分解以及最小范數(shù)進行解算獲得唯一解[19]。

2.2 SAR 數(shù)據(jù)

本次研究使用的數(shù)據(jù)為歐洲航天局Sentinel-1A衛(wèi)星在TOPS 成像模式下的寬幅干涉SLC 數(shù)據(jù)，為目前常用的免費中分辨率SAR 數(shù)據(jù)，數(shù)據(jù)產(chǎn)品詳細(xì)參數(shù)見表1。根據(jù)監(jiān)測周期，共收集了2020 年8 月—2022 年11 月共58 景升軌數(shù)據(jù)，軌道號為113/126。

表1 Sentinel-1A IW 產(chǎn)品基本參數(shù)Table 1 Basic parameters of Sentinel-1A IW product

本次采用的DEM 數(shù)據(jù)為ALOS Global Digital Surface Model“ALOS World 3D-30 m”（AW3D30），是由日本宇宙航空研究開發(fā)機構(gòu)（JAXA）2015 年5 月免費提供的高精度全球數(shù)字地表模型數(shù)據(jù)，水平分辨率為30 m（1"），高程精度5 m，是目前世界上最精確的3D 地圖，覆蓋全球所有的土地尺度，工作區(qū)的DEM 數(shù)據(jù)如圖1 所示。精密軌道數(shù)據(jù)采用成像21 d之后發(fā)布的POD 精密軌道數(shù)據(jù)。

3 結(jié)果與分析

3.1 InSAR 形變分析

工作區(qū)年平均形變速率圖如圖4 所示。形變部分主要分布在①、②、③三個區(qū)域，沿著工作面均表現(xiàn)為長條狀分布，主要覆蓋的工作面有30112 工作面、30114 工作面、30116 工作面、30115 工作面、30117工作面、30301 工作面和30302 工作面。三個形變中心可探測的最大形變速率分別約為-52 mm/a、-39 mm/a 和-39 mm/a，形變最大的區(qū)域位于區(qū)域①，分別對三個形變區(qū)進行時間序列和剖線形變分析。

圖4 工作區(qū)年平均形變速率圖Fig.4 Annual average deformation rate of the working area

圖5 為區(qū)域①的時間序列形變圖。由圖5 可知，從2020 年8 月開始，30112 工作面處于開采狀態(tài)，圖5（a）為基準(zhǔn)面，在圖5（b）中地表開始出現(xiàn)形變。隨著開采面的推進，地表形變范圍逐漸擴大，直至監(jiān)測期結(jié)束，最終形變主要覆蓋的工作面為30112 工作面、30114 工作面和30116 工作面。其中，30114 工作面地表形變量級及范圍都較大，最大累積形變量約為-130 mm，為整個工作區(qū)的可探測到的最大形變量，在30114 工作面上（圖4）畫一條剖線AA′以及AA′剖線上的形變最值特征點a 進行時間序列分析。

圖5 區(qū)域①時間序列形變圖Fig.5 Time series deformation of region ①

圖6 為剖線AA′的時間序列剖線圖，其中主要選擇了6 個時間段顯示，共形成了4 個形變漏斗。由圖6 可知，靠近A 點的第一個形變漏斗的變化相對較為平穩(wěn)；第二個形變漏斗在2021 年1 月后發(fā)生突變；第三個漏斗在2021 年5 月后發(fā)生突變；第四個形變漏斗在2021 年8 月發(fā)生突變，這三個突變情況與工作面的開采時間都較為吻合。其中，第四個形變漏斗的形變量遠(yuǎn)遠(yuǎn)超過其他區(qū)域的形變量。圖7 為剖線AA′上的形變最值點a 的時間序列形變折線圖。由圖7 可知，2021 年8 月，該點發(fā)生了突變，而后速率逐漸減小。2022 年5 月前后，再次出現(xiàn)較小幅度的突變，可能是受到30116 工作面開采的影響。

圖6 AA'時間序列剖線圖Fig.6 Time series deformation of AA' profile

圖7 a 點時間序列圖Fig.7 Time series deformation of point a

圖8 為區(qū)域②的時間序列形變圖。由圖8 可知，從2020 年8 月開始，30115 工作面處于開采狀態(tài)。由于初期的形變較小，絕對量值在20 mm 之內(nèi)，在圖8（d）中才顯示出來。最終形變主要覆蓋30115 工作面和30117 工作面。在30117 工作面上（圖4）畫一條剖線BB′以及BB′剖線上的形變最值特征點b 進行時間序列分析。

圖8 區(qū)域②時間序列形變圖Fig.8 Time series deformation in region ②

圖9 為剖線BB′的時間序列剖線圖，其中，選擇了6 個時間段顯示，共形成多個形變漏斗。由圖9 可知，靠近B 點的形變突變基本都發(fā)生在2021 年1 月之后，后面時間段的形變變化相對都較為平穩(wěn)。圖10 為剖線BB′上的形變最值點b 的時間序列形變折線圖。由圖10 可知，在2021 年3 月前后，b 點發(fā)生了突變，與工作面開采的時間基本吻合，4 月之后速率逐漸減小。

圖9 BB'時間序列剖線圖Fig.9 Time series deformation of BB' profile

圖10 b 點時間序列圖Fig.10 Time series deformation of point b

圖11 為區(qū)域③的時間序列形變圖。由圖11 可知，2021 年此區(qū)域的30301 工作面才開始開采；2021年9 月，工作面的形變才開始顯示。最終，形變主要覆蓋的工作面為30301 工作面和30302 工作面。在30302 工作面上（圖4）畫一條剖線CC′以及CC′剖線上的形變最值特征點c 進行時間序列分析。

圖11 區(qū)域③時間序列形變圖Fig.11 Time series deformation in region ③

圖12 為剖線CC′的時間序列剖線圖，其中，選擇了6 個時間段顯示，只形成1 個形變漏斗，漏斗中心基本靠近剖線CC′的中間位置。由于剖線的位置靠近30301 工作面，在2021 年8 月和2022 年1 月，剖線上的形變先后發(fā)生兩次突變，第一次突變是受到30301 工作面開采的影響，第二次突變是受到30302工作面開采的影響。圖13 為剖線CC′上的形變最值點c 的時間序列形變折線圖。由圖13 可知，在2021年9 月前后及2022 年3 月前后，c 點發(fā)生了兩次突變，與剖線整體的形變特征基本吻合。

圖12 CC'時間序列剖線圖Fig.12 Time series deformation of CC' profile

圖13 c 點時間序列圖Fig.13 Time series deformation of point c

隨著時間的推移，礦區(qū)內(nèi)主要形變區(qū)域的形變愈加明顯，時序InSAR 監(jiān)測結(jié)果數(shù)據(jù)對比反映了礦區(qū)內(nèi)局部地區(qū)的地表變形情況。從獲取的結(jié)果來看，礦區(qū)的形變區(qū)域是時刻發(fā)生變化的，這是因為礦區(qū)開采造成的形變隨著時間的推移而逐漸演化，可以清楚地看到礦區(qū)形變的緩慢變化過程，礦區(qū)其他區(qū)域地表大部分處于穩(wěn)定狀態(tài)。

3.2 精度評定

精度評定是驗證獲取地表形變量精度的一種方法，一般分為內(nèi)符合精度評定和外符合精度評定。內(nèi)符合精度評定一般是獲取的形變的中誤差（標(biāo)準(zhǔn)差）、形變量分布直方圖或者不同數(shù)據(jù)源之間形變值的相互比較。外符合精度評定一般是借助外部的形變數(shù)據(jù)，如GPS 數(shù)據(jù)或水準(zhǔn)數(shù)據(jù)，與InSAR 獲取的形變數(shù)據(jù)進行對比驗證。

由于本次試驗沒有收集到外部數(shù)據(jù)，且SAR 數(shù)據(jù)源只有Sentinel-1A。因此，通過形變量值的標(biāo)準(zhǔn)差的計算，進行內(nèi)符合精度評定。圖14 為形變速率像元分布圖。由圖14 可知，85.26%的像元的形變速率的絕對值小于10 mm。圖15 為工作區(qū)形變速率標(biāo)準(zhǔn)差分布圖。由圖15 可知，大部分區(qū)域的像元的標(biāo)準(zhǔn)差均在6 mm 范圍內(nèi)，在形變量級大的區(qū)域的像元的標(biāo)準(zhǔn)差相對較大，在6～15 mm 范圍之間。

圖14 形變速率像元分布圖Fig.14 Pixel distribution of deformation rate

圖15 形變速率標(biāo)準(zhǔn)差分布圖Fig.15 Distribution map of standard deviation of deformation rate

3.3 光學(xué)影像解譯

通過無人機正射飛行獲取到區(qū)域①的高分辨率光學(xué)影像，將InSAR 形變與光學(xué)影像疊加到一起進行光學(xué)解譯，如圖16 所示。礦區(qū)地面形變的表現(xiàn)特征一般為地面塌陷造成的地裂縫和塌陷坑等。

圖16 區(qū)域①光學(xué)影像與形變疊加圖Fig.16 Optical image and deformation superposition in region ①

由圖16 可知，區(qū)域①的形變主要為（a）（b）（c）三部分，地表形變最主要的特征表現(xiàn)為地裂縫展布，部分區(qū)域分布有多個塌陷坑，而形變量級較小的區(qū)域地表的光學(xué)特征不甚明顯。結(jié)果顯示，光學(xué)影像特征與InSAR 監(jiān)測的形變結(jié)果的分布較為吻合，但InSAR 技術(shù)更適合小量級的形變監(jiān)測，尤其是光學(xué)影像上無明顯特征的區(qū)域。

4 結(jié)論

本文基于SBAS-InSAR 方法對楊伙盤煤礦在2020 年8 月至2022 年10 月間的地表進行了形變監(jiān)測，并利用高分辨率光學(xué)影像進行了特征解譯。

1）礦區(qū)范圍內(nèi)共分布了三處形變區(qū)域，均呈長條狀分布，與工作面開采范圍高度吻合。整個區(qū)域可探測的最大形變值達(dá)到-130 mm，最大年平均沉陷速率約為-52 mm/a。

2）通過時間序列形變分析，每處形變區(qū)在對應(yīng)或相鄰的工作面開采時發(fā)生突變，持續(xù)時間在1～2個月之間，隨即開始緩慢的持續(xù)性形變。

3）通過高分辨率光學(xué)影像和InSAR 形變結(jié)果的綜合分析，在形變量級較大的區(qū)域，在光學(xué)影像上均展布多條地裂縫，部分區(qū)域還分布有塌陷坑，表明InSAR 監(jiān)測結(jié)果的可靠性。而形變量級較小的區(qū)域，光學(xué)影像上無明顯的特征，表明InSAR 技術(shù)在地表小量級形變監(jiān)測方面的優(yōu)勢。