張家口明長城景觀廊道Sentinel-1影像SBAS形變監(jiān)測(cè)示范研究

何海英，陳彩芬，陳富龍，唐攀攀

(1.中國科學(xué)院空天信息創(chuàng)新研究院，北京 100094； 2.中國科學(xué)院大學(xué)，北京 101408； 3.北京聚才振邦企業(yè)管理顧問有限公司，北京 100038)

0 引言

張家口明長城墻體遺產(chǎn)裸露于地表并可受自然與人類過程諸多因素影響[1]。盡管長城沿線地區(qū)在一定周期和范圍內(nèi)未發(fā)生地震等大型自然災(zāi)害，長時(shí)間緩慢地表形變?nèi)钥沙蔀橛绊憦埣铱诿鏖L城文化景觀整體穩(wěn)定性的重要因素之一?？紤]到張家口明長城部分墻段緊鄰采礦工業(yè)區(qū)和自然不穩(wěn)定坡體，加之冬奧會(huì)修建場(chǎng)館的影響(例如，位于崇禮區(qū)東部的北歐中心跳臺(tái)滑雪場(chǎng)和越野滑雪場(chǎng)選址規(guī)劃緊鄰明長城遺址)； 因此亟須以微形變?yōu)榈湫投恐笜?biāo)，監(jiān)測(cè)并評(píng)估地質(zhì)活動(dòng)和人為擾動(dòng)等綜合因子對(duì)長城可持續(xù)化保護(hù)的影響； 以通過形變危害識(shí)別，指導(dǎo)墻體的保護(hù)修復(fù)措施和支撐長城景觀廊道的整體科學(xué)保護(hù)。

由于張家口明長城及周邊地區(qū)地形、地貌復(fù)雜，實(shí)地觀測(cè)較為困難； 而遙感因宏觀、客觀和遠(yuǎn)距離探測(cè)具備獨(dú)一無二優(yōu)勢(shì)。受制于時(shí)空失相干、大氣延遲等因素影響，常規(guī)差分雷達(dá)干涉技術(shù)(differential interferometric synthetic aperture Radar，D-InSAR)在自然場(chǎng)景的長城文化景觀時(shí)序形變監(jiān)測(cè)并不適用[2-3]。近年來，為了改進(jìn)D-InSAR技術(shù)的缺陷，獲取時(shí)序形變信息，科研工作者們提出了小基線集(small baseline subsets InSAR，SBAS-InSAR)方法[4-5]。該方法可在抑制時(shí)空去相干的同時(shí)，利用長時(shí)間序列影像獲取自然場(chǎng)景區(qū)雷達(dá)視線向形變場(chǎng)[6-9]。為了彌補(bǔ)長城大型線性遺產(chǎn)系統(tǒng)性形變監(jiān)測(cè)的方法與實(shí)踐空白，本研究基于已有研究成果，利用Sentinel-1升降軌數(shù)據(jù)開展SBAS-InSAR形變反演，并經(jīng)投影變換獲取升降軌垂直向形變速率場(chǎng)。研究結(jié)果可有效探測(cè)并甄別長城顯著形變熱點(diǎn)地區(qū)，為張家口明長城宏觀監(jiān)測(cè)保護(hù)提供科學(xué)數(shù)據(jù)和技術(shù)支撐。

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

1.1 研究區(qū)概況

實(shí)驗(yàn)區(qū)位于張家口市區(qū)以北，崇禮城區(qū)以東，其地理坐標(biāo)范圍為N40.75°～41.06°，E115.17°～115.60°，整個(gè)實(shí)驗(yàn)區(qū)還包括礦坑和不穩(wěn)定坡體以及多個(gè)為2022年張家口冬奧會(huì)修建的場(chǎng)館及配套設(shè)施。實(shí)驗(yàn)區(qū)影像覆蓋范圍如圖1所示，黃線框?yàn)閷?shí)驗(yàn)區(qū)范圍，紅線框?yàn)镾entinel-1升軌影像范圍，藍(lán)線框內(nèi)為降軌影像范圍。

圖1 實(shí)驗(yàn)區(qū)范圍示意圖Fig.1 Schematic diagram of test area

1.2 數(shù)據(jù)源

