蔚縣礦區地面沉陷InSAR 多維形變監測

杜建濤，閆 麗，趙超英

(1.長安大學地質工程與測繪學院，陜西 西安 710054；2.中煤科工集團西安研究院有限公司，陜西 西安 710077)

我國是一個煤炭能源大國，已探明的煤儲量位居世界前列，而礦產資源的過度開采容易破壞礦區地質結構，進而引發地表開裂、塌陷甚至滑坡等地質災害[1]。因此，礦區沉陷監測是煤礦安全開采和可持續發展的重要工作，對預防潛在地質災害具有重要意義[2]。合成孔徑雷達干涉測量(Interferometric Synthetic Aperture Rader，InSAR)因其高精度、大覆蓋范圍及可重復觀測等技術優勢在采礦沉陷監測中得到廣泛應用，極大地克服了傳統測量手段空間分辨率低、周期長、成本高等局限。Zhao Chaoying等[3]、Yang Chengsheng 等[4]對山西大同礦區大量地表塌陷進行InSAR 監測；Zheng Meinan 等[5]、Du Sen等[6]基于多種先進InSAR 技術分析了河北峰峰礦區的沉陷機理；Yang Zefa 等[7]、Ma Chao 等[8]、賀衛中等[9]針對陜西榆林礦區地面塌陷采用InSAR 技術進行大范圍、多維形變監測。然而，目前針對采礦沉陷的InSAR 時空演化分析和精細化監測研究不足，且該技術在河北蔚縣礦區沉陷監測的應用也較少。為了揭示該礦區的地表形變演化特征，本文利用Sentinel-1A/B 數據對采礦沉陷開展InSAR 監測，獲取了礦區升降軌數據沿視線向(LOS)形變監測結果，并在此基礎上重點分析了西細莊礦近1 a 內地表二維形變的時空演化特征。

1 研究區概況

蔚縣位于河北省張家口市最南邊，境內南部為高山地區；中部地勢平坦，壺流河蜿蜒而過；北部為丘陵區，溝壑密布，地下蘊藏著豐富的煤炭資源。蔚縣礦區總面積264 km2，距離蔚縣縣城約10 km，是河北省最大的地下采煤區之一。礦區井田分布如圖1 所示，崔家寨井田、單侯井田、南留莊井田和西細莊井田為目前主要開采區域；鄭溝灣礦、興源礦及周邊部分礦井已開采殆盡。受小煤窯破壞，崔家寨井田西翼四采區目前處于停采狀態。礦區含煤地層為侏羅系下花園組，分為上、下兩段，上段為雜色含煤地層，可靠性較差，均不可采；下段含煤8 層，位于底部的1 號煤層(西細莊礦主采層，平均厚度3.16 m)和中部的5 號、6 號煤層(崔家寨礦主采層，平均厚度分別為0.9 m 和2.04 m)，層位穩定，厚度和分布面積均較大，為礦區主采煤層，埋藏深度為143.6～336.5 m。煤層頂板巖性為泥巖、粉砂巖、細砂巖，底板巖性為粉砂巖、細砂巖及泥巖。地表第四系松散層厚140～245 m，基巖厚度30～55 m，松散層下段砂礫石含水層厚7.5～34.3 m，巖性為粗砂、礫石及礫卵石[10]。

圖1 河北蔚縣地理位置及SAR 影像覆蓋范圍Fig.1 Location of Yuxian County and SAR image coverage

2 數據采集與處理

本文采用61景Sentinel-1A/B干涉寬幅(Interferometric Wide swath，IW)模式數據進行礦區形變監測，包括33 景降軌數據(軌道號T120，時間為2017 年7 月至2018 年9 月)和28 景升軌數據(軌道號T142，時間為2017 年8 月至2018 年8 月)，影像具體參數如表1 所示。此外，本文還利用精密軌道數據AUX_POEORB文件來精化影像初始軌道狀態矢量，利用30 m 分辨率的SRTM DEM 數據來模擬和去除地形相位。

表1 Sentinel-1A/B SAR 數據參數Table 1 SAR data and their imaging parameters

對獲取的61 景升降軌Sentinel-1A/B 數據分別進行小基線集(Small Baseline Subsets，SBAS)干涉處理[11]，然后采用融合升降軌數據的多維形變時序估計方法，生成垂向和東西向的二維形變時間序列[12]。關鍵處理步驟如下：

