双极化Sentinel-1数据极化时序InSAR技术地表形变监测—以上海市浦东机场为例

冯瀚， 赵峰， 汪云甲， 闫世勇， 彭锴， 王腾，张念斌， 徐东彪

1. 中国矿业大学 环境与测绘学院, 徐州 221116;

2. 国土环境与灾害监测国家测绘地理信息局重点实验室, 徐州 221116;

3. 黄河勘测规划设计研究院有限公司, 郑州 450003

1 引 言

时 序InSAR 技 术PSI （Persistent scatterer Interferometry）在地表变形监测上具有监测范围广、监测精度高等优势，目前已被广泛用于大区域地表形变探测（郭华东和张露，2019）。PSI 技术最早由Ferretti等（2001）提出，该技术作为差分干涉合成孔径雷达DInSAR（Differential Interferometric Synthetic Aperture Radar）的延伸，为减小差分干涉中的失相干及大气延迟影响，对多幅SAR 影像中的永久散射体PS（Persistent Scatterers）进行时序分析，其适用范围和形变监测精度均优于传统DInSAR 技术（朱建军 等，2019）。之后，研究学者相继推出SBAS（Berardino 等，2002）、干涉点目标分析（Werner 等，2003a 和2003b）、StaMPS（Hooper 等，2004）、时域相干点目标分析（Zhang等，2011）等时序InSAR 技术，并在各类变形监测中进行了应用。

但是，时序InSAR 技术在人工建筑物较少的非城市区域，可获取的PS 监测点密度低，难以准确反映地表变形，对相干性较低的分布式散射体区域的形变探测存在不足（Cao等，2015；Crosetto等，2016；林珲 等，2017）。因此，国内外研究学者基于经典PSI 算法提出了更高级的分布式DSI（Distributed Scatterer Interferometry） 技术，如联合永久散射体（PS）与分布式散射体（DS）的SqueeSAR 技术（Ferretti 等，2011）、SqueeSAR 的各类变体技术（Parizzi 和Brcic，2011；Wang 等，2012）以及CAESAR 技术（Fornaro 等，2015），很好地弥补了传统PSI 对DS 目标监测的不足，推动了PSI的进一步发展。但是上述方法只使用了单极化数据，对多极化SAR 数据信息挖掘研究相对较少。然而，各类地表散射体对不同极化的响应不同，利用不同极化通道的丰富极化信息可提高PSI技术地表形变监测能力（Lee和Pottier，2009）。长期以来，长时序极化数据获取难度大，基于多极化数据进行时序地表形变监测较为困难。但是，随着近年来具备多极化SAR 数据卫星传感器的相继发射，如Radarsat-2、TerraSAR-X、ALOS-2、Sentinel-1 等，已经具备了将PSI 技术扩展到多极化的现实基础。2009 年，Pipia 等（2009）基于全极化地基SAR 数据，最早提出了极化时序InSAR的PolPSI（Polarimetric PSI）并首次成功将极化优化算法应用于PSI技术；此后许多学者利用星载极化SAR 数据开展了大量的研究工作，获得了较好的研究成果（Navarro-Sanchez 等，2010；Navarro-Sanchez 和Lopez-Sanchez，2014；Iglesias 等，2014；Zhao和Mallorqui，2019a和2019b）。

Sentinel-1 数据作为一种可免费获取的多极化SAR 数据，现有相关研究对其极化信息的利用不足，已有的研究案例只是简单对比两种极化数据的PSI 监测效果（熊佳诚 等，2019），或者是基于单一相位质量指标（例如DA或相干系数）对所有像元进行统一的干涉相位优化（Shamshiri 等，2018；Azadnejad 等，2020），并未针对PS 与DS 目标不同特性进行自适应的极化优化，Sentinel-1 极化数据在地表形变中的应用潜力还需进一步挖掘和研究。

有鉴于此，为充分利用Sentinel-1 数据的双极化信息，本文借鉴DSI技术和相关研究思想，针对PS 与DS 目标各自的极化特性，提出了一种基于双极化Sentinel-1 数据的PS 与DS 目标自适应极化优化的PolPSI 方法。以上海浦东国际机场为研究区，对比分析本文方法较传统PSI技术优势，并对相关结果进行讨论。

2 研究区概况及数据介绍

浦东国际机场为4F 级民用机场，是中国3 大门户复合枢纽之一。该机场以南北向姿态坐落于长江三角洲冲积平原，地处上海市浦东新区，距上海市中心约30 km。据资料显示，浦东机场旅客吞吐量居全国第二，仅次于首都机场，目前已拥有4 条运行跑道、两个航站楼以及在建的第五跑道，其中第四跑道于2015 年建成并正式投入使用（刘冬明，2018）。该机场所在区域土层软弱，机场由西向东地质条件差异大：其西边场地成陆时间早，东边场地为围涂驻塘后形成的滩地，地质条件较为脆弱复杂。图1为研究区域示意图：浦东机场位于左图中白框区域；右图为机场的Sentinel-1强度影像时序平均图，图中黄色圆形与红色方形标志分别标示典型PS、DS 目标位置，将用于后续结果分析。

