基于双极化时序InSAR技术的地表形变监测

2022-05-26陈佳炜吴希文陈炳杰

广东工业大学学报 2022年3期

陈佳炜，吴希文，王华，陈炳杰

（1. 广东工业大学土木与交通工程学院, 广东广州 510006；

2. 广东工业大学大湾区城市环境安全与绿色发展教育部重点实验室, 广东广州 510006）

近年来，我国城市地面沉降灾害频发[1]，严重影响社会的稳定运行。为了防治和减小地面沉降及其带来的影响，精确而高效的地面沉降监测是必要手段。现阶段，精密水准测量和全球定位导航测量等传统大地测量方法虽能提供较高精度的测量结果，但在大范围的沉降监测中需要耗费巨大的人力物力，存在一定的局限性[2]。

作为空间大地测量和遥感技术结合体的星载合成孔径雷达干涉测量技术(Interferometric Synthetic Aperture Radar, InSAR)能提供高精度、高分辨率、大尺度的形变监测结果，是现今国内外监测地面沉降的重要手段和研究方向[3-5]。其中，永久散射体干涉测量(Persistent Scatterer Interferometry，PSI)技术作为时序InSAR技术的一种，被广泛应用于地表形变监测[6-7]。

时序InSAR技术的原理是通过获取多景SAR影像中的高相干点，并对高相干点的相位信息进行大气延迟、轨道、高程等误差的改正，最终获取地表形变结果[8-11]。然而，由于该技术是以影像中的高相干点为基础进行的，因此在相干性差的区域容易出现由于高相干点密度过低或无高相干点导致的形变监测结果差甚至无形变结果的情况。由此，提高高相干点的密度和质量成为了该技术领域的研究热点。

