徐州礦區2016—2018年地面沉降監測與分析

李夢夢， 范雪婷， 陳 超， 李倩楠， 楊 錦

(江蘇省基礎地理信息中心，南京 210013)

0 引言

隨著科學的進步和經濟的發展，地面沉降問題日益突出，已成為21世紀的主要地質災害之一。地面沉降可破壞基礎設施，改變土壤環境和水文環境，引起沉降附近社區的居民恐慌，影響其未來居住規劃等[1-2]。造成地面沉降的原因錯綜復雜，其中固體礦產的開采是影響地面沉降的原因之一[3-4]。被譽為“華東地區的煤炭海洋”的徐州，煤炭資源豐富，但在煤炭開采的同時也造成了地面沉降、地面塌陷等地質環境問題，嚴重制約了該地區的可持續發展[5-6]。傳統的地面沉降監測主要以水準測量和全球衛星導航定位系統(global navigation satellite system，GNSS)測量為主，這些方法雖然單點精度高，但是外業周期長、連續性差、成本較高，且不能獲得整個監測區域的整體沉降趨勢[7]。

礦區本身的地質條件比較特殊，且現場測量安全隱患大，水準作業難度較高，傳統地面沉降監測方法對不定期發生的沉降和塌陷難以做到實時監測[8]，但時間序列合成孔徑雷達干涉測量(interferometric synthetic aperture Radar，InSAR)技術可以安全高效地獲取大范圍、高精度地面沉降數據[9-10]，因此基于InSAR技術來監測礦區地面沉降是較佳選擇。

Liu等[5]基于Sentinel-1影像，監測徐州沛縣部分礦區地面沉降狀況，監測結果表明： ZSL煤礦和SHJ煤礦沉降值與時間呈線性關系，并基于Landsat8影像的分類結果，研究沛縣煤礦的沉降與土地利用分類的關系，發現農村地區比城市地區沉降嚴重，但影像的分類結果還需進一步提高，以便更精確地研究沉降與土地利用分類的關系； 李達等[11]基于小基線集(small baseline subset，SBAS)技術對13景TerraSAR-X影像進行時間序列處理，并基于合成孔徑雷達差分干涉測量(differential InSAR，D-InSAR)技術，驗證SBAS技術得到的監測結果，發現二者監測結果相近，但SBAS技術得到的監測結果更精確，由于缺乏研究區域的實測數據，監測結果的縱向正確性還有待驗證； Chatterjee等[12]基于6對多頻C和L波段的SAR影像，通過D-InSAR技術研究礦區的地面沉降，監測結果表明，C波段可以監測到緩慢或中度沉降的礦區，L波段可以監測到快速沉降的礦區，但影像數量較少，基線長度相對較大； Grzovic等[13]基于時間序列永久散射體干涉測量技術(permanent scatter interferometry，PSI)和SBAS技術處理ERS-1/2和ALOS PALSAR數據，得到Springfield地區煤礦開采區的地面沉降監測結果，由于ALOS PALSAR影像數量較少，其監測結果有失真現象； Bateson等[14]基于間斷相干小基線集(intermittent small baseline subset，ISBAS)技術處理ERS-1/2數據，提高了監測點密度，監測結果表明，礦區存在地面抬升現象，抬升速率達1 cm/a，其原因可能是相較于開采期，停采期礦區水位恢復，但缺乏礦區開采數據，無法準確確定開采初期的礦井水位，其結果分析只能根據假定的開采數據來進行； Gupta等[15]利用InSAR技術獲取研究區域數字高程模型(digital elevation model,DEM)數據來分析礦區地面高程值變化，高程精度偏差為23.07 m，均方根誤差為±2.31 m； 陸燕燕等[16]利用ALOS PALSAR影像對徐州張雙樓煤礦進行地表形變監測，得到2011-01-16—2011-03-03的地表形變分布圖，監測結果表明張雙樓煤礦出現3處沉降漏斗，最大沉降量達420 mm，且沉降漏斗與礦區分布一致，由于實驗數據比較少，只能采用D-InSAR方法來監測礦區沉降變化，并未考慮時空去相關以及大氣非均勻性影響，使得地表形變的監測精度受到限制。

