基于Sentinel TOPS模式Stacking技術(shù)監(jiān)測(cè)淮南礦區(qū)沉降

張曉博， 趙學(xué)勝， 葛大慶， 劉斌， 張玲， 李曼， 王艷

(1.防災(zāi)科技學(xué)院生態(tài)環(huán)境學(xué)院，廊坊 101601；2.中國(guó)礦業(yè)大學(xué)(北京) 地球科學(xué)與測(cè)繪 工程學(xué)院，北京 100083；3.中國(guó)國(guó)土資源航空物探遙感中心，北京 100083)

0 引言

當(dāng)巷道上方巖石本身的重量超過(guò)了其最大支撐能力，在采煤過(guò)程中或采煤完成以后容易發(fā)生崩塌。煤礦過(guò)度開采引起地面沉降的空間尺度和量級(jí)大小一般取決于覆蓋深度、上覆地層特征、煤柱尺寸和地形地貌等多種因素[1]。但是，礦區(qū)地面沉陷危害性極強(qiáng)，不僅破壞地面建筑、交通和水利等硬件設(shè)施，同時(shí)還可能危害生態(tài)環(huán)境。因此，地下煤礦開采引起的地面沉降是采煤企業(yè)、政府及社會(huì)各界尤需關(guān)注的問(wèn)題。

在淮南煤礦區(qū)地面沉降監(jiān)測(cè)的相關(guān)研究中，李楠等[2]利用2景PALSAR影像對(duì)該區(qū)進(jìn)行形變監(jiān)測(cè)發(fā)現(xiàn)，采動(dòng)區(qū)均處在活躍階段；張飛[3]和劉文倩[4]利用ENVISAT影像分別監(jiān)測(cè)了淮南礦區(qū)2006年11月—2010年3月和2009年12月—2010年3月期間的地面沉降，結(jié)果表明該區(qū)各個(gè)煤礦開采力度都比較大，所造成的地面塌陷比較明顯。除此之外，未見利用合成孔徑雷達(dá)干涉技術(shù)(interferometric synthetic aperture Radar，InSAR)監(jiān)測(cè)該地區(qū)近期地面沉降的相關(guān)研究。

Stacking技術(shù)基于多景影像解決常規(guī)InSAR面臨的干涉失相干和大氣延遲的問(wèn)題，對(duì)影像數(shù)量要求不太嚴(yán)苛，適用于礦區(qū)沉降速率較快并且存檔數(shù)據(jù)較少的情況[5-6]。針對(duì)Sentinel-1數(shù)據(jù)新型TOPS(terrain observation with progressive scans)成像技術(shù)，本研究基于數(shù)字高程模型(digital elevation model,DEM)配準(zhǔn)和burst重疊區(qū)配準(zhǔn)的方法提高配準(zhǔn)精度，并通過(guò)模擬多項(xiàng)式去除趨勢(shì)性相位，獲取淮南礦區(qū)的沉降速率，分析淮南礦區(qū)的沉降特征，以及Sentinel-1數(shù)據(jù)在該礦區(qū)沉降監(jiān)測(cè)中的適用性和局限性。

1 Sentinel-1數(shù)據(jù)

Sentinel-1是歐空局為哥白尼全球地面觀測(cè)計(jì)劃研制的首個(gè)空間組件，由2顆C波段合成孔徑雷達(dá)(synthetic aperture Radar,SAR)成像極軌衛(wèi)星組成星座[7]。其中，Sentinel-1A于2014年4月發(fā)射，Sentinel-1B于2 a后成功發(fā)射，目前2顆星均可正常工作，雙星重訪周期為6 d。Sentinel-1計(jì)劃的構(gòu)建基于歐洲遙感衛(wèi)星計(jì)劃(European remote-sensing satellite，ERS)和ENVISAT-ASAR傳感器，保證了C波段SAR數(shù)據(jù)的連續(xù)性，其應(yīng)用領(lǐng)域包括土地利用變化、地表形變、水文地質(zhì)和災(zāi)害監(jiān)測(cè)等[8]。