本次實(shí)驗(yàn)使用歐空局提供的Sentinel-1 SAR數(shù)據(jù)，時(shí)間跨度從2017年5月—2018年7月，包括33景升軌和34景降軌數(shù)據(jù)，成像模式為干涉寬幅模式(IW)，數(shù)據(jù)參數(shù)如表1所示。DEM數(shù)據(jù)采用美國宇航局(NASA)提供的SRTM1數(shù)據(jù)，空間分辨率為30 m×30 m，用于針對(duì)TOPS (terrain observation with progressive)模式的增強(qiáng)譜分集(enhanced spectral diversity，ESD)配準(zhǔn)、模擬地形相位及地形相位去除[10]。

表1 Sentinel-1影像數(shù)據(jù)參數(shù)Tab.1 Parameters of Sentinel-1 data

2 研究方法與數(shù)據(jù)處理

2.1 研究方法

SBAS-InSAR技術(shù)原理為： 根據(jù)時(shí)空基線閾值，對(duì)配準(zhǔn)后的N幅影像進(jìn)行干涉組合生成由M幅干涉圖組成得若干子集。在去除地形相位后，由tB時(shí)刻的主影像與tA時(shí)刻的從影像生成的第i幅干涉圖，在坐標(biāo)點(diǎn)(x，r)點(diǎn)上的干涉相位為[2，11]：

式中：φ表示干涉相位；i表示干涉圖景號(hào)，i=1，2，…，M；λ表示中心波長；d(tB,x,r)和d(tA,x,r)分別表示坐標(biāo)點(diǎn)(x，r)在tB和tA時(shí)刻相對(duì)于起始時(shí)刻的雷達(dá)視線方向(line of sight,LOS)形變累積量；φdem表示外部數(shù)字高程模型(digital elevation model,DEM)帶來的殘余高程誤差相位；φatm(tB,x,r)和φatm(tA,x,r)表示坐標(biāo)點(diǎn)(x，r)在tB和tA時(shí)刻的大氣變化引起的大氣延遲相位；Δni表示噪聲相位。

去除誤差相位后，式(1)可寫成以下形式：

Aφ=δφ，

(2)

式中A為一個(gè)M×N的矩陣。

當(dāng)M≥N，且A為滿秩矩陣時(shí)，可利用最小二乘準(zhǔn)則求解式(2)得到：

φ=(ATA)-1ATδφ。

(3)

若A為秩虧矩陣，式(3)中方程無唯一解，可利用奇異值分解得到A的廣義矩陣，從而得到φ的最小范數(shù)解。

SAR衛(wèi)星沿軌道飛行并進(jìn)行觀測(cè)，通過聯(lián)合地面點(diǎn)的相位值與SAR幅度值，可經(jīng)SBAS-InSAR形變反演得到LOS向的地表形變量Dlos[12]。根據(jù)雷達(dá)成像幾何可知，LOS形變矢量可根據(jù)投影關(guān)系分解為垂直向、東西向及南北向形變矢量，用DU，DE，DN分別表示，投影關(guān)系如下式所示[13-14]：

(4)

式中：θ表示衛(wèi)星入射角；α為方位角(正北方向與衛(wèi)星飛行方向的順時(shí)針夾角)。

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

本次實(shí)驗(yàn)基于開源InSAR處理軟件GMTSAR和GIAnT進(jìn)行SBAS時(shí)序處理，并利用ArcGIS軟件做結(jié)果優(yōu)化與專題制圖。主要步驟包括InSAR干涉處理和SBAS時(shí)序反演。

針對(duì)覆蓋研究區(qū)的升(降)軌影像，選擇20171210(20171209)作為升(降)軌數(shù)據(jù)的公共主影像，基于GMTSAR軟件對(duì)影像進(jìn)行包括數(shù)據(jù)預(yù)處理、增強(qiáng)譜分集配準(zhǔn)、生成干涉圖、去地形相位、干涉圖濾波、相位解纏等處理。其中，設(shè)定時(shí)間基線閾值為48 d，空間基線閾值為200 m，升降軌干涉影像數(shù)據(jù)集共產(chǎn)生110個(gè)和115個(gè)干涉對(duì)。多視處理時(shí)采用方位向視數(shù)為1，距離向視數(shù)為5。相位濾波聯(lián)合采用40 m Gauss與7×7窗口Goldtein自適應(yīng)濾波器以抑制噪聲相位，并進(jìn)而提高相位解纏的可靠性。實(shí)驗(yàn)升降軌時(shí)空基線分布如圖2所示。