基于Sentinel-1A數據，對在采礦區與停產礦區長時間序列沉降狀況監測與對比分析的研究較少，結合地理國情監測數據研究高相干點分布的研究也比較少，本文利用Sentinel-1A影像數據基于多主影像相干目標小基線InSAR方法(multiple master-image coherent target small-baseline interferometric SAR，MCTSB-InSAR)，結合江蘇省2016年地理國情監測數據，開展徐州部分在采礦區與停采礦區2016—2018年的地面沉降監測研究與對比分析。

1 研究區概況與實驗數據

1.1 研究區概況

徐州市地處蘇魯豫皖四省接壤地區，屬溫帶季風氣候，四季分明，夏無酷暑，冬無嚴寒，是資源富集且組合條件優越的地區之一。其中煤炭儲量約69億t，是江蘇省唯一的煤炭產地，也是江蘇省重要的能源基地。本文根據礦區停采時間，共選擇6個礦區作為研究對象，礦區基本信息見表1，影像、水準點和礦區分布見圖1，圖中白色框為礦區放大影像。

表1 礦區基本信息Tab.1 Basic information of mining area

圖1 影像、水準點和礦區分布Fig.1 Images, bench mark and mining map

1.2 實驗數據

本文采用的實驗數據是Sentinel-1A雷達衛星影像。Sentinel-1是由歐洲委員會和歐洲航天局針對哥白尼全球對地觀測項目研制的一個載有C波段合成孔徑雷達的雙星星座。Sentinel-1的A星于2014年4月3日發射成功，B星于2016年4月25日發射成功，單星重訪周期為12 d，兩星同時工作其重訪周期為6 d。Sentinel-1A衛星有4種工作模式： 條帶模式(strip map，SM)、干涉寬幅模式(interferometric wide swath，IW)、波模式(wave，WV)和超寬幅模式(extra wide swath，EW)，本文使用IW模式的數據。IW模式幅寬為250 km，空間分辨率為5×20 m，采用了TOPS(terrain observation with progressive scans SAR)成像技術。TOPS技術不僅能在保證空間分辨率的情況下增大地面覆蓋范圍，還解決了ScanSAR的圖像不均勻問題，但是該成像模式使影像在方位向的多個脈沖序列(burst)之間易產生多普勒中心不一致問題[17]，引入的干涉相位偏差可表示為：

Δφ=2πfηcΔη

(1)

式中： Δφ為干涉相位偏差；fηc為burst多普勒中心頻率變化值； Δη為方位向主輔影像配準偏移量。干涉測量通常要求干涉相位偏差Δφ不超過3°[18]，則影像配準精度需達到0.001個像素[18-19]。

本研究共使用58景研究區范圍內的Sentinel-1A升軌數據，影像獲取時間跨度為2015年11月—2018年12月。

2 數據處理

基于GDEMSI5.0軟件采用MCTSB-InSAR技術實現對徐州礦區的時間序列監測。MCTSB-InSAR融合了永久散射體InSAR技術(permanent scatterer InSAR，PS-InSAR)、SBAS和相干目標等方法的優勢，針對Sentinel-1A數據的TOPS成像模式特性，數據處理主要包括： 數據預處理、線性形變反演和非線性形變反演3個部分。

2.1 數據預處理

SAR數據預處理主要是SAR影像差分干涉處理的過程。針對Sentinel-1A影像特點，預處理主要包括軌道加密、精確配準、干涉相對選擇、生成差分干涉圖等，技術路線如圖2所示。

圖2 技術路線圖Fig.2 Technique flow chart

