张家口明长城景观廊道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数据的三维形变反演将是未来工作的又一重要方向。