Sentinel-1成像模式有stripmap，IW (interferometric wide swath)，EW (extra wide swath)和wave4種。其中，IW模式是陸地監(jiān)測(cè)的默認(rèn)模式，幅寬250 km，地面空間分辨率5 m×20 m，采用新型TOPS成像技術(shù)。TOPS是一種通過(guò)周期性切換多個(gè)相鄰條帶之間的天線波束獲取數(shù)據(jù)的ScanSAR成像技術(shù)，通過(guò)減弱扇形邊效果提供寬幅和強(qiáng)輻射性影像[9]。TOPS模式旨在取代傳統(tǒng)的ScanSAR模式，既能達(dá)到ScanSAR同樣的影像空間分辨率和覆蓋范圍，又能得到信噪比均一的影像。TOPS技術(shù)除了在距離向控制天線波束，還在方位向上通過(guò)電控天線波束轉(zhuǎn)向避免同一子條帶中不同burst影像質(zhì)量不均勻的問(wèn)題。

2 Stacking技術(shù)

合成孔徑雷達(dá)差分干涉測(cè)量技術(shù)(differential interferometry SAR，DInSAR)是利用衛(wèi)星在不同時(shí)刻通過(guò)同一地區(qū)時(shí)拍攝的2幅影像，根據(jù)衛(wèi)星飛行參數(shù)和相位差計(jì)算cm級(jí)地表形變[10]。隨后發(fā)展起來(lái)的一些基于多景影像計(jì)算形變信息的分析技術(shù)，包括Stacking、永久散射體合成孔徑雷達(dá)干涉測(cè)量技術(shù)(persistent scatterer interferometric SAR，PSInSAR)和短基線集技術(shù)(small baseline subset，SBAS)等，都是通過(guò)組合多個(gè)干涉圖解決常規(guī)InSAR面臨的干涉失相干和大氣延遲的問(wèn)題，以實(shí)現(xiàn)區(qū)域高精度測(cè)量[11-13]。由于PSInSAR技術(shù)需要使用大量數(shù)據(jù)(至少25景影像)估計(jì)和修正大氣延遲，因此，對(duì)影像數(shù)量要求不嚴(yán)苛的Stacking技術(shù)適用性更廣。Stacking技術(shù)通過(guò)生成的多景差分干涉解纏圖估算線性相位速率，實(shí)際是基于最小二乘法(least squares，LS)對(duì)N組觀測(cè)值線性回歸的過(guò)程，估算公式為

(1)

式中：Δtj為第j組干涉圖的時(shí)間基線；φj為第j組解纏的差分干涉相位；ph_rate為線性相位速率。

利用SAR數(shù)據(jù)提取地表形變信息即是差分干涉相位各分量(形變、地形和大氣延遲等)分離的過(guò)程，傳統(tǒng)的Stacking處理流程主要包括影像配準(zhǔn)、干涉對(duì)選擇、干涉圖生成和形變速率估算等步驟。Sentinel-1影像成像方式與以往衛(wèi)星不同，需要對(duì)傳統(tǒng)的處理方法進(jìn)行改進(jìn)，處理流程如圖1所示。

圖1 數(shù)據(jù)處理流程

TOPS影像干涉最大的問(wèn)題是同一burst存在較大的多普勒中心變化。SAR衛(wèi)星側(cè)視成像，聚焦后在方位向和距離向均出現(xiàn)相位跳變的現(xiàn)象，因此失配準(zhǔn)會(huì)導(dǎo)致沿軌道和垂直軌道方向存在相位跳變，然而在同樣的配準(zhǔn)誤差水平下距離向的相位跳變比方位向小，垂直軌道方向的相位跳變可以忽略不計(jì)[14-15]。方位向上TOPS每個(gè)burst都有不同的線性相位跳變，其斜率取決于多普勒中心頻率，方位向失配準(zhǔn)引入的TOPS干涉相位偏差公式為[16]

