張曉博, 趙學勝, 葛大慶, 劉斌, 張玲, 李曼, 王艷
(1.防災科技學院生態環境學院,廊坊 101601;2.中國礦業大學(北京) 地球科學與測繪 工程學院,北京 100083;3.中國國土資源航空物探遙感中心,北京 100083)
當巷道上方巖石本身的重量超過了其最大支撐能力,在采煤過程中或采煤完成以后容易發生崩塌。煤礦過度開采引起地面沉降的空間尺度和量級大小一般取決于覆蓋深度、上覆地層特征、煤柱尺寸和地形地貌等多種因素[1]。但是,礦區地面沉陷危害性極強,不僅破壞地面建筑、交通和水利等硬件設施,同時還可能危害生態環境。因此,地下煤礦開采引起的地面沉降是采煤企業、政府及社會各界尤需關注的問題。
在淮南煤礦區地面沉降監測的相關研究中,李楠等[2]利用2景PALSAR影像對該區進行形變監測發現,采動區均處在活躍階段;張飛[3]和劉文倩[4]利用ENVISAT影像分別監測了淮南礦區2006年11月—2010年3月和2009年12月—2010年3月期間的地面沉降,結果表明該區各個煤礦開采力度都比較大,所造成的地面塌陷比較明顯。除此之外,未見利用合成孔徑雷達干涉技術(interferometric synthetic aperture Radar,InSAR)監測該地區近期地面沉降的相關研究。
Stacking技術基于多景影像解決常規InSAR面臨的干涉失相干和大氣延遲的問題,對影像數量要求不太嚴苛,適用于礦區沉降速率較快并且存檔數據較少的情況[5-6]。針對Sentinel-1數據新型TOPS(terrain observation with progressive scans)成像技術,本研究基于數字高程模型(digital elevation model,DEM)配準和burst重疊區配準的方法提高配準精度,并通過模擬多項式去除趨勢性相位,獲取淮南礦區的沉降速率,分析淮南礦區的沉降特征,以及Sentinel-1數據在該礦區沉降監測中的適用性和局限性。
Sentinel-1是歐空局為哥白尼全球地面觀測計劃研制的首個空間組件,由2顆C波段合成孔徑雷達(synthetic aperture Radar,SAR)成像極軌衛星組成星座[7]。其中,Sentinel-1A于2014年4月發射,Sentinel-1B于2 a后成功發射,目前2顆星均可正常工作,雙星重訪周期為6 d。Sentinel-1計劃的構建基于歐洲遙感衛星計劃(European remote-sensing satellite,ERS)和ENVISAT-ASAR傳感器,保證了C波段SAR數據的連續性,其應用領域包括土地利用變化、地表形變、水文地質和災害監測等[8]。
Sentinel-1成像模式有stripmap,IW (interferometric wide swath),EW (extra wide swath)和wave4種。其中,IW模式是陸地監測的默認模式,幅寬250 km,地面空間分辨率5 m×20 m,采用新型TOPS成像技術。TOPS是一種通過周期性切換多個相鄰條帶之間的天線波束獲取數據的ScanSAR成像技術,通過減弱扇形邊效果提供寬幅和強輻射性影像[9]。TOPS模式旨在取代傳統的ScanSAR模式,既能達到ScanSAR同樣的影像空間分辨率和覆蓋范圍,又能得到信噪比均一的影像。TOPS技術除了在距離向控制天線波束,還在方位向上通過電控天線波束轉向避免同一子條帶中不同burst影像質量不均勻的問題。
合成孔徑雷達差分干涉測量技術(differential interferometry SAR,DInSAR)是利用衛星在不同時刻通過同一地區時拍攝的2幅影像,根據衛星飛行參數和相位差計算cm級地表形變[10]。隨后發展起來的一些基于多景影像計算形變信息的分析技術,包括Stacking、永久散射體合成孔徑雷達干涉測量技術(persistent scatterer interferometric SAR,PSInSAR)和短基線集技術(small baseline subset,SBAS)等,都是通過組合多個干涉圖解決常規InSAR面臨的干涉失相干和大氣延遲的問題,以實現區域高精度測量[11-13]。由于PSInSAR技術需要使用大量數據(至少25景影像)估計和修正大氣延遲,因此,對影像數量要求不嚴苛的Stacking技術適用性更廣。Stacking技術通過生成的多景差分干涉解纏圖估算線性相位速率,實際是基于最小二乘法(least squares,LS)對N組觀測值線性回歸的過程,估算公式為
(1)
式中:Δtj為第j組干涉圖的時間基線;φj為第j組解纏的差分干涉相位;ph_rate為線性相位速率。
利用SAR數據提取地表形變信息即是差分干涉相位各分量(形變、地形和大氣延遲等)分離的過程,傳統的Stacking處理流程主要包括影像配準、干涉對選擇、干涉圖生成和形變速率估算等步驟。Sentinel-1影像成像方式與以往衛星不同,需要對傳統的處理方法進行改進,處理流程如圖1所示。