增強譜分集配準方法雖然可以達到很高的配準精度，但是有2點缺陷： ①若初始影像之間的配準精度低于0.05個像元，則該方法容易出現迭代死循環的現象，無法滿足配準要求； ②若2幅影像時間基線較長，容易造成時間失相干。又因為基于精密軌道和DEM數據的配準不受影像相干性影響，可與基于相關性配準互補，使粗配準精度達0.05個像元。因此本文采用基于精密軌道數據和DEM數據、相關性配準和遞進式增強譜分集的多級配準方法。其中遞進式增強譜分集就是選擇與其相近的，能夠與第一主影像配準的影像作為第二主影像，配準后再與第一主影像配準，避免時間失相干[20]。

生成差分干涉圖分3步，首先基于典型小區域(一般大小為500像素×500像素)快速生成全部(兩兩組合)小干涉圖； 然后根據小干涉圖質量確定初選干涉像對； 最后根據初選像對生成大區域干涉圖并進行篩選。針對整景影像，先快速生成1 653個小區域差分干涉圖，然后對所有小干涉圖的相干性進行定量評估，按平均相干系數由高到低，篩選出影像數3～5倍的干涉相對，最終保留315幅差分干涉圖，基線對組合如圖3所示。

圖3 小基線對基線時空分布Fig.3 Spatial and temporal distribution of small baseline pair

2.2 線性形變反演

線性形變反演主要包括： 點目標提取、線性形變速率和高程誤差求解等處理。

本文使用幅度離差、平均幅度、平均相干系數“三閾值”法來提取高相干點目標，分別設為0.35，1.1，0.54，本實驗整景影像共提取高相干點目標3 046 739個。

對高相干點目標進行三角網連接，三角網邊上兩頂點之間的二次差分相位差包括線性形變、非線性形變、高程誤差、大氣影響、噪聲5個部分貢獻，具體公式為：

(2)

Δω=ΔRnon-liner+ΔφAPS+Δφnoise

(3)

式中：B⊥和T分別為干涉像對的空間基線和時間基線；λ為波長；R為傳感器到目標的距離；θ為雷達入射角； Δv為線性形變速率； Δh為高程誤差； ΔRnon-liner為非線性形變相位； ΔφAPS為大氣延遲相位； Δφnoise為噪聲相位； Δω為殘余相位。通常|Δω|<π，可通過整體相干性模型來估算Δv和Δh，即

(4)

式中：γ為相干因子值；M為干涉圖數量；i為第i幅干涉圖；j表示單位復數。利用二維周期圖[21]或者空間搜索[22-23]等優化方法，通過模型解算出Δh和Δv。

通常選擇|γ|≥0.7的邊作為可靠的連接，剔除不連接的高相干點目標。最后從某一參考點開始對連接邊上解算出的速率和高程誤差增量進行積分，將得到各相干點上形變速率和高程誤差的絕對量[24]。

2.3 非線性形變反演

從差分干涉相位中去除線性形變相位和高程誤差相位，得到殘余相位，包括大氣相位、非線性形變相位和噪聲相位。根據殘余相位3個分量的不同時空頻譜特征，利用時間和空間濾波的方法將三者分離出來[21]。最后疊加線性形變信息和非線性形變信息得到累計地面沉降量。

垂直方向的沉降結果可根據雷達視線方向沉降結果得到，公式為：

(5)

式中： Δh⊥為垂直方向沉降結果； Δr為雷達方向沉降結果。基于GDEMSI5.0軟件可直接得到垂直方向的形變結果。

3 精度評定

為驗證MCTSB-InSAR方法得到的監測結果，收集江蘇省“十二五”和“十三五”2套水準數據，比較實測水準點得到的年均沉降速率與鄰近InSAR點目標監測值的差異來評估InSAR沉降監測的精度。本文數據處理基于江蘇省域數據處理，監測結果由覆蓋江蘇全省的多個Frame數據拼接得到，數據驗證基于最后的全省拼接結果進行。根據80 m鄰近原則，全省共篩選出35個實測水準點用來對監測結果進行精度評定，即： 2期水準數據的高程差值比上日期間隔，得到基于水準數據的點位沉降速率，然后與該水準點80 m范圍內最近的InSAR點目標的沉降速率對比分析，如圖4所示，年均沉降速率二者差值最大為8.2 mm/a，最小為0.1 mm/a，均方根誤差為4.0 mm/a，差值絕對值小于1 mm/a的水準點占比為40%，滿足規范[25]的要求，保證監測結果的可靠性。