利用干涉處理獲取的時(shí)序相干系數(shù)圖，通過Matlab工具包處理得到升(降)軌數(shù)據(jù)的平均相干系數(shù)圖，并以此作為高相干點(diǎn)篩選的依據(jù)。設(shè)定平均相干系數(shù)閾值為0.2，生成平均相干掩模圖與解纏圖、相干系數(shù)圖作為GIAnT軟件的同步輸入數(shù)據(jù)，通過對(duì)解纏圖進(jìn)行平均相干掩模以提高自然場(chǎng)景高相干點(diǎn)空間分布密度?？紤]到張家口長城段隸屬山區(qū)，其與地形相關(guān)的大氣效應(yīng)較為嚴(yán)重，傳統(tǒng)的時(shí)空濾波不能滿足該地區(qū)的大氣相位校正需求[16]。因此亟須引入外部大氣建模數(shù)據(jù)，進(jìn)行大氣系統(tǒng)誤差模擬和糾正，以提升干涉圖質(zhì)量?？紤]到GIAnT軟件附帶的ECMWF(European centre for medium-range weather forecasts)ERA-intrim氣象數(shù)據(jù)當(dāng)天每隔6 h更新，升、降軌影像成像時(shí)間與當(dāng)天最近發(fā)布的同化分析數(shù)據(jù)可分別相差1小時(shí)47分和1小時(shí)41分； 相較于GACOS(generic atmospheric correction online service for InSAR)提供的近實(shí)時(shí)(一分鐘更新)氣象數(shù)據(jù)其時(shí)間分辨率較低。此外GIAnT軟件附帶的ERA-intrim數(shù)據(jù)空間分辨率為0.7°，而GACOS數(shù)據(jù)空間分辨率更高，為0.125°，因此在地形復(fù)雜地區(qū)相較于ECMWF模型有更好的適應(yīng)性[17-19]。綜上所述，本研究選用GACOS氣象數(shù)據(jù)來估計(jì)并校正大氣延遲系統(tǒng)誤差。

此外，考慮到大氣校正后可能引入的趨勢(shì)性相位斜坡，研究基于GIAnT反演線性系統(tǒng)來精確估算每個(gè)SAR影像的軌道參數(shù)，并對(duì)干涉圖進(jìn)行趨勢(shì)校正； 進(jìn)而可利用聯(lián)合DEM誤差估計(jì)的SBAS反演算法估算形變時(shí)序信息?；趐ython工具包，利用線性回歸模型擬合LOS向年形變速率圖，經(jīng)地理編碼得到WGS84坐標(biāo)系下的LOS向年形變速率圖。然后根據(jù)投影關(guān)系轉(zhuǎn)換為垂直向年形變速率圖，投影轉(zhuǎn)換公式為：

(5)

式中：vU表示垂直向形變速率矢量；vlos表示LOS向形變速率矢量； A(D)表示升(降)軌。

3 結(jié)果與分析

3.1 升降軌長城景觀廊道地表形變分析

通過上述SBAS時(shí)序處理，獲取升降軌垂直向形變場(chǎng)，根據(jù)長城墻體設(shè)置250 m的緩沖區(qū)并疊加幅度圖以提供地形信息，得到升降軌長城景觀廊道垂直向形變速率場(chǎng)，如圖3所示。

(a)升軌長城廊道垂直形變場(chǎng)

(b)降軌長城廊道垂直形變場(chǎng)圖3 升降軌長城廊道垂直形變場(chǎng)Fig.3 Vertical deformation field of ascending and descending Great Wall corridor