a.Sentinel-1 預處理 由于Sentinel-1 數據的特殊成像模式，首先需對原始影像進行預處理，將其轉化成可以干涉的單視復數(SLC)格式，具體包括采用AUX_POEORB 軌道文件進行基線精化，對主、從影像進行去斜和配準，使方位向配準精度達到千分之一像元。

b.SLC 影像干涉 設置時間和空間基線閾值生成小基線差分干涉圖。為提高干涉圖的相干性，并抑制噪聲，對差分干涉圖進行自適應頻率域濾波處理，采用基于Delaunay 三角網的最小費用流算法進行相位解纏。為抑制時空去相干的影像，本文時間基線閾值取36 d，空間基線閾值取200 m。干涉組合的時空基線分布如圖2 所示。

圖2 干涉組合時空基線分布Fig.2 Baseline of InSAR combination

c.生成形變圖 挑選高質量干涉組合，估計并去除趨勢項和大氣延遲相位，構建相干像元的形變相位解算模型，采用最小二乘法進行參數求解，得到SAR 影像獲取時刻間的平均形變速率。

d.計算二維形變時間序列 融合多軌道SAR 數據的多維形變時序估計技術，利用奇異值分解法來計算多源SAR 數據在重疊時間段內的二維形變速率和形變時間序列。為了避免觀測方程的病態，需要對設計矩陣進行正則化處理[13]。觀測方程表示為式(1)。

式中A=(-cosαsinθΔt,cosθΔt)為設計矩陣；α為衛星飛行方向的方位角；θ為雷達波入射角；Δt為相鄰SAR影像時間間隔；ξ為正則化參數；I為單位矩陣；V=(VE,VU)T，VE和VU為待求參數，分別表示相鄰SAR 影像獲取時刻的東西向和垂向形變速率；D為形變相位矩陣。待求參數估計值可表示為：

最后，對相鄰SAR 獲取時刻間的形變速率在時間域積分，即可得到地表東西向和垂向的形變時間序列：

式中d為累積形變時間序列；n為不同軌道影像時間跨度重合的總數量。

3 結果分析

3.1 礦區沉陷監測數據

將獲取的升降軌SAR 數據按上述流程處理，生成蔚縣礦區LOS 向年平均形變速率圖，重采樣到地理坐標系下，如圖3 所示，其中，圖3a 為降軌形變監測結果，圖3b 為升軌形變監測結果。

圖3 河北蔚縣礦區2017—2018 年地表形變場Fig.3 Surface deformation site of Yuxian mining area of Hebei Province from 2017 to 2018

從圖3 可以看出，升、降軌數據獲取的礦區形變場基本一致，整個礦區在監測期間均出現不同程度的地表形變現象，且沉陷區普遍呈近橢圓形的漏斗狀分布。礦區崔家寨井田、單侯井田和西細莊井田內采煤造成的地面沉陷較為嚴重，南留莊井田在監測期間形變相對不明顯。其中，崔家寨井田沉陷面積最大，且沉降速率達-30 cm/a。資料顯示，礦區在2015 年調整井田邊界和優化資源配置，提高了各礦井生產能力，其中，崔家寨煤礦生產能力180 萬t/a，多個立井開拓；單侯煤礦生產能力150 萬t/a，3 個立井開拓；西細莊礦生產能力60 萬t/a，一對立井開拓；南留莊礦生產能力30 萬t/a，一對立井開拓[14]。地下煤炭資源的長期高強度開采破壞開采區周圍原始應力平衡，從而導致覆巖移動并向上波及地表，使地表發生形變和沉陷。

為了定量分析礦區地面沉陷災害情況，本文對不同量級年沉陷速率區域的面積進行統計分析。圖4 為礦區沉陷等值線圖，等值線由內向外依次為-20、-15、-10、-5 cm/a，從圖4 可以看出，崔家寨、麥子破、回回墓等多個村莊受到采礦沉陷影響。表2 為不同量級年沉陷速率區域面積統計結果，其中年沉陷速率超過-10 cm/a 的區域達到2.16 km2，年沉陷速率超過-20 cm/a 的區域達到0.10 km2。大面積的地表沉陷對礦區的基礎建設及地表生態環境造成了嚴重的影響[15]。

圖4 河北蔚縣礦區沉陷等值線(單位：cm/a)Fig.4 Subsidence contours in Yuxian mining area in Hebei Province

