蔚縣礦區(qū)地面沉陷InSAR 多維形變監(jiān)測

杜建濤，閆 麗，趙超英

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

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

1 研究區(qū)概況

蔚縣位于河北省張家口市最南邊，境內(nèi)南部為高山地區(qū)；中部地勢平坦，壺流河蜿蜒而過；北部為丘陵區(qū)，溝壑密布，地下蘊藏著豐富的煤炭資源。蔚縣礦區(qū)總面積264 km2，距離蔚縣縣城約10 km，是河北省最大的地下采煤區(qū)之一。礦區(qū)井田分布如圖1 所示，崔家寨井田、單侯井田、南留莊井田和西細莊井田為目前主要開采區(qū)域；鄭溝灣礦、興源礦及周邊部分礦井已開采殆盡。受小煤窯破壞，崔家寨井田西翼四采區(qū)目前處于停采狀態(tài)。礦區(qū)含煤地層為侏羅系下花園組，分為上、下兩段，上段為雜色含煤地層，可靠性較差，均不可采；下段含煤8 層，位于底部的1 號煤層(西細莊礦主采層，平均厚度3.16 m)和中部的5 號、6 號煤層(崔家寨礦主采層，平均厚度分別為0.9 m 和2.04 m)，層位穩(wěn)定，厚度和分布面積均較大，為礦區(qū)主采煤層，埋藏深度為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 數(shù)據(jù)采集與處理

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

表1 Sentinel-1A/B SAR 數(shù)據(jù)參數(shù)Table 1 SAR data and their imaging parameters

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

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

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

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

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

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

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

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

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

3 結果分析

3.1 礦區(qū)沉陷監(jiān)測數(shù)據(jù)

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

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

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

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

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

表2 河北蔚縣礦區(qū)沉陷面積統(tǒng)計Table 2 Subsidence area statistics in Yuxian mining area of Hebei Province

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

3.2 二維形變時序

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

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

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

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

圖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 技術獲取了河北蔚縣采礦區(qū)大范圍地面沉陷的空間分布特征，且部分區(qū)域年形變速率達到-20 cm/a。該結果顯示出InSAR 技術用于礦區(qū)地表形變監(jiān)測的優(yōu)越性，為礦區(qū)安全開采和生態(tài)恢復提供科學依據(jù)。

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

c.基于時間函數(shù)模型參數(shù)的采礦沉陷預測和礦區(qū)形變?nèi)S分解將是下一步研究計劃。