已知該長城段總長度為85.1 km，由圖3可知，大部分長城段年形變速率在-10 mm/a到10 mm/a之間，但在(E115°27′，N40°44′)附近存在較大的沉降區(qū)，鄰近長城景觀廊道升軌InSAR監(jiān)測(cè)沉降最大值為-34.5 mm/a，降軌InSAR監(jiān)測(cè)沉降最大值為-55.2 mm/a，如圖4(a)和(b)所示。結(jié)合SAR幅度圖及DEM數(shù)據(jù)可知，該處位于山脊，推測(cè)可能存在不穩(wěn)定自然坡體，導(dǎo)致該區(qū)段廊道存在顯著形變。同時(shí)在(E115°13′，N40°47′)附近存在采礦工業(yè)區(qū)，其相鄰景觀廊道升軌InSAR監(jiān)測(cè)沉降最大值為-35.8 mm/a，降軌InSAR監(jiān)測(cè)沉降最大值為-64.5 mm/a，如圖4(c)和(d)所示，表明采礦對(duì)鄰近長城景觀廊道地表穩(wěn)定性有一定影響。實(shí)地考察這2個(gè)主要沉降區(qū)，如圖5所示，定性分析了沉降區(qū)段的主導(dǎo)驅(qū)動(dòng)力，即人類采礦活動(dòng)及自然滑坡風(fēng)險(xiǎn)。考慮到冬奧會(huì)場(chǎng)館等建設(shè)活動(dòng)對(duì)鄰近長城遺址的影響，選擇鄰近冬季2項(xiàng)場(chǎng)館中心明長城景觀廊道，監(jiān)測(cè)結(jié)果發(fā)現(xiàn)升軌InSAR沉降最大值為-41.6 mm/a，降軌InSAR沉降最大值為-44.7 mm/a，如圖4(e)和(f)所示，揭示人工建設(shè)活動(dòng)對(duì)長城廊道周邊地表形變的擾動(dòng)和觸發(fā)作用。

3.2 升降軌形變交叉互檢

由于缺少外部地面實(shí)測(cè)數(shù)據(jù)(水準(zhǔn)或GNSS)，研究利用相同觀測(cè)周期的升降軌SBAS-InSAR形變測(cè)量值的交叉互檢來評(píng)價(jià)形變反演精度及可靠性?？紤]實(shí)驗(yàn)區(qū)包括山地和平地地貌，其中山地占主導(dǎo)。為不失典型性與普適性，根據(jù)地形圖分別選取采礦山區(qū)和龍觀鎮(zhèn)平地區(qū)的長城廊道做剖線，開展升降軌SBAS-InSAR沉降速率精度定量評(píng)價(jià)； 對(duì)應(yīng)剖線的相關(guān)統(tǒng)計(jì)信息，如圖6所示。其中，剖線Ⅰ和Ⅱ位于長城景觀廊道的采礦山區(qū)； 而剖線Ⅲ和Ⅳ位于景觀廊道的平地區(qū)。剖線上選取的相干點(diǎn)沉降速率測(cè)量值，如表2所示。

表2 升降軌沉降速率剖線測(cè)量值Tab.2 Measurements of ascending and descending vertical deformation profiles (mm·a-1)

(續(xù)表)

研究以升軌沉降速率為參考，降軌沉降速率為對(duì)照，采用均方根誤差和最大偏離度來評(píng)估兩者測(cè)量數(shù)據(jù)的一致性，即

(6)

式中：y表示降軌沉降形變值；x表示升軌沉降形變值；N表示升(降)軌沉降形變值個(gè)數(shù)，i=1,2,…,N。

結(jié)果表明(表3)，升降軌4條采樣剖線的形變測(cè)量差異性指標(biāo)： 均方根誤差分別為9.3 mm/a，3.4 mm/a，1.9 mm/a，1.4 mm/a(計(jì)算獲得均方根誤差平均值4.0 mm/a)，對(duì)應(yīng)最大偏離度為3.0～18.1 mm/a； 間接驗(yàn)證了SBAS-InSAR形變反演的定量精度。研究同時(shí)發(fā)現(xiàn)，剖線在沉降量大的山區(qū)可能存在少量不可靠監(jiān)測(cè)點(diǎn)(例如剖線I最大偏離度可達(dá)18.1 mm/a), 如圖6(a)所示； 揭示了在加密山區(qū)場(chǎng)景(地形與地貌相對(duì)復(fù)雜)相干點(diǎn)測(cè)量空間密度的同時(shí)(選用0.2平均空間相干系數(shù)閾值)，不可避免引入了少量隨機(jī)觀測(cè)誤差。

表3 升降軌沉降速率剖線交叉對(duì)比測(cè)度Tab.3 Measurement indices in the cross comparison of vertical deformation profiles between ascending and descending results (mm·a-1)