φerr=2πfdcΔt

(2)

式中：fdc為多普勒中心變化量；Δt為由配準(zhǔn)誤差引入的干涉信號(hào)的時(shí)間偏移。對(duì)Sentinel-1而言，TOPSAR方位向天線掃描的多普勒中心變化量約為5.5 kHz，方位向行掃描時(shí)間約為0.002 s。根據(jù)公式(2)，當(dāng)配準(zhǔn)精度小于0.001個(gè)像元時(shí)，配準(zhǔn)誤差引入的干涉相位偏差小于4°。

為避免不同burst之間會(huì)出現(xiàn)明顯的相位跳變現(xiàn)象，利用三步配準(zhǔn)法進(jìn)行主輔影像高精度配準(zhǔn)。首先基于地形和軌道參數(shù)進(jìn)行初步配準(zhǔn)。影像為三維空間信息在二維平面的展示，成像過(guò)程中會(huì)出現(xiàn)拉伸、扭轉(zhuǎn)等變換，采用單一的多項(xiàng)式進(jìn)行配準(zhǔn)在邊緣誤差較大，因此可利用軌道參數(shù)和地形數(shù)據(jù)提煉主輔影像配準(zhǔn)的查找表。然后，通過(guò)強(qiáng)度互相關(guān)最大化的方法估算距離向和方位向的殘余偏移量，配準(zhǔn)精度達(dá)到0.01個(gè)像元。精確配準(zhǔn)后生成的差分干涉圖出現(xiàn)圖2(a)的現(xiàn)象，即burst之間的條紋變化不連續(xù)，需要進(jìn)一步基于頻譜多樣性進(jìn)行重疊區(qū)配準(zhǔn)，最終得到圖2(b)所示的burst相位跳變消除后的差分干涉圖，一個(gè)色周表示差分干涉相位從0到2π的變化。

(a) 精確配準(zhǔn)后生成的差分干涉圖 (b) burst相位跳變消除后的差分干涉圖

在干涉處理之前，所有影像都添加了精密軌道參數(shù)，但差分干涉圖依然存在殘余相位，該相位包括軌道誤差引入的趨勢(shì)性相位及對(duì)流層和電離層延遲不均勻分布等因素導(dǎo)致的隨機(jī)相位誤差[17]。Stacking技術(shù)通過(guò)多差分條紋圖疊加求平均，將大氣信號(hào)作為隨機(jī)誤差進(jìn)行平差去除大氣影響引入的相位誤差。利用基于二次多項(xiàng)式擬合的方法，減小軌道誤差對(duì)形變結(jié)果精度的影響，模型為

φtrend=a0+a1x+a2y+a3xy+a4x2+a5y2。

(3)

式中：φtrend為趨勢(shì)性相位；a0，a1，…，a5為二次多項(xiàng)式系數(shù)。

基于小基線原則的Stacking技術(shù)，將空間基線閾值設(shè)置為100 m，時(shí)間基線閾值為36 d，得到12個(gè)干涉對(duì)。對(duì)濾波后的差分干涉圖利用最小費(fèi)用流算法進(jìn)行相位解纏，然后基于多景解纏的差分干涉圖根據(jù)公式(1)計(jì)算像元的線性相位速率。

3 淮南礦區(qū)沉降監(jiān)測(cè)

3.1 礦區(qū)概況

淮南礦區(qū)位于安徽省北部，橫跨淮河兩岸，其區(qū)域地質(zhì)構(gòu)造屬華北板塊南緣，東起郯廬斷裂帶，西至阜陽(yáng)斷層，北接蚌埠隆起，南以老人倉(cāng)—壽縣斷層與合肥坳陷相鄰[18]。礦區(qū)有鐵路支線與淮南鐵路(蚌埠—裕溪口)和津浦鐵路接軌，水陸運(yùn)輸都很方便。該礦區(qū)由淮南和潘謝2塊煤田構(gòu)成，煤炭?jī)?chǔ)量豐富，總儲(chǔ)量占安徽省的74%，占華東地區(qū)的50%，且品位優(yōu)良，被譽(yù)為綠色能源。研究區(qū)東西長(zhǎng)約76 km，南北約35 km，分布了顧橋礦、謝橋礦和張集礦等多個(gè)礦井。礦區(qū)在植被茂盛的季節(jié)相干性較低，因此收集了2015年11月3日—2016年3月14日期間的9景影像。