圖1 數據處理流程
TOPS影像干涉最大的問題是同一burst存在較大的多普勒中心變化。SAR衛星側視成像,聚焦后在方位向和距離向均出現相位跳變的現象,因此失配準會導致沿軌道和垂直軌道方向存在相位跳變,然而在同樣的配準誤差水平下距離向的相位跳變比方位向小,垂直軌道方向的相位跳變可以忽略不計[14-15]。方位向上TOPS每個burst都有不同的線性相位跳變,其斜率取決于多普勒中心頻率,方位向失配準引入的TOPS干涉相位偏差公式為[16]
φerr=2πfdcΔt
(2)
式中:fdc為多普勒中心變化量;Δt為由配準誤差引入的干涉信號的時間偏移。對Sentinel-1而言,TOPSAR方位向天線掃描的多普勒中心變化量約為5.5 kHz,方位向行掃描時間約為0.002 s。根據公式(2),當配準精度小于0.001個像元時,配準誤差引入的干涉相位偏差小于4°。
為避免不同burst之間會出現明顯的相位跳變現象,利用三步配準法進行主輔影像高精度配準。首先基于地形和軌道參數進行初步配準。影像為三維空間信息在二維平面的展示,成像過程中會出現拉伸、扭轉等變換,采用單一的多項式進行配準在邊緣誤差較大,因此可利用軌道參數和地形數據提煉主輔影像配準的查找表。然后,通過強度互相關最大化的方法估算距離向和方位向的殘余偏移量,配準精度達到0.01個像元。精確配準后生成的差分干涉圖出現圖2(a)的現象,即burst之間的條紋變化不連續,需要進一步基于頻譜多樣性進行重疊區配準,最終得到圖2(b)所示的burst相位跳變消除后的差分干涉圖,一個色周表示差分干涉相位從0到2π的變化。

(a) 精確配準后生成的差分干涉圖 (b) burst相位跳變消除后的差分干涉圖
在干涉處理之前,所有影像都添加了精密軌道參數,但差分干涉圖依然存在殘余相位,該相位包括軌道誤差引入的趨勢性相位及對流層和電離層延遲不均勻分布等因素導致的隨機相位誤差[17]。Stacking技術通過多差分條紋圖疊加求平均,將大氣信號作為隨機誤差進行平差去除大氣影響引入的相位誤差。利用基于二次多項式擬合的方法,減小軌道誤差對形變結果精度的影響,模型為
φtrend=a0+a1x+a2y+a3xy+a4x2+a5y2。
(3)
式中:φtrend為趨勢性相位;a0,a1,…,a5為二次多項式系數。
基于小基線原則的Stacking技術,將空間基線閾值設置為100 m,時間基線閾值為36 d,得到12個干涉對。對濾波后的差分干涉圖利用最小費用流算法進行相位解纏,然后基于多景解纏的差分干涉圖根據公式(1)計算像元的線性相位速率。
淮南礦區位于安徽省北部,橫跨淮河兩岸,其區域地質構造屬華北板塊南緣,東起郯廬斷裂帶,西至阜陽斷層,北接蚌埠隆起,南以老人倉—壽縣斷層與合肥坳陷相鄰[18]。礦區有鐵路支線與淮南鐵路(蚌埠—裕溪口)和津浦鐵路接軌,水陸運輸都很方便。該礦區由淮南和潘謝2塊煤田構成,煤炭儲量豐富,總儲量占安徽省的74%,占華東地區的50%,且品位優良,被譽為綠色能源。研究區東西長約76 km,南北約35 km,分布了顧橋礦、謝橋礦和張集礦等多個礦井。礦區在植被茂盛的季節相干性較低,因此收集了2015年11月3日—2016年3月14日期間的9景影像。
淮南礦區的地面沉降主要受煤礦開采活動影響。利用Stacking技術獲取該礦區監測期間的沉降速率如圖3所示。圖中沉降速率為負值表示地面下沉。一些礦井由于長期開采造成地面下沉在低凹積水,后向散射信息較弱;另外,當礦區沉降中心沉降速率過大時,會出現失相干的現象,導致沉降中心往往只能提取到部分信息。沉降等值面圖是通過點插值的方法利用相干性較好點的控制編制的。結果表明,研究區有多個明顯的沉降中心,主要分布于西部和北部,沉降速率從邊緣到中心逐漸增加,形成橢圓形或圓形。從圖3可以看出,最為顯著的沉降中心分布在東北方向的潘四井附近,最大沉降速率在80~90 cm/a,該區段內沉降速率大于10 cm/a的區域約1.82 km2;潘四井東南方向2個沉降區相連接,最大沉降速率在70~80 cm/a,沉降速率超過10 cm/a的區域約2.35 km2;位于朱集井西部的沉降中心最大沉降速率超過了70 cm/a,沉降速率大于10 cm/a的區域約1.54 km2。研究區南部淮河沿岸出現3個沉降中心,新莊孜礦最大沉降速率約50 cm/a,沉降速率大于10 cm/a的區域約2.25 km2;望峰崗井沉降區的最大沉降速率在40~50 cm/a,2個連通區沉降速率大于10 cm/a的面積約2.50 km2;謝橋礦、張集礦、顧橋井、新集礦、潘一礦及潘三礦區域出現多個沉降中心,有的沉降區連成一片,沉降漏斗的邊緣最大沉降速率在40~50 cm/a。將研究區2 622 km2范圍按照沉降速率進行面積統計,如圖4所示。