由以上研究可知，平地區(qū)升降軌沉降速率吻合度高(剖線III和IV)，其均方根誤差在2 mm/a以內(nèi)，表征SBAS-InSAR技術(shù)在平地區(qū)形變監(jiān)測(cè)可達(dá)mm級(jí)； 而在采礦山區(qū)(剖線I和II)，升降軌沉降速率趨勢(shì)總體一致，但仍可表征為較為顯著的均方根誤差： 剖線Ⅰ均方根誤差為9.3 mm/a，剖線Ⅱ?yàn)?.4 mm/a。進(jìn)一步研究發(fā)現(xiàn)，均方根誤差可隨著沉降速率強(qiáng)度的增大而增大。升降軌沉降速率交叉互檢在山區(qū)產(chǎn)生了較大均方根誤差，分析可由地表稀疏植被擾動(dòng)、存在南北-東西向形變、以及InSAR LOS測(cè)量在山區(qū)(坡向、疊掩和陰影等)固有局限性等綜合原因共同決定。

3.3 長城廊道形變風(fēng)險(xiǎn)制圖

綜合形變顯著性水平和平均均方根誤差測(cè)度(4.0 mm/a), 均值融合升降軌長城廊道垂直形變場(chǎng)并以10.0 mm/a為閾值進(jìn)行地表相對(duì)穩(wěn)定和顯著變化專題分類，得到專題風(fēng)險(xiǎn)圖(圖7)[20]，圖7中A，B和C分別標(biāo)示了采礦區(qū)、不穩(wěn)定坡體和冬奧會(huì)場(chǎng)館位置。對(duì)長城墻體左右兩側(cè)各250 m緩沖帶的地表形變趨勢(shì)進(jìn)行統(tǒng)計(jì)，結(jié)果表明，因地表破碎、不穩(wěn)定坡體、人工開礦以及冬奧會(huì)場(chǎng)館修建等綜合因素影響，景觀廊道地表形變速率絕對(duì)值大于10 mm/a閾值的長城段，其長度約為17.5 km，即占比觀測(cè)段總長的20.5%。對(duì)照而言，79.5%占比的長城段沉降速率較小，景觀廊道地表相對(duì)穩(wěn)定。形變風(fēng)險(xiǎn)制圖為后續(xù)明長城景觀廊道及其墻體遺產(chǎn)潛在病害重點(diǎn)勘查提供了靶區(qū)，便利長城大型線性遺產(chǎn)的整體性規(guī)劃與保護(hù)。

圖7 長城廊道地表形變穩(wěn)定性專題分類Fig.7 Thematic deformation riskingmapping of the Great Wall corridor

4 結(jié)論

本文利用Sentinel-1升降軌影像，基于SBAS-InSAR技術(shù)開展張家口明長城景觀廊道2017年5月—2018年7月地表形變前沿示范研究?；贕ACOS大氣相位系統(tǒng)校正、聯(lián)合Gauss與Goldstein濾波的SBAS-InSAR方法獲得了張家口85.1 km明長城景觀廊道統(tǒng)計(jì)意義上mm級(jí)的年形變速率場(chǎng)。研究結(jié)果表明，對(duì)照79.5%的相對(duì)穩(wěn)定段，20.5%的明長城景觀廊道存在較為顯著形變(年形變速率大于10 mm/a)； 為后期長城建筑遺產(chǎn)形變危害識(shí)別、靶區(qū)定位和整改維修等保護(hù)措施的規(guī)劃與落實(shí)提供了定量監(jiān)測(cè)數(shù)據(jù)和全新監(jiān)測(cè)手段。研究表明，基于相干目標(biāo)的SBAS-InSAR技術(shù)在自然場(chǎng)景地區(qū)具備較好適用性； 技術(shù)可望推廣至長城、運(yùn)河等大型線性遺產(chǎn)景觀廊道的整體宏觀監(jiān)測(cè)與動(dòng)態(tài)評(píng)估。隨機(jī)監(jiān)測(cè)誤差的自動(dòng)識(shí)別與噪聲去除是今后算法改進(jìn)方向； 同時(shí)考慮到張家口明長城及周邊地表可能存在南北-東西向形變，聯(lián)合升降軌InSAR數(shù)據(jù)的三維形變反演將是未來工作的又一重要方向。