联合InSAR与PCA的北京平原地面沉降时空分析

何 旭, 何 毅*，张立峰，陈 毅，蒲虹宇，陈宝山

1. 兰州交通大学测绘与地理信息学院，甘肃 兰州 730070 2. 地理国情监测技术应用国家地方联合工程研究中心，甘肃 兰州 730070 3. 甘肃省地理国情监测工程实验室，甘肃 兰州 730070

引 言

地面沉降形成时间缓慢、 形成原因复杂、 造成后果严重、 治理难度巨大，在世界各国正以不同严重程度蔓延[1]。 地面沉降监测已在全世界重点区域开展数年。 传统的监测技术获取沉降点的空间密度低、 成本高、 周期长，难以获取大范围、 高密度的沉降数据，给整个区域沉降的时空演变特征分析带来不便，为地面沉降的有效治理带来一定困难。 近年来已有众多学者将合成孔径雷达干涉(interferometric synthetic aperture radar，InSAR)技术，如： PS-InSAR(persistent scatterer)技术与SBAS-InSAR(small baseline subset)技术应用于大范围、 高精度、 连续地面沉降监测[2]。 其中SBAS-InSAR技术避免了因时空基线过长引起的失相干，已在城市、 港口、 冻土区等地表沉降监测领域得到了广泛的应用。

北京作为中国的首都，其平原区域承载着高大密集的建筑群体、 立体发达的交通网络、 孕育着高度集中的人口，不均匀的地表沉降有可能会夺走其孕育的部分生命、 摧毁其承载的大量建筑，带来不可预估的损失与伤害。 因此，对北京平原地表沉降进行实时监测分析，将一切可能存在的隐患提前遏制，是一项必要且长久的举措。 已有众多学者利用不同的监测技术对北京平原地表沉降进行监测，如利用分层标、 基岩标、 水准测量等对北京平原2006年—2019年沉降较为严重的区域进行了监测并分析[3]； 利用PS-InSAR技术监测了2011年—2018年北京平原的地表沉降[4]。 传统监测方法仅研究了平原局部地区的地面沉降，而InSAR技术获取的沉降时间序列大部分研究仅持续到2018年，对于2018年之后的研究较少，并且已有研究仅从时间或者空间的角度对北京平原地表沉降的演化特征进行分析，结合时间和空间上相关性的考虑较少，使得数据中潜在的时空特征以及可能存在的某些规律未被充分挖掘。

主成分分析方法(principal component analysis，PCA)在没有先验约束的情况下，将一组相关变量经过数学变换为一组互不相关的变量，以此表示时间和空间变化的变形模式。 PCA在分析夏季降水、 北太平洋气候的时空特征时，具有适用性、 可靠性。 在地表沉降监测分析中，PCA曾被用于分析加利福尼亚州某山谷处地面沉降的时空演变特征，并成功分离出嵌入在沉降时间序列中随时间变形的时空模式[5]。 目前一些学者使用SBAS-InSAR技术获取火山喷发前的地表形变数据，利用PCA分析其时空演变特征[6]，挖掘了深层时空规律。

本研究使用Sentinel-1A数据，以SBAS-InSAR技术为支撑，对北京平原2017年5月—2020年5月地表形变进行监测，获取地表形变速率以及时序累积沉降量，并进行内部精度验证； 利用PCA对2017年—2020年北京平原地表InSAR时序累积沉降量进行时空特征分析，以更深入地挖掘北京平原地面沉降的时空演变过程，研究结果在北京地面沉降控制管理方面具有一定的科学价值。

1 实验部分

1.1 研究区