伴随合成孔径雷达(Synthetic Aperture Radar，SAR)传感器技术的发展和新一代卫星的发射，SAR数据实现了由单极化向双/多极化的转变。双/多极化数据的出现为解决上述技术难题提供了新的方向。1998年，Cloude等[12]提出采用Pauli向量转换而成的简化散射矢量对全极化SAR图像中的每个像元的像素值进行表示。2009年，Pipia等[13]成功利用地基SAR所拍摄的双同向极化数据进行时序InSAR处理，获取了西班牙萨连特村的地面形变结果。2010年，Navarro-Sanchez等[14]在Pipia的研究基础上，成功将星载SAR的双同向极化数据应用于时序InSAR技术中。2014年，Iglesias等[15]提出依据选取像素最优化通道，从而得到像素的最优相位的极化优化方法。2019年，Zhao等[16]基于平均时间相干性/振幅离差指数最优原则，提出了结合相干矩阵分解法(Coherency Matrix Decomposition, CMD）和最优极化通道法（BEST）的CMD-BEST方法。基于上述研究，本文依据散射体的散射性质对Sentinel-1A卫星影像的两种极化模式(VV＆VH)数据进行组合，并基于像素振幅离差指数最优的原则进行数据优化，将优化后的数据进行时序InSAR处理；将双极化数据处理结果与单极化数据处理结果进行对比分析，验证该方法的可行性和可靠性；对实验结果中的重点沉降区域进行分析，确定其沉降规律和成因。

1 研究区域概况及数据源

1.1 研究区域及概况

高栏港经济区位于珠海市金湾区西南沿海，总面积约为380 km2，是依托华南沿海主枢纽港高栏港而设立的经济功能区。区内遍布海运、石油、化工等多种产业，是广东海洋经济中最具活力和潜力的地区之一。

本文所选研究区域的经纬度范围为113.16°E～113.25°E、21.9°N～22°N，属于填海区，面积约为145.75 km2，位于高栏港经济区的中部。研究发现，该地区地质条件极不稳定，除土壤含水量高外，区域中大量淤泥为欠固结软土，存在极高的沉降风险[17-18]。同时该区域内快速的经济发展建设和土地开发更加剧了沉降灾害发生的风险。因此，该研究区域具有较高的监测价值和研究意义。

1.2 数据源

实验选取覆盖研究区域的47景双极化Sentinel-1A卫星影像作为数据源，影像参数如表1所示。数字高程模型(Digital Elevation Model, DEM)选用航天飞机雷达地形测绘使命(Shuttle Radar Topography Mission, SRTM)的数据。

表1 雷达影像参数Table 1 The parameters of radar satellite data

2 原理与方法

PSI技术的技术路线是从一系列SAR影像数据中选出高相干性的像素点，即永久性散射体，通过对其相位组成成分进行时空分析和一系列的误差剔除，最终获得高精度的地表形变结果[8-9]。在PSI技术中，永久散射体的选取通常以振幅离差指数DA作为标准[9]，振幅离差指数DA表示为

式中： σa为单个像素的振幅的标准差，a¯为单个像素振幅的均值，N表示影像的景数，S为像素的复数信号值。为保证所选像素点具有较高的相干性，通常选择0.4作为阈值，即当像素的振幅离差指数DA≤0.4时，将被视为PS候选点。

在全极化SAR影像数据中，每一个像素值可以用简化散射矢量K进行表示[12]。

式中：S表示像素的复数信号值，HH、VV、HV表示各个极化方向。对于双极化的Sentinel-1A数据，简化散射矢量K表示为[19]

在双极化数据中，采用复酉矢量 ω对散射体的散射特性进行表示，其表达式为[14,20]

式中： α表示对应散射体的散射机制， ψ为散射相位角。将每个像素的目标散射矢量K投影到复酉矢量ω上求得其对应的散射系数μ，其计算公式为[21]

式中： *表示共轭运算，散射系数μ 等同于单极化数据中像素的复数信号值。因此，通过将式(3)～(5)代入到式(1)，得到双极化数据模式下振幅离差指数DA的计算公式为

在时序InSAR技术中，振幅离差指数DA反映像素相位信息的可靠程度，像素的振幅离差指数DA越小，该像素的相位信息越可靠[8-9]。基于此理论，本实验在散射参数 α 和 ψ的定义范围内，以0.005 rad为间隔，采用穷举法计算所有散射参数组合所对应的像素振幅离差指数DA，进而求得像素的最小振幅离差指数DA及其对应的散射参数组合，最终求得每个像素的最佳散射系数μ。在完成所有像素的最佳散射系数 μ的计算后，选取振幅离差指数DA≤0.4的像素点作为PS候选点，将其极化信息和相位信息利用GEOS-PSI软件进行时序InSAR处理，筛选出PS点，并利用PS点获取研究区域的地表形变结果，实验流程如图1所示。

图1 数据处理流程Fig.1 Data processing flow

3 实验结果及分析

3.1 高相干点的数量和密度

为验证本实验方法的可行性，本文将对双极化数据和同向极化VV/交叉极化VH数据所处理得到的PS候选点和PS点的数量和密度进行对比，结果如表2所示。

由表2可得，在PS候选点数量和密度方面，采用双极化数据相比于采用单极化(VV/VH)数据提高了6.1倍和7.8倍。在PS点数量和密度方面，采用双极化数据相比于单极化(VV/VH)数据提高了2.1倍和2.7倍。相较于PS候选点，PS点的数量和密度提高的幅度较小，其原因是大量像素点经过相位优化后被选为PS候选点，但其时间相干性仍然较差，无法经过筛选成为PS点，由此造成PS点的密度和数量的提高幅度远低于PS候选点的情况。以上实验结果表明：采用振幅离差指数最优的方法对双极化数据进行处理能有效提高InSAR特征点的数量和密度。

表2 PS候选点及PS点的数量和密度Table 2 The number and density of PSC and PS points

图2为研究地区中部分区域的PS点分布情况。如图所示，区域Ⅰ和区域Ⅱ中PS点的分布密度大小均依次为：双极化＞同向极化VV＞交叉极化VH。此外，相较于单极化(VV/VH)数据的PS点获取结果，采用双极化数据所获得的PS点分布范围更广。综上可得：采用双极化数据能够有效提高PS点的分布密度和扩大PS点的分布范围，从而能够有效地提高时序InSAR技术的监测水平。

图2 PS点分布情况Fig.2 The distribution of PS points

3.2 地面形变结果验证及分析

本文利用双极化时序InSAR方法对实验区域的形变情况进行探究，得到了研究区域的地表形变结果。实验结果显示，研究区域的年平均形变速率范围为-69.4～18.3 mm/a。为对实验结果进行验证，本文将双极化数据与单极化(VV/VH)数据的处理结果中的研究区域整体平均沉降速率及其标准差进行对比(见表3)。依表3可知，3种数据的沉降结果在整体上保持一致。三者的平均沉降速率差均不超过0.2 mm/a，其沉降速率标准差的差值在1.1 mm/a以内。同时，双极化数据的处理结果与同向极化VV数据的处理结果更为接近，其主要原因是实验过程中所采用的极化加权法，对相干性更好的同向极化VV数据赋予了更高的权重。图3为利用双极化数据处理得到的研究区域地表形变结果。从图3可知，研究区域内存在多处沉降；大部分沉降区的平均沉降速率超过了30 mm/a，严重沉降区形成了沉降漏斗。为对双极化数据处理的沉降结果进一步验证，本文将其与赖波等处理2017～2018年的Sentinel-1A卫星影像得到的沉降结果[22]、李治斌等处理2017～2019年的Sentinel-1A卫星影像所得到的沉降结果[23]进行对比，发现双极化数据处理得到的沉降区域与上述两篇文献的沉降区域在分布和范围上基本一致；同时，对应沉降区域的沉降速率十分相近。综上，通过与两种单极化(VV/VH)数据的实验结果及该区域已有监测结果进行对比，验证了本实验中双极化数据处理结果的可靠性。

表3 双极化、单极化VV/VH沉降结果比较Table 3 Comparison of displacement result between dual-pol and single-pol VV/VH

图3 实验区域沉降速率图Fig.3 Displacement of the experimental area

为探究其沉降成因和规律，实验选取严重沉降区中的特征点进行时序分析，位置如图3所示，其经纬度分别为：A (113.223 4°N, 21.985 7°E)、B (113.234 1°N, 21.993 6°E)、C (113.238 3°N, 21.954 4°E)。图4为所选特征点A、B、C的沉降时间序列。由图可得，3个特征点均呈线性沉降的趋势，且沉降速率均超过40 mm/a。其中A点的沉降情况最为严重，其沉降速率超过60 mm/a。经调查，3个严重沉降区在监测时段内并未进行大型的施工建设，因此排除了由于房屋建设导致地面沉降的可能。与此同时，通过对其地质条件进行调查，发现该地区地处海陆交界，故土质多为淤泥类软弱土层，且大部分软弱土层为未固结土，含水量超过60%，具有极高的压缩系数。在这种极不稳定的地质条件下，一旦受到较大外力作用，容易发生地表沉降灾害。此外，调查资料显示该研究区域在2008年到2016年间进行了大规模的填海造陆，并在其施工过程中选用了先分隔围堰，后陆上抛石挤淤进行回填的填海方式，该方法导致填土层需要经过由未固结态到固结态的转变，在此过程中势必出现地表的沉降。综上可得，该地区的地质条件和填海工程是导致该地区发生严重地表沉降的主要原因。鉴于该地区出现了较为严重的沉降现象，需要对该区域进行持续的监测并采取必要的防治措施，避免造成经济和人员的损失。