3.2 沉降結(jié)果

淮南礦區(qū)的地面沉降主要受煤礦開采活動(dòng)影響。利用Stacking技術(shù)獲取該礦區(qū)監(jiān)測(cè)期間的沉降速率如圖3所示。圖中沉降速率為負(fù)值表示地面下沉。一些礦井由于長(zhǎng)期開采造成地面下沉在低凹積水，后向散射信息較弱；另外，當(dāng)?shù)V區(qū)沉降中心沉降速率過(guò)大時(shí)，會(huì)出現(xiàn)失相干的現(xiàn)象，導(dǎo)致沉降中心往往只能提取到部分信息。沉降等值面圖是通過(guò)點(diǎn)插值的方法利用相干性較好點(diǎn)的控制編制的。結(jié)果表明，研究區(qū)有多個(gè)明顯的沉降中心，主要分布于西部和北部，沉降速率從邊緣到中心逐漸增加，形成橢圓形或圓形。從圖3可以看出，最為顯著的沉降中心分布在東北方向的潘四井附近，最大沉降速率在80～90 cm/a，該區(qū)段內(nèi)沉降速率大于10 cm/a的區(qū)域約1.82 km2；潘四井東南方向2個(gè)沉降區(qū)相連接，最大沉降速率在70～80 cm/a，沉降速率超過(guò)10 cm/a的區(qū)域約2.35 km2；位于朱集井西部的沉降中心最大沉降速率超過(guò)了70 cm/a，沉降速率大于10 cm/a的區(qū)域約1.54 km2。研究區(qū)南部淮河沿岸出現(xiàn)3個(gè)沉降中心，新莊孜礦最大沉降速率約50 cm/a，沉降速率大于10 cm/a的區(qū)域約2.25 km2；望峰崗井沉降區(qū)的最大沉降速率在40～50 cm/a，2個(gè)連通區(qū)沉降速率大于10 cm/a的面積約2.50 km2；謝橋礦、張集礦、顧橋井、新集礦、潘一礦及潘三礦區(qū)域出現(xiàn)多個(gè)沉降中心，有的沉降區(qū)連成一片，沉降漏斗的邊緣最大沉降速率在40～50 cm/a。將研究區(qū)2 622 km2范圍按照沉降速率進(jìn)行面積統(tǒng)計(jì)，如圖4所示。

圖3研究區(qū)沉降速率

Fig.3Subsidencevelocityoftheresearcharea

圖4 研究區(qū)沉降面積統(tǒng)計(jì)

從圖4可知，2015年11月—2016年3月期間，沉降速率為10～30 cm/a的面積約70.63 km2；在30～60 cm/a之間的沉降面積約10.92 km2；大于60 cm/a的嚴(yán)重沉降區(qū)面積為0.61 km2，約占研究區(qū)總面積的3.13%。研究結(jié)果表明礦區(qū)開采沉陷具有幅度大、沉降范圍小的顯著特點(diǎn)。

煤礦區(qū)地面沉降一般為非線性的過(guò)程，將9景影像按時(shí)間順序配對(duì)成8組干涉對(duì)，分析不同時(shí)間段地下采煤引起的地面沉降特征。圖5為研究區(qū)東北部較嚴(yán)重沉降區(qū)的差分干涉條紋圖。