北京平原(39°27′—40°29′N，115°39′—117°19′E)位于北京市的东南部(图1)，西面、 北面均为群山环绕，属典型的大陆季风性气候。 平原面积为6 930 km2，西北较高，东南较低，主要由北运河、 永定河、 潮白河、 大清河、 蓟运河冲洪积扇组成[7]。 冲洪积扇顶部，为单一砂卵砾石构成的含水层，埋深较浅，渗透性较好，有利于大气降水补给； 中部地区是由二到三层砂卵砾石层、 砂卵砾石与黏土层、 粉质黏土层互层组成的承压含水层结构，富水条件较好； 下部及平原区则变成由粗砂、 中砂、 细砂、 粉细砂共同组成的多层承压含水层，即： 潜水层(25 m)、 第一承压含水层(80～100 m)、 第二承压含水层(100～180 m)、 第三承压含水层(200～300 m)，蕴含着丰富的地下水资源[8]。 由于沉积环境较为复杂，导致第四系含水层岩性分带也比较复杂。 由平原至山前，第四系沉积物厚度逐渐减小，含水层颗粒由细变粗，地下水埋深由浅变深。 平原年内降水分布不均，大约60%～80%的降水集中在6月—9月，而农作物生长期为每年的4月—6月，此时间段内，农业用水需求量增加，且降雨量不足、 蒸发量较大，使得4月—6月北京地区地下水开采量较大[9]。

1.2 数据源