(a) 水準點與InSAR點沉降速率差值統計 (b) 不同差值范圍內點數統計圖4 水準點與InSAR點沉降速率差值及不同差值范圍內點數統計Fig.4 Statistics of subsidence rate difference between level point and InSAR point, and point for different difference ranges

4 結果與分析

4.1 礦區地面沉降速率與累計沉降量分析

監測結果顯示6個礦區都有沉降發生，各個礦區的最大沉降速率和最大累計沉降量如圖5所示，各個礦區的地面沉降速率和累計地面沉降量空間分布如圖6所示。其中權臺煤礦、張集煤礦和拾屯煤礦的最大沉降速率均介于20～30 mm/a之間，旗山煤礦略大于30 mm/a，三河尖煤礦和張雙樓煤礦均大于100 mm/a。權臺煤礦、張集煤礦和拾屯煤礦的最大累計沉降量均不超過100 mm； 旗山煤礦略大于100 mm，三河尖煤礦和張雙樓煤礦均大于300 mm。

(a) 礦區最大沉降速率 (b) 礦區最大累計沉降量圖5 礦區最大沉降速率及最大累計沉降量Fig.5 Maximum subsidence rate and value of mining area

(a) 權臺煤礦地面沉降速率 (b) 權臺煤礦累計地面沉降量圖6-1 礦區地面沉降速率和累計地面沉降量空間分布Fig.6-1 Spatial distribution of land subsidence rate and accumulated land subsidence in mining area

(c) 旗山煤礦地面沉降速率 (d) 旗山煤礦累計地面沉降量

(e) 張集煤礦地面沉降速率 (f) 張集煤礦累計地面沉降量

(g) 拾屯煤礦地面沉降速率 (h) 拾屯煤礦累計地面沉降量

(i) 三河尖煤礦地面沉降速率 (j) 三河尖煤礦累計地面沉降量圖6-2 礦區地面沉降速率和累計地面沉降量空間分布Fig.6-2 Spatial distribution of land subsidence rate and accumulated land subsidence in mining area

(k) 張雙樓煤礦地面沉降速率(l) 張雙樓煤礦累計地面沉降量圖6-3 礦區地面沉降速率和累計地面沉降量空間分布Fig.6-3 Spatial distribution of land subsidence rate and accumulated land subsidence in mining area

統計各個礦區內的高相干點，如圖7所示，結果表明： 權臺煤礦沉降速率大于0 mm/a的點數占比為68.47%。張集煤礦、拾屯煤礦和旗山煤礦沉降速率介于0～10 mm/a的點數占比最大，分別為： 93.90%，69.48%和79.50%。三河尖煤礦沉降速率介于10～30 mm/a的點數占比為92.18%，張雙樓煤礦沉降速率介于30～50 mm/a的點數占比為40.15%，介于10～30 mm/a的點數占比為39.68%。結合江蘇省2016年地理國情監測數據，把地類分為房屋及道路、水體、植被和裸露地表4類，如圖8所示，由于地理國情數據的原因，三河尖煤礦覆蓋區域不完全。從圖8中可看出，三河尖煤礦內房屋及道路面積明顯大于張雙樓煤礦，三河尖煤礦有2 844個高相干點落入房屋及道路內，占該礦區總高相干點數的73.66%，張雙樓煤礦有672個高相干點落入房屋及道路內，占該礦區總高相干點數的63.33%。由此也說明，高相干點目標提取較為可靠。

圖7 礦區內不同沉降速率范圍高相干點個數統計Fig.7 Statistic of high-coherence points in different subsidence rate ranges in the mining area

(a) 三河尖煤礦(b) 張雙樓煤礦圖8 礦區土地利用分類及高相干點分布Fig.8 Land use classification and high-coherence distribution map of mining area

4.2 礦區地表時序沉降空間分析

