張家口明長城景觀廊道Sentinel-1影像SBAS形變監測示范研究

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

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

0 引言

張家口明長城墻體遺產裸露于地表并可受自然與人類過程諸多因素影響[1]。盡管長城沿線地區在一定周期和范圍內未發生地震等大型自然災害，長時間緩慢地表形變仍可成為影響張家口明長城文化景觀整體穩定性的重要因素之一。考慮到張家口明長城部分墻段緊鄰采礦工業區和自然不穩定坡體，加之冬奧會修建場館的影響(例如，位于崇禮區東部的北歐中心跳臺滑雪場和越野滑雪場選址規劃緊鄰明長城遺址)； 因此亟須以微形變為典型定量指標，監測并評估地質活動和人為擾動等綜合因子對長城可持續化保護的影響； 以通過形變危害識別，指導墻體的保護修復措施和支撐長城景觀廊道的整體科學保護。

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

1 研究區概況與數據源

1.1 研究區概況

實驗區位于張家口市區以北，崇禮城區以東，其地理坐標范圍為N40.75°～41.06°，E115.17°～115.60°，整個實驗區還包括礦坑和不穩定坡體以及多個為2022年張家口冬奧會修建的場館及配套設施。實驗區影像覆蓋范圍如圖1所示，黃線框為實驗區范圍，紅線框為Sentinel-1升軌影像范圍，藍線框內為降軌影像范圍。

圖1 實驗區范圍示意圖Fig.1 Schematic diagram of test area

1.2 數據源

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

表1 Sentinel-1影像數據參數Tab.1 Parameters of Sentinel-1 data

2 研究方法與數據處理

2.1 研究方法

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

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

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

Aφ=δφ，

(2)

式中A為一個M×N的矩陣。

當M≥N，且A為滿秩矩陣時，可利用最小二乘準則求解式(2)得到：

φ=(ATA)-1ATδφ。

(3)

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

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

(4)

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

2.2 數據處理

本次實驗基于開源InSAR處理軟件GMTSAR和GIAnT進行SBAS時序處理，并利用ArcGIS軟件做結果優化與專題制圖。主要步驟包括InSAR干涉處理和SBAS時序反演。

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

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

此外，考慮到大氣校正后可能引入的趨勢性相位斜坡，研究基于GIAnT反演線性系統來精確估算每個SAR影像的軌道參數，并對干涉圖進行趨勢校正； 進而可利用聯合DEM誤差估計的SBAS反演算法估算形變時序信息。基于python工具包，利用線性回歸模型擬合LOS向年形變速率圖，經地理編碼得到WGS84坐標系下的LOS向年形變速率圖。然后根據投影關系轉換為垂直向年形變速率圖，投影轉換公式為：

(5)

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

3 結果與分析

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

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

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

(b)降軌長城廊道垂直形變場圖3 升降軌長城廊道垂直形變場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′)附近存在較大的沉降區，鄰近長城景觀廊道升軌InSAR監測沉降最大值為-34.5 mm/a，降軌InSAR監測沉降最大值為-55.2 mm/a，如圖4(a)和(b)所示。結合SAR幅度圖及DEM數據可知，該處位于山脊，推測可能存在不穩定自然坡體，導致該區段廊道存在顯著形變。同時在(E115°13′，N40°47′)附近存在采礦工業區，其相鄰景觀廊道升軌InSAR監測沉降最大值為-35.8 mm/a，降軌InSAR監測沉降最大值為-64.5 mm/a，如圖4(c)和(d)所示，表明采礦對鄰近長城景觀廊道地表穩定性有一定影響。實地考察這2個主要沉降區，如圖5所示，定性分析了沉降區段的主導驅動力，即人類采礦活動及自然滑坡風險。考慮到冬奧會場館等建設活動對鄰近長城遺址的影響，選擇鄰近冬季2項場館中心明長城景觀廊道，監測結果發現升軌InSAR沉降最大值為-41.6 mm/a，降軌InSAR沉降最大值為-44.7 mm/a，如圖4(e)和(f)所示，揭示人工建設活動對長城廊道周邊地表形變的擾動和觸發作用。

3.2 升降軌形變交叉互檢

由于缺少外部地面實測數據(水準或GNSS)，研究利用相同觀測周期的升降軌SBAS-InSAR形變測量值的交叉互檢來評價形變反演精度及可靠性。考慮實驗區包括山地和平地地貌，其中山地占主導。為不失典型性與普適性，根據地形圖分別選取采礦山區和龍觀鎮平地區的長城廊道做剖線，開展升降軌SBAS-InSAR沉降速率精度定量評價； 對應剖線的相關統計信息，如圖6所示。其中，剖線Ⅰ和Ⅱ位于長城景觀廊道的采礦山區； 而剖線Ⅲ和Ⅳ位于景觀廊道的平地區。剖線上選取的相干點沉降速率測量值，如表2所示。

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

(續表)

研究以升軌沉降速率為參考，降軌沉降速率為對照，采用均方根誤差和最大偏離度來評估兩者測量數據的一致性，即

(6)

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

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

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

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

3.3 長城廊道形變風險制圖

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

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

4 結論

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