(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)外，其他干涉圖的時(shí)間基線均為12 d。差分干涉圖中，每個(gè)條紋表示地表視線向移動(dòng)量為2.8 cm。與地面沉降速率圖相同，最大沉降中心位于潘四井，圖5(h)期間的形變量最大達(dá)11.2 cm，而干涉圖5(b)、(c)和(e)期間形變量較小，約5.6 cm。分布在朱集井西部和潘四井東側(cè)的2個(gè)沉降區(qū)，在整個(gè)觀測(cè)周期內(nèi)相干性較好。朱集井西部的沉降區(qū)在2幅時(shí)間基線較長(zhǎng)的干涉圖5(d)和(f)中形變量較大，分別為8.4 cm和5.6 cm，其他時(shí)間段內(nèi)的形變量均在2.8 cm以內(nèi)；潘四井東側(cè)沉降區(qū)在干涉圖5(f)中形變量最大達(dá)11.2 cm，并且出現(xiàn)多個(gè)沉降中心連片的現(xiàn)象。另外，潘一礦和潘三礦附近的部分沉降區(qū)受水體影響，相干性較差，干涉條紋均不完整，僅得到沉降中心邊緣的形變信息。

4 結(jié)論

首次利用新型Sentinel TOPS模式升軌數(shù)據(jù)基于Stacking技術(shù)監(jiān)測(cè)了淮南礦區(qū)地面沉降，在研究過(guò)程中采用三步法精確配準(zhǔn)影像，并利用多項(xiàng)式擬合來(lái)消除差分干涉圖中的趨勢(shì)性相位，最后基于最小二乘法線性回歸得到了研究區(qū)的沉降速率圖。主要結(jié)論為：研究區(qū)出現(xiàn)了多個(gè)沉降中心，主要分布于研究區(qū)的西部和北部，沉降速率空間分布不均一，從邊緣到中心逐漸增加，形成橢圓形或圓形；礦井地面沉降非線性特征顯著，即沉降中心在不同時(shí)間段的形變量不同；Sentinel影像的burst重疊區(qū)配準(zhǔn)和差分干涉圖的趨勢(shì)性相位去除效果較好；時(shí)間基線為36 d時(shí)研究區(qū)內(nèi)有些沉降區(qū)失相干現(xiàn)象較嚴(yán)重，因此礦區(qū)InSAR監(jiān)測(cè)應(yīng)根據(jù)研究區(qū)沉降情況設(shè)置合適的時(shí)間基線閾值。

另外，根據(jù)本文的監(jiān)測(cè)結(jié)果對(duì)基于Sentinel-1影像監(jiān)測(cè)礦區(qū)沉降的適用性和局限性進(jìn)行總結(jié)，包括以下幾點(diǎn)：

1)煤礦區(qū)一般位于山區(qū)或城市郊區(qū)，冬季植被較少相干性較好，而夏季植被茂盛相干目標(biāo)較少，故利用星載SAR數(shù)據(jù)更適合研究冬季礦區(qū)的地面沉降情況。

2)相比于TerraSAR和Radarsat等衛(wèi)星，Sentinel-1雙星重訪周期大大縮短，較高的時(shí)間采樣頻率可降低失相干的影響；影像幅寬較大，空間分辨率適中，利于對(duì)井田范圍較大的礦區(qū)進(jìn)行沉降監(jiān)測(cè)或區(qū)域性的礦區(qū)開采調(diào)查。

3)Stacking技術(shù)主要生成差分干涉圖和沉降速率圖，差分干涉圖的條紋數(shù)量直接與主輔影像獲取時(shí)間段內(nèi)的形變量有關(guān)，條紋越多形變量越大；沉降速率圖是對(duì)多景差分干涉圖解纏相位線性回歸的結(jié)果，可反映礦區(qū)開采沉陷的嚴(yán)重程度及沉陷范圍。

4)形變量結(jié)果為雷達(dá)視線向形變量，這是InSAR技術(shù)固有的局限性，因此利用單一軌道解算沉降速率的前提是假設(shè)視線向形變量?jī)H為垂直向形變量沿視線向的分量，即通過(guò)投影關(guān)系計(jì)算沉降速率。