圖3研究區沉降速率
Fig.3Subsidencevelocityoftheresearcharea

圖4 研究區沉降面積統計
從圖4可知,2015年11月—2016年3月期間,沉降速率為10~30 cm/a的面積約70.63 km2;在30~60 cm/a之間的沉降面積約10.92 km2;大于60 cm/a的嚴重沉降區面積為0.61 km2,約占研究區總面積的3.13%。研究結果表明礦區開采沉陷具有幅度大、沉降范圍小的顯著特點。
煤礦區地面沉降一般為非線性的過程,將9景影像按時間順序配對成8組干涉對,分析不同時間段地下采煤引起的地面沉降特征。圖5為研究區東北部較嚴重沉降區的差分干涉條紋圖。

(a) 2015/11/03—2015/11/15 (b) 2015/11/15—2015/11/27 (c) 2015/11/27—2015/12/09

(d) 2015/12/09—2016/01/14 (e) 2016/01/14—2016/01/26 (f) 2016/01/26—2016/02/19

(g) 2016/02/19—2016/03/02 (h) 2016/03/02—2016/03/14
除干涉圖5(d)和(f)外,其他干涉圖的時間基線均為12 d。差分干涉圖中,每個條紋表示地表視線向移動量為2.8 cm。與地面沉降速率圖相同,最大沉降中心位于潘四井,圖5(h)期間的形變量最大達11.2 cm,而干涉圖5(b)、(c)和(e)期間形變量較小,約5.6 cm。分布在朱集井西部和潘四井東側的2個沉降區,在整個觀測周期內相干性較好。朱集井西部的沉降區在2幅時間基線較長的干涉圖5(d)和(f)中形變量較大,分別為8.4 cm和5.6 cm,其他時間段內的形變量均在2.8 cm以內;潘四井東側沉降區在干涉圖5(f)中形變量最大達11.2 cm,并且出現多個沉降中心連片的現象。另外,潘一礦和潘三礦附近的部分沉降區受水體影響,相干性較差,干涉條紋均不完整,僅得到沉降中心邊緣的形變信息。
首次利用新型Sentinel TOPS模式升軌數據基于Stacking技術監測了淮南礦區地面沉降,在研究過程中采用三步法精確配準影像,并利用多項式擬合來消除差分干涉圖中的趨勢性相位,最后基于最小二乘法線性回歸得到了研究區的沉降速率圖。主要結論為:研究區出現了多個沉降中心,主要分布于研究區的西部和北部,沉降速率空間分布不均一,從邊緣到中心逐漸增加,形成橢圓形或圓形;礦井地面沉降非線性特征顯著,即沉降中心在不同時間段的形變量不同;Sentinel影像的burst重疊區配準和差分干涉圖的趨勢性相位去除效果較好;時間基線為36 d時研究區內有些沉降區失相干現象較嚴重,因此礦區InSAR監測應根據研究區沉降情況設置合適的時間基線閾值。
另外,根據本文的監測結果對基于Sentinel-1影像監測礦區沉降的適用性和局限性進行總結,包括以下幾點:
1)煤礦區一般位于山區或城市郊區,冬季植被較少相干性較好,而夏季植被茂盛相干目標較少,故利用星載SAR數據更適合研究冬季礦區的地面沉降情況。
2)相比于TerraSAR和Radarsat等衛星,Sentinel-1雙星重訪周期大大縮短,較高的時間采樣頻率可降低失相干的影響;影像幅寬較大,空間分辨率適中,利于對井田范圍較大的礦區進行沉降監測或區域性的礦區開采調查。
3)Stacking技術主要生成差分干涉圖和沉降速率圖,差分干涉圖的條紋數量直接與主輔影像獲取時間段內的形變量有關,條紋越多形變量越大;沉降速率圖是對多景差分干涉圖解纏相位線性回歸的結果,可反映礦區開采沉陷的嚴重程度及沉陷范圍。
4)形變量結果為雷達視線向形變量,這是InSAR技術固有的局限性,因此利用單一軌道解算沉降速率的前提是假設視線向形變量僅為垂直向形變量沿視線向的分量,即通過投影關系計算沉降速率。