由阿拉斯加合成孔径雷达设施(Alaskan SAR Facility，ASF)官网获取了2017年5月—2020年5月覆盖北京的39景升轨Sentinel-1A单视复数影像。 其中，成像方式采用IW(interferometric wide swath)，幅宽240 km，采用C波段，波长为5.6 cm，空间分辨率为5 m×20 m，入射角为38.85°，VV同极化方式； 使用由美国航空航天局(National Aeronautics and Space Administration，NASA)所获取30 m分辨率的SRTM-DEM(shuttle radar topography mission-digital elevation model)数据作为外部参考DEM，以估计和去除地形相位； 并从中国气象数据网(http: //data.cma.cn/)获取了北京54511号站点的月降水数据，用于后文分析。

1.3 研究方法

1.3.1 SBAS-InSAR技术

SBAS-InSAR技术最先是由Berardino等于2002年提出，采用多主影像的时间序列InSAR技术进行地面沉降监测[10]。 在植被覆盖区以及地表覆盖变化地区，较短的干涉对时空基线，能够有效抑制失相干以及减少DEM误差等带来的不利影响，使结果更加可靠。

本研究将从ASF网站获取的Sentinel-1A数据，利用SARscape5.2软件进行处理，处理流程如图2所示： 设置时间基线阈值为120 d，空间基线阈值为最大临界基线的45%，生成小基线对(图3)，超级主影像的日期为2019年9月19日； 根据生成的小基线对对影像进行逐像素配准，配准精度要求达到1/8个像素； 为提高信噪比以提供更可靠的相干性以及提高计算效率，将配准后影像的方位向与距离向的多视视数设置为4∶1，计算相干性并进行复共轭相乘得到干涉图，生成的干涉图包含了地形相位、 形变相位、 大气相位、 轨道相位、 噪声相位； 采用参考DEM去除地形相位； 使用Goldstein自适应滤波算法对干涉相位进行滤波，去除噪声相位； 设置解缠相干性阈值为0.2，相干性阈值大于0.2的相位点组建Delaunay三角网，应用最小费用流算法进行相位解缠； 选取质量较好的干涉图，编辑连接对剩余37景影像参与后续反演； 选取相干性良好、 远离形变区域的GCP(ground control points)点进行轨道精炼与重去平，以估算和去除解缠后依然存在的恒定相位和相位跃变； 通过时间高通滤波、 空间低通滤波去除大气相位； 采用最小二乘法和奇异值分解(singular value decomposition，SVD)法获取沿视线向的形变； 采用DEM为参考进行地理编码，得到研究区沉降速率以及时序累积沉降量。

图2 流程图Fig.2 flow chart

图3 时空基线图(a): 空间基线； (b): 时间基线Fig.3 Time-space baseline map(a): Spatial baseline map; (b): Time baseline map

1.3.2 PCA方法

PCA，也称为经验正交函数分析(empirical orthogonal function，EOF)，通过分析数据中的结构特征，以提取数据中主要的特征信号。 在本研究中，假设获取了m个高相干点，n为影像数据，构成矩阵Xm×n，具体算法流程如下：

(1)首先对矩阵Xm×n同一时间观测数据进行去均值处理，见式(1)

(1)

(2)

(3)计算协方差矩阵Cm×m的特征值(λ1, …,λm)和特征向量Vm×m，见式(3)

(3)

将得到的特征值按照从大到小的顺序排列，对应的特征向量也按照相同的顺序进行排列； 特征向量为主成分在空间上的响应，即空间特征。

(4)求解第k个主成分的方差贡献率，见式(4)

(4)

(5)根据方差贡献率求取前j个主成分，见式(5)

(5)

所得到的PCj×n中每行数据就是主成分在时间上的响应，即时间特征。

2 结果与讨论

2.1 时序InSAR监测结果与精度验证

2.1.1 2017年—2020年北京平原地面沉降时空分布

采用SBAS-InSAR技术获取了北京平原沿视线方向(line of sight，LOS)的形变量，地面沉降通常呈垂直向，因此将获取的LOS形变量转换到垂直于地面方向形变量，其垂直向地面平均形变速率如图4所示。 图4(a)表示北京平原2017年5月—2020年5月的地面平均形变速率，可以看出在研究时段内东城、 西城、 石景山、 丰台等区域的地面沉降速率较小，基本呈稳定状态； 昌平东南部、 海淀西北部、 朝阳东部、 通州西北部、 顺义西南部、 大兴南部等区域均存在地面沉降，地面平均形变速率在77.6～-114.9 mm·yr-1之间，其中海淀西北部、 朝阳东部、 通州西北部为地面沉降重灾区，而朝阳东部金盏区域地面沉降速率达到114.9 mm·yr-1。

选取2017年5月—2020年5月Sentinel-1A数据，2017年与2020年数据均不为完整年份，统计年平均地面形变速率会存在误差，如发生季节性沉降等情况，本工作仅统计2018与2019年年均地面形变速率[图4(b，c)]。 图4(b)可以看出，2018年最大沉降速率为126.4 mm·yr-1； 2019年北京平原地面最大沉降速率(109.4 mm·yr-1)相比2018年有所减缓，并且地面沉降区域空间分布也相对减少，2018年朝阳金盏地面最大沉降速率比2019年(105.2 mm·ys-1)大，与文献[3]中水准测量的结果较为一致。 选取了三个不同位置对比2018年与2019年地面形变速率时空特征，如图4(a)中粉色虚线框所示区域。 图4中红色虚线框为三个位置的地面形变速率放大图。 从图中可以看出，位置1与位置3处2019年相比于2018年地面沉降速率超过50 mm·yr-1的区域显著增加； 在位置2处2019年沉降速率超过50 mm·yr-1的范围相比于2018年有较为明显的减小。

图4 北京平原垂直形变速率图(a): 2017年5月—2020年5月； (b): 2018年; (c): 2019年Fig.4 Vertical deformation rate map of Beijing plain(a): 2017.5—2020.5； (b): 2018; (c): 2019

为了详细地揭示2017年—2020年北京平原地面形变时空演变特征，绘制了该区域垂直方向上地面累积形变量的时间序列图(图5)。 从图5中可以看出，2017年5月20日开始至2017年12月22日朝阳金盏地面沉降量达到了91.00 mm，截止到2020年5月4日，朝阳金盏地面最大累积沉降量达到345.90 mm，为北京平原区的最大沉降量所在地。 从时空演变看出(图5)，北京平原地面累积沉降量随时间推移不断增大，影响范围也在持续扩大。 已有学者对其成因进行了论述，但对其时空变化诱因的特征分离，还少有文献对其详细阐述。 因此，将采用PCA方法对北京平原地面沉降时空演变特征的成因进行详细讨论。

图5 北京平原累积形变量图Fig.5 Cumulative deformation map of Beijing plain

2.1.2 精度验证

(1) 内部精度验证

为确保北京平原地面时序形变结果的准确性，根据SBAS-InSAR技术采用均值相干系数判定稳定点为依据，对北京平原地面时序沉降结果进行内部检验。 首先，从筛选得出的时序干涉对中选取不同时段时间间隔接近的时序相干性图进行对比，如图6所示。 图中越亮的区域相干系数值越大，相干性越好，干涉相位结果越可靠，其中，相干性系数可根据式(6)解出

(6)

式(6)中：M和N分别为计算相干性的数据块大小；k和l分别为数据块的行列号；μ1和μ2是干涉对数据中，坐标为(k，l)处的复数值，⊗表示复共轭相乘。

图6 不同时段的相干性图及相干系数Fig.6 Coherence diagrams and coherence coefficients in different time periods

从图6中可以看出几乎相同时间间隔，不同季节的相干性及平均相干系数(见各图右下角)存在较大差异。 相干性由高到底的季节依次为冬季、 初春、 晚秋、 初夏，相干性呈阶段变化的这种情况被认为是植被覆盖变化所产生的结果[11]。 在夏季，植被生长茂盛，使得相干性受到严重的影响。 尽管受到夏季植被生长的影响，研究区内的平均相干系数也保持在0.5及以上，确保了研究结果的准确性[12]。

(2)外部精度验证

杨艳等曾得出朝阳金盏、 黑庄户两地2018年和2019年水准监测的地面最大沉降速率[3]，利用这一数据对结果的精度进行验证(图7)。 图7中实线表示朝阳金盏，虚线表示黑庄户； 蓝色表示InSAR监测结果，红色表示水准监测结果。 朝阳金盏2018年和2019年InSAR监测结果与水准监测结果差值分别为7.62和8.34 mm·yr-1； 黑庄户2018年和2019年InSAR监测结果与水准监测结果差值分别为4.47和-16.76mm·yr-1。 通过对比发现，除2019年黑庄户InSAR监测值与水准监测结果相差较大，其余监测结果的绝对误差较小，与水准监测值具有较好的吻合度。

图7 InSAR监测结果与水准测量结果对比图Fig.7 Comparison chart of InSAR monitoring resultsand leveling measurement results

2.2 基于PCA的北京平原地面沉降时空演变特征成因分析

采用SBAS-InSAR技术获取了2017年—2020年北京平原5 245 675个监测点的时间序列数据，全部参与运算数据量太大，因此以等间隔抽取5 677个监测点数据进行PCA计算。 由于第一期InSAR监测点数据所有值均为0，因此去除第一期InSAR监测值，研究时段共有36期沉降时间序列，构成矩阵X5 677×36，对矩阵X5 677×36进行PCA分析，保留的主成分数量反映了北京平原地表形变不同的时空特征，其空间特征变形模式采用克里金插值法进行插值，深蓝色表示与时间特征呈极大正相关，是重点分析部分。 基于PCA分析发现，前三个主成分便能保留数据集中99.11%的特征[图8(a)]，因此，仅对前三个主成分进行分析。

主成分分析的第一主成分(PC1)解释了数据集96.48%的特征。 时间特征[图8(b)]显示了2017年6月—2020年5月的长期线性沉降过程； 其空间特征[图8(c)]与时间特征呈正相关的区域与沉降速率[图4(a)]以及累积沉降量(图5)的空间位置基本一致，证明PC1解释了2017年—2020年北京平原地面沉降的主要成因。 周朝栋等研究结果表明，北京平原地下水长期过度开采，使得地下水埋深下降明显，有效应力增加，含水层系统固结，造成不均匀的、 长期的地面沉降[13]。 由于PC1的长期线性特征与地下水长期过度开采而引起的地面长期沉降具有相同特征，因此认为PC1揭示了因地下水过度开采引起的沉降。

图4中2018年沉降速率大于2019年，此处以2018年12月最后一景影像为临界，对时间特征进行线性拟合，其中2018年1月—12月的时间特征斜率为-2.6，如图中浅蓝色虚线所示； 2019年1月—12月，时间特征的斜率为-1.9，如图中绿色虚线所示。 2019年时间特征斜率的绝对值比2018年小，这表明2019年整个平原大部分区域的沉降速率在减缓，与InSAR测量结果一致(图4)。 由北京市统计年鉴得知，从2018年至2019年地下水供水量由16.3亿立方米减少到15.5亿立方米[14]； 并且，自南水开始进京后，北京平原地下水开采量逐年减少，平原中部地区(位置2)含水层得到补给，地下水埋深存在小幅度上升，使得土壤孔隙水压力增加，沉降速率减缓[15]。

图8 方差贡献率及第一主成分的时空特征图(a): 方差贡献率； (b): 第一主成分时间特征； (c): 第一主成分空间特征Fig.8 Variance contribution rate and spatio-temporal feature map of the first principal component(a): Variance contribution rate; (b): The first principal component teme feature;(c): The first principal component spatial feature

第二主成分(PC2)解释了数据集2.11%的特征。 其时间特征[图9(a)]分为2017年6月—2018年2月、 2018年9月—2019年2月两个相对平稳期，2018年3月—2018年8月、 2019年2月—2019年7月两个快速沉降期，斜率分别为-8.5和-9.2，两个快速沉降期斜率的绝对值后者大于前者，表明此特征反映区域整体沉降速率在增大； 空间特征[图9(b)]响应较强区域与图4中位置1、 位置3分布基本相同，通过已有文献所统计的北京平原第四系沉积厚度图[图9(c)]，发现PC2空间特征较强的区域与北京平原可压缩层厚度分布较为相似[16]，因此认为其空间特征与可压缩层厚度有关。 研究中查看了位置1与位置3的光学影像[图9(d，e)]，发现空间特征较强的区域广泛分布耕地，每年4月—6月为农作物生长期，农业用水需求量增加，且降雨量不足、 蒸发量较大，使得4月—6月耕地区域地下水开采量增加，在6月之后，随着农业用水的减采、 停采，地下水开采量减少，沉降速率减小。 然而，时间特征的平稳期开始于8月或9月，相比6月之后地下水的减采，存在2～3个月的时间滞后，分析认为农业用水对地下水的索取大部分来自深层承压含水层(第二与第三承压含水层)，埋深为100～300 m，由下而上响应至地面，存在一定的时间滞后。

第三主成分(PC3)解释了数据集0.52%的特征，从时间特征[图10(a)]可以看出明显的周期性变化，从中国气象数据网获取研究期间的月降雨数据，如图10(a)所示，2017年11月降雨停止，在2018年2月时间特征出现拐点，斜率由正变负，此时为地面沉降期； 2018年3月开始出现降雨，时间特征在同年6月出现拐点，斜率由负变正，此时为地面隆起期，但随着5月降雨的骤减，时间特征的斜率在2018年7月变为负值，又使得地面进入短暂的沉降期； 自2018年7月开始，平原地区大量降雨，时间特征信号在同年9月达到峰值，此时斜率由负变正，地面进入隆起阶段； 但11月降雨停止后，时间特征在2019年2月达到峰值，斜率变为负值，代表地面又进入沉降阶段，随后开始循环此过程。 从这个过程中可以看出降雨开始的2～4个月后，地面入渗补给地下水，使土壤孔隙水压力增加，地表轻微隆起，PC3时间特征的斜率为正； 降雨停止2～4个月后，由降雨提供的地面入渗补给也相应停止，土壤孔隙水压力减小，有效应力增加，地表发生沉降，PC3的时间特征为负。 其空间特征[图10(b)]响应出现在平原西部山前冲洪积扇以及平原中部。 平原西部山前冲洪积扇地质岩性主要为砂卵砾石，入渗补给能力强。 降雨对山前冲洪积扇补给，使得此处空间特征较强。 在平原中部，空间特征较强区域分布与地面沉降严重区域较为相似，可能是因为此处地下水常年超采形成了地下水降落漏斗，降雨入渗后，此处地下水位埋深较低，受侧向补给的作用出现了

图9 第二主成分时空特征以及可压缩层厚度图(c图引自文献[16])(a): 第二主成分时间特征； (b): 第二主成分空间特征； (c): 可压缩层厚度； (d): 位置1处光学影像； (e): 位置3处光学影像Fig.9 Spatiotemporal characteristics of the second principal component and the thickness of the compressible layer(a): The second principal component time feature; (b): The second principal component spatial feature;(c): The thickness of the compressible layer; (d): Optical image map at position 1; (e): Optical image map at position 3

较强的空间特征。

引起地面沉降的原因众多，如地下资源超采、 城市高大密集的建筑物荷载、 不同地质结构、 土地利用类型等，这些因素可能独立驱动地面沉降，也可能存在一定的相关性，共同作用于地面沉降。 在共同作用于地面时，需要分析各因素对地面沉降影响的重要性，以确保科学治理地面沉降。 用于分离混合在地面沉降时间序列中不同特征的传统方法依赖于具有先验模型的函数，而PCA是在没有先验约束的情况下，将一组相关性的变量，通过正交变化，转换为一组线性不相关的变量，以表征不随时间变化且不相关的空间特征部分以及仅依赖时间变化且不相关的时间特征部分，并根据特征值的大小确定各分量对地面沉降贡献率的大小。 本研究根据北京平原的InSAR时间序列数据，利用PCA进行分析，成功分离出三种主要的时空变形特征，并根据特征分析成因，以达到数据深层规律发掘的目的，为未来大范围地表沉降数据分析提供一定的参考。

3 结 论

利用SBAS-InSAR技术，对Sentienl-1A的时间序列数据进行处理，获取了北京平原2017年5月—2020年5月的地面沉降速率以及累积沉降量，并对时间序列沉降进行主成分分析，主要结论如下：

(1)2017年—2020年北京平原最大沉降速率出现在朝阳金盏地区，为114.9 mm·yr-1； 最大累积沉降量达到345.9 mm。 2019年平原中部沉降严重的区域相比于2018年有所减少； 昌平东南部、 海淀西北部、 大兴南部的沉降速率超过50 mm·yr-1的范围不断扩大。

(2)利用主成分(PCA)分析了北京平原地面沉降时间序列，得到了其时间特征以及相应的空间特征。 第一主成分的时间特征描述了沉降的线性过程，其空间特征与沉降速率空间分布一致，解释了地面沉降时空变化的主要趋势，在2019年前后的斜率变化表明南水进京缓解了北京平原地下水的超采状况； 第二主成分的时间特征描述了以年际变化为周期的沉降过程，其空间特征与可压缩层厚度以及土地利用类型在空间上具有较强的相关性，与2018年相比，空间响应较强区域2019年沉降影响范围在持续增大，证明这些区域的承压地下水仍被过度开采； 根据第三主成分的时空特征判断其主要为由降雨调控以年际为周期的弹性过程，其时间特征滞后于降雨2～4个月。

主成分分析通过对广泛分布于空间的时序监测点进行分析，寻求方差最大投影方向，以提取数据中主要的时空特征信号。 但主成分分析只对原始数据做了二阶统计，得到的特征信号仅仅为不相关，并未使其相互独立，使得分离的时空特征与真实特征可能存在差异，为避免这种情况，未来可以利用高阶统计方法，如独立成分分析对数据进行四阶统计分析。 本研究分离的时空特征是通过方差贡献率解释其重要性，并无实际量纲，未来将研究获取具有时空量纲的特征。