图1 研究区域及SAR强度影像图Fig. 1 Study area and SAR intensity image

本文选取34 景双极化Sentinel-1 数据进行PolPSI时序形变监测。图2为Sentinel-1影像时空基线分布图，监测时间从2016年5月至2017年12月，以2017年5月22日为主影像。

图2 SAR影像时空基线分布图Fig. 2 SAR images’ spatial and temporal baseline distribution

3 方法原理

本文借鉴DS-InSAR 思想，提出了PolPSI 技术。不同于常见DS-InSAR 技术对DS 目标基于其同质像元平均来进行优化，本研究利用最小均方差MMSE（Minimum Mean Square Error）滤波方法（Lee 等，1999）对DS目标进行处理。该方法是一种同质像元加权平均技术，可更好保持地物细节信息。提出PolPSI技术整体方法流程如图3所示。

首先通过统计检验识别同质像元（Jiang 等，2017），并设置同质像元阈值划分DS 与PS 目标（Ferretti 等，2011）。在划分目标过程中，对DS 目标再次使用DA准则（DA设为0.25）进行剔除，剔除的像元划归为PS（图3）。

图3 PolPSI方法流程图Fig. 3 Flowchart of PolPSI

之后，针对不同目标，为保持PS 目标的空间分辨率，使用DA准则对PS目标进行极化优化；对DS目标，为最大程度提升其相位质量，先进行DS 目标MMSE滤波，之后使用相干性准则进行极化优化。

最后，基于生成的优化差分干涉图，通过TPC方法选取高质量像元（Zhao 和Mallorqui，2019c）进行PSI 处理（Hooper 等，2004），得到PolPSI 形变监测结果。本文TPC阈值选为0.92（对应相位标准差约为25°）。

3.1 DS目标MMSE自适应滤波

DS 目标MMSE 自适应滤波主要包括同质像元识别和基于同质像元识别结果的极化相干矩阵MMSE 滤波两个步骤。本研究使用蒋弥（Jiang 等，2017）提出的HTCI 方法，以15×15 的窗口识别同质像元，结果如图4（b）所示：结合Google Earth影像分析（图4（a）），识别结果在机场跑道、绿植草地等地区同质像元个数较多，在建筑物等区域的同质像元个数较少，符合实际情况，识别效果较好。

图4 浦东机场影像Fig. 4 Pudong airport image

MMSE滤波能很好地保留影像的极化特征和边缘清晰度（Lee 等，1999；Lee 和Pottier，2009）。该滤波方法基于识别的同质像元进行自适应加权平均滤波，较常见DSI方法中使用的同质像元平均滤波具有更好的干涉图滤波效果。

式（1）—（3）给出了VV 极化与VH 极化哨兵数据的Pauli 基向量ki与极化相干矩阵T4矩阵的构建方法（Zhao和Mallorqui，2019b）：

式（1）—（3）中，S 代表SLC 影像数据，下标i、j表示两个不同时刻，上标T表示矩阵的转置运算，H表示矩阵的共轭转置。极化相干矩阵T4的滤波权重b按式（5）计算，并通过式（6）完成MMSE 滤波（Lee等，1999）。

式（4）—（6）中，tr表示矩阵的迹，下标SHPS表示同质像元，下标mean 代表窗口内同质像元的平均，coef 值取0.51（Lee 等，1999；Lee 和Pottier，2009）。影像滤波前后假彩色显示如图5 所示，图中用于合成为假彩色影像的R、G、B 三个通道分别 为R =|SVV|2，G =|SVH|2，B = |SVV+SVH|2，其 中第二行对应于第一行白框区域的放大显示，黄色圆形标示PS1目标，红色方形标示DS1目标，后续讨论了两种典型目标的极化优化空间参数搜索结果。从图5中可看到MMSE 滤波较好保持了强散射体的空间分辨率，如城市区域强散射体的旁瓣信号细节在滤波后被保留下来；滤波后同质像元区域斑点噪声明显降低，同质区更加平滑。

图5 假彩色影像Fig. 5 The false color image

需要说明的是，利用MMSE 自适应滤波处理之后的极化相干矩阵Tfilter，得到优化后VV 极化差分干涉图（干涉图中PS 目标保持不变，DS 目标基于同质像元进行自适应滤波处理），进一步利用常见DSI 时序分析技术可得地表形变监测结果。该DSI 形变监测结果（MMSE 方法结果）将用于后续同本文提出PolPSI方法结果对比分析。