以2015年11月27日為起始時間，其他時間相對起始時間的研究礦區時序累計沉降量空間分布如圖9所示。從圖9中可以看出，從起始時間到2016年6月期間，張雙樓煤礦、三河尖煤礦、旗山煤礦、拾屯煤礦有10 mm以上的沉降量，張集煤礦、權臺煤礦沒有發生明顯沉降現象，直到2017年6月，權臺煤礦開始出現明顯沉降現象，張集煤礦稍有沉降現象。隨著時間推移，張雙樓煤礦和三河尖煤礦的沉降范圍迅速擴大，沉降逐漸加深，對周邊影響也比較嚴重。其余煤礦的沉降均發生在礦區范圍內，且沉降范圍較穩定，沒有明顯擴張趨勢，對其周邊地區影響較小。總體來說，礦區在停采一段時間后沉降仍然會繼續發生，不同停采礦區的沉降情況不同，但是停采礦區的沉降情況比在采礦區要輕微得多。

(a) 2015年11月27日

(b) 2016年6月6日

(c) 2016年12月27日

(d) 2017年6月25日圖9-1 礦區累計沉降量時序變化Fig.9-1 Time-series variation diagram of accumulated subsidence in mining area

(e) 2017年12月22日

(f) 2018年7月2日

(g) 2018年12月5日圖9-2 礦區累計沉降量時序變化Fig.9-2 Time-series variation diagram of accumulated subsidence in mining area

4.3 礦區地表沉降時間特征分析

為分析研究區域的沉降時間特征，提取各礦區最大沉降點的時序沉降值，利用簡單的線性擬合模型，得到沉降值和時間的關系，如圖10所示。從圖10中可以看出，礦區的沉降量與時間基本都呈線性關系，R2均大于0.9。線性擬合斜率絕對值從大到小為： 張雙樓煤礦(0.365)、三河尖煤礦(0.277)、旗山煤礦(0.097)、拾屯煤礦(0.065)、權臺煤礦(0.057)、張集煤礦(0.056)，張雙樓煤礦和三河尖煤礦的線性規律比其他煤礦要強，權臺煤礦的線性規律較弱。對權臺煤礦的時序沉降量做簡單的非線性擬合，發現其沉降規律更符合非線性變化，在監測末期，沉降接近平穩狀態。雖然旗山煤礦的停采時間與權臺煤礦較一致，但旗山煤礦的時序沉降量沒有表現出明顯的非線性變化，由于影響礦區沉降的因素多種多樣，停采礦區的沉降規律還有待進一步研究。

(a) 權臺煤礦(b) 旗山煤礦(c) 張集煤礦

(d) 拾屯煤礦(e) 三河尖煤礦(f) 張雙樓煤礦圖10 礦區累計沉降量變化統計Fig.10 Statistical chart of accumulative subsidence in mining are

5 結論

本文基于Sentinel-1A雷達影像數據利用MCTSB-InSAR方法監測徐州部分礦區2016—2018年間地面沉降狀況，經實測水準數據驗證，二者差值的均方根誤差為4.0 mm/a，符合技術規范要求，監測結果較準確、可靠。

監測結果表明，2016—2018年間張雙樓煤礦沉降最為嚴重，其次為三河尖煤礦，二者的最大年均沉降速率均超過100 mm/a，最大累計沉降量均超過300 mm。結合江蘇省2016年的地理國情監測數據發現，本研究提取的高相干點目標較大部分對應房屋、道路等人工地物，揭示了人工地物具有高相干特性。通過分析沉降量的時序特征發現，除權臺煤礦外，其他礦區的沉降有明顯線性規律，在采礦區的線性規律比停采礦區的線性規律要強一些，權臺煤礦呈非線性沉降，具體沉降規律還有待進一步研究。

研究結果表明，針對礦區特殊的地質條件，基于Sentinel-1A數據和MCTSB-InSAR方法，在安全、高效監測礦區沉降方面有較好的應用前景，其監測結果可為生態環境治理、礦區沉降預測、安全隱患區監測等提供數據支撐和技術服務。