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

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

0 引言

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

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

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

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

1 研究區(qū)概況與實驗數(shù)據(jù)

1.1 研究區(qū)概況

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

表1 礦區(qū)基本信息Tab.1 Basic information of mining area

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

1.2 實驗數(shù)據(jù)

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

Δφ=2πfηcΔη

(1)

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

本研究共使用58景研究區(qū)范圍內(nèi)的Sentinel-1A升軌數(shù)據(jù)，影像獲取時間跨度為2015年11月—2018年12月。

2 數(shù)據(jù)處理

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

2.1 數(shù)據(jù)預(yù)處理

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

圖2 技術(shù)路線圖Fig.2 Technique flow chart

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

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

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

2.2 線性形變反演

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

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

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

(2)

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

(3)

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

(4)

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

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

2.3 非線性形變反演

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

垂直方向的沉降結(jié)果可根據(jù)雷達視線方向沉降結(jié)果得到，公式為：

(5)

式中： Δh⊥為垂直方向沉降結(jié)果； Δr為雷達方向沉降結(jié)果?；贕DEMSI5.0軟件可直接得到垂直方向的形變結(jié)果。

3 精度評定

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

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

4 結(jié)果與分析

4.1 礦區(qū)地面沉降速率與累計沉降量分析

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

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

(a) 權(quán)臺煤礦地面沉降速率 (b) 權(quán)臺煤礦累計地面沉降量圖6-1 礦區(qū)地面沉降速率和累計地面沉降量空間分布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 礦區(qū)地面沉降速率和累計地面沉降量空間分布Fig.6-2 Spatial distribution of land subsidence rate and accumulated land subsidence in mining area

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

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

圖7 礦區(qū)內(nèi)不同沉降速率范圍高相干點個數(shù)統(tǒng)計Fig.7 Statistic of high-coherence points in different subsidence rate ranges in the mining area

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

4.2 礦區(qū)地表時序沉降空間分析

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

(a) 2015年11月27日

(b) 2016年6月6日

(c) 2016年12月27日

(d) 2017年6月25日圖9-1 礦區(qū)累計沉降量時序變化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 礦區(qū)累計沉降量時序變化Fig.9-2 Time-series variation diagram of accumulated subsidence in mining area

4.3 礦區(qū)地表沉降時間特征分析

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

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

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

5 結(jié)論

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

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

研究結(jié)果表明，針對礦區(qū)特殊的地質(zhì)條件，基于Sentinel-1A數(shù)據(jù)和MCTSB-InSAR方法，在安全、高效監(jiān)測礦區(qū)沉降方面有較好的應(yīng)用前景，其監(jiān)測結(jié)果可為生態(tài)環(huán)境治理、礦區(qū)沉降預(yù)測、安全隱患區(qū)監(jiān)測等提供數(shù)據(jù)支撐和技術(shù)服務(wù)。