3.2 干涉相位的时序极化优化

时序极化优化通过构建极化空间ω，遍历搜索使像元质量最优的空间投影方向，随后将影像重投至该方向，获取最优极化干涉相位。对双极化哨兵影像，空间投影矢量ω构建如下（Zhao 和Mallorqui，2019b）：

振幅离差为SAR 强度信息的统计，相干性可以用于评价干涉质量。通常利用振幅离差来识别永久散射体（Ferretti 等，2001），用相干性评价干涉图相位质量（Lee和Pottier，2009）。因此，搜寻最优极化也就是寻找使振幅离差最小或相干性最高的空间投影方向，并求解出α与ψ参数。极化数据的相干性（Shamshiri 等，2018）公式如式（8）和式（9）所示：

式（8）和（9）中k表示第k幅干涉图，Nintf表示干涉图总数量，γ表示干涉图之间的相干性，Ωij与Tii为式（3）极化相干矩阵中的分量。

极化数据的振幅离差DA（Shamshiri 等，2018）公式如式（10）和式（11）：

式中，N表示第N幅SLC 影像，ki表示式（1）里的Pauli 基向量。需指出PS 目标优化是基于SLC的，即最优投影空间ω是使SLC像元质量最好的矢量方向，后续需要基于优化SLC 重新构建差分干涉图。

图6 展示了PS1 目标（图5 中黄色圆形标志）与典型DS1目标（图5中红色方形标志）的极化空间参数搜索结果。可发现PS与DS目标的极化参数在极化空间中呈平滑变化，并在某一位置取到最优（图中白色星号所示），说明优化方法的有效性。

图6 典型散射体空间参数搜索结果Fig. 6 Parameter searching of typical scatters

4 结果与讨论

为分析本文提出的PolPSI 技术性能，将其与传统单极化PSI 以及基于MMSE 滤波的DSI 方法（只利用VV 极化干涉图，既MMSE 方法）结果进行比较分析。此外，从极化散射特性角度，对本文PolPSI优化结果进行讨论。

4.1 差分干涉图结果分析

选取时间基线最长的第一幅差分干涉图进行质量分析。显示原始VV 极化干涉图、MMSE 滤波后VV 极化干涉图与时序极化优化干涉图如图7 所示。图7中可以看出滤波、极化优化后的干涉图条纹比原始干涉图更为清晰，时序极化优化子图7放大的机场航站楼（黄色虚线圈出）条纹最平滑，轮廓最清晰，说明时序极化优化较好地保持了分辨率、相位优化质量最好；进一步选取相位梯度均值、相位标准偏差、残差点数目均值作为相位优化的评价指标，这3种评价指标越小，对应干涉图噪声水平越低、优化效果越好。如表1 所示，在3 项指标上，MMSE 滤波相比于原始干涉图分别下降了10.76%、49.46%、22.81%，时序极化优化（表中黑色字体加粗显示）相比原始干涉图分别下降了19.97%、57.47%、41.75%，优化效果最好。

表 1 干涉图质量评价结果Table 1 Interferogram quality evaluation results

图7 差分干涉图Fig. 7 Differential Interferogram phase

此外，计算时序平均相位标准偏差如图8 所示，图中黄色圆形与红色方形标志分别标示典型PS、DS 目标位置，后续统计散射体优化参数表格（表2）。图8 中可看出本文方法在机场跑道区域（红色实线表示区域）以及建筑物区域取得了更好的优化效果，其时序平均相位标准偏差最小；同时，统计不同方法时序平均相位标准偏差直方图如图9所示，可明显看出两种优化方法都降低了时序干涉图的相位标准偏差，说明两种方法对时序干涉图都起到了优化效果。其中，本文的极化优化方法效果最好。

图8 时序平均相位标准偏差图Fig. 8 Mean phase standard deviation of Interferograms

图9 不同方法时序平均相位标准偏差直方图Fig. 9 Standard deviation histograms of different methods

4.2 极化散射特性讨论

为比较典型散射体最优极化优化效果，分别统计了3个典型PS目标（图1中红色圆形标示）与DS 目标（图1 中红色方形标示）的相位标准偏差与优化空间参数如表2所示。从表中可以看出极化优化最大程度降低了两种典型散射体的相位标准差。根据Azadnejad 等（2020）的研究，极化空间参数α 值趋于0°对应二面散射体、垂直极散射体等，VV 极化的贡献度较大；α趋于90°对应随机体散射，VH极化能有较好响应；而对于空间参数ψ，趋于±180°时表示目标对垂直极化响应较好，越趋于0°则表示水平极化具有较好的响应。表2中，PS散射体α 值趋于0°，ψ值趋于-180°，DS 则表现相反，两种典型目标满足此规律。

表 2 典型散射体优化参数统计Table 2 Statistics of optimization parameters for typical scatterers