表2 河北蔚縣礦區沉陷面積統計Table 2 Subsidence area statistics in Yuxian mining area of Hebei Province

對比圖3a、圖3b 可以發現，升、降軌數據獲取的監測結果中部分區域在形變范圍和量級上存在一些差異。為此，本文提取跨越西細莊礦塌陷區的剖線AA′，對其形變空間特征進行分析(圖5)。從圖5 可以看出，升、降軌數據獲取的西細莊礦最大年形變速率相差約5 cm/a，且沉陷中心位置也具有一定偏差。究其原因，可能是由于升、降軌數據具有不同的軌道方向和入射角度，其對地表真實形變的敏感度不同。另外，采礦沉陷的LOS 向形變中可能包含水平形變的影響。

3.2 二維形變時序

圖5 蔚縣西細莊煤礦年形變速率剖線Fig.5 Annual deformation rate along profile A—A′ of Xixizhuang mine in Yuxian County

由于西細莊井田的干涉圖整體相干性較高，便于時序分析，因此，在升、降軌一維形變分析基礎上，仍以該礦為試驗區，采用融合多軌道SAR 數據的多維形變估算法，研究礦區二維形變時空演化特征。將統一到相同地理坐標系下的升降軌形變圖，按式(1)構建方程組，求解不同時間段內的二維形變速率，然后依據式(3)，獲取西細莊礦2017 年8月至2018 年8 月期間垂向與東西向的形變時間序列，如圖6 所示。分析圖6 可以發現，西細莊礦地表在監測期間以垂向形變為主，也具有明顯的東西向水平形變，且東、西水平形變方向相反，量級也存在一定差異。垂向形變在空間上呈規則的橢圓形，東西向水平形變則表現為兩個半圓形，二者形變范圍均隨時間逐漸擴大，形變量級也相應增加。

為定量分析礦區地表形變的時空演化特征，共選取4 個特征點P1、P2、P3、P4(圖6)，進行二維形變時間序列分析(圖7)，其中P1、P2 點分別位于沉陷中心區域及采礦工作面傾向邊緣，P3、P4 點位于采礦工作面走向邊緣。

由圖7 可以看出，在2017 年8 月至2018 年8月期間，P1 點的垂向累積形變最大，達–30 cm，且具有持續下沉的趨勢；P3、P4 點位于東、西向水平形變最大的區域，形變量分別達到–10 cm 和7 cm，其形變方向相反，其中正值表示向東移動，負值表示向西移動，且形變都指向沉陷漏斗中心，同時二者均具有較大的垂向形變，形變量分別達–22 cm 和–13 cm?？梢钥闯?，P3 點的垂向和東西向形變量都較點P4 大，這與工作面的開采方向有關，P2 點的垂向形變達–20 cm，但和P1 點都具有水平方向保持不變的特點。此外，不同于城市地面沉降的線性特征，采礦沉降的形變過程與時間呈現出非線性關系，從圖中P1 點的垂向形變時間序列曲線(圖7a)可以看出，該點處的地表形變具有先緩慢下沉后加速沉降最后趨于穩定的特征，這一形變演化過程對應于礦區地表單點沉降的初始沉降、主要沉降和殘余沉降3 個時期，與“S”型增長的Logistic 時間函數模型能夠很好地吻合[16-17]，更能反映采礦沉陷的真實規律，因此，利用InSAR 數據進行多維形變估算獲得的地表沉陷值，與實際情況較為吻合，可為沉陷的分析、治理提供參考數據。

圖6 蔚縣西細莊礦二維形變時間序列Fig.6 Time sequence of 2D deformation of Xixizhuang mine in Yuxian County

圖7 蔚縣西細莊礦特征點二維形變時間序列Fig.7 Time sequence of 2D deformation of feature points of Xixizhuang mine in Yuxian County

4 結論

a.采用InSAR 技術獲取了河北蔚縣采礦區大范圍地面沉陷的空間分布特征，且部分區域年形變速率達到-20 cm/a。該結果顯示出InSAR 技術用于礦區地表形變監測的優越性，為礦區安全開采和生態恢復提供科學依據。

b.通過對研究區西細莊井田地表形變的精細化監測與分析，獲得采礦沉陷的時空演化特征，即礦區地表動態沉降過程符合Logistic 時間函數模型，該特征揭示了真實的采礦沉陷規律。

c.基于時間函數模型參數的采礦沉陷預測和礦區形變三維分解將是下一步研究計劃。