为进一步说明最优极化与地表散射体之间的联 系（Cloude 和Papathanassiou，1998；Navarro-Sanchez 等，2010），计算所有PS 与DS 目标极化优化前后时序相位标准偏差均值，优化前分别为0.51、1.73，优化后为0.33、0.81；同时，统计研究区内所有PS与DS目标最优空间参数，绘制统计直方图如图10所示。

图10 最优参数统计直方图Fig. 10 The optimal parameter histograms

从图10（a）和10（c）所有PS 目标的统计图中可以看出，其α 最优值基本集中在0°附近，ψ最优值集中分布在-180°附近，表明PS 目标对VV 极化的响应更好，因此在使用哨兵数据进行城市区域PS 目标形变监测时，利用VV 单通道可得到较好的监测效果；对于DS 目标而言（图10（b）和10（d）），α 最优值从0°到90°均有覆盖，ψ最优值从-180°到180°近似均匀分布，表明VH 极化对于DS 目标优化有一定的贡献度，从侧面说明VH 通道中蕴含一定的DS相位信息。

4.3 形变监测结果

VV 极化PSI、MMSE 滤波DSI、时序PolPSI 这3 种方法时序平均形变结果如图11 所示。传统PSI算法获得29961 个高质量监测点；基于MMSE 滤波的DSI 方法得到46511 个高质量监测点，相比传统PSI 增加了55%；本文PolPSI 获取到60829 个高质量监测点，相对传统PSI 则增长近103.0%，相对MMSE的DSI方法增加30.8%。图11 中东北角沿海岸边的线性目标为修建的沿海堤坝，在针对该沿海堤坝的变形监测上，传统PSI算法仅能选取到零星的高质量监测点；基于MMSE 滤波的DSI 结果对堤坝上的沉降有所体现；本文PolPSI 方法结果，该线性堤坝上点位有明显增加。

图11 地表形变监测结果（白框区域放大图位于最右下角）Fig. 11 Ground deformation velocity maps obtained （The close-up area of white rectangle are locating in the bottom right corner）

除了沿海堤坝工程，选取第四跑道南北两端的典型区域进行放大展示：传统PSI方法在跑道同质区基本没有监测出沉降信息；而基于MMSE 的DSI技术能监测出一些形变点，但是密度不足，不能很好反映跑道上的形变情况，甚至在跑道南端区域没能探测出形变信息（图11（B2））；本文PolPSI在跑道上的点位数量明显增多，沉降信息更为突显，跑道南部区域（白框C1，C2 区域）探测出了小范围的沉降（图11（C2））；此外，PolPSI方法在其它区域的表现也更好。本文研究区监测结果与前人利用高分辨率SAR 数据研究结果吻合（廖明生 等，2020）。

为进一步说明3种方法形变结果可靠性，统计了不同方法中共同监测点的形变速率差异直方图，如图12 所示。从图12 中可看出，本文PolPSI 方法与MMSE 滤波的DSI 方法和传统PSI 的沉降速率差异都很小，基本在±4 mm 之内，说明PSI、MMSE滤波的DSI 和本文PolPSI 这3 种方法监测结果差异小，方法可靠。

图12 不同监测结果之间共同监测点形变速率差异直方图Fig. 12 Histogram of deformation rate difference of corresponding-points between different monitoring results

选取3 种方法稳定和形变区域的共同监测点P1 与P2（图11（a）中红色标示），绘制3 种方法得到时序形变，如图13 所示。两个特征点3 种方法的时序形变一致性很高，进一步说明了本文PolPSI监测方法的可靠及有效性。

图13 P1点（a）和P2点（b）时序形变图Fig. 13 Time series deformation of points P1 （a） and P2 （b）

5 结 论

为充分挖掘Sentinel-1 双极化数据信息用于地表形变监测，本文提出一种基于双极化Sentinel-1数据的PolPSI 地表形变监测技术方法。以上海浦东国际机场为例，对比分析本文PolPSI 方法的监测效果，得到以下结论：（1）本文PolPSI 技术通过对PS进行DA准则的极化优化，对DS进行自适应滤波与极化优化，能很好地保持PS 目标分辨率，提升DS 目标相位质量，在降低干涉图相位噪声的同时较好保留地物细节；（2）本文PolPSI 技术能显著提升高质量监测点密度，相比于VV 极化传统PSI 方法，监测点增加103.0%，相比基于MMSE 自适应滤波的DSI 方法，监测点增加了30.8%；在一些机场跑道区域，由于哨兵数据分辨率较低，传统方法未能得到形变监测点，而本文PolPSI 技术能探测出这些区域形变；（3）对研究区SAR 像元最优极化参数的统计分析显示，PS目标对VV极化响应较好，VH 极化对DS 目标干涉相位极化优化有较大贡献。