基于奇异值分解的DS-InSAR相位优化

基于奇异值分解的DS-InSAR相位优化 DS-InSAR phase optimization based on singular value decomposition 1 2 first-author 彭锴，研究方向为InSAR理论及沉降监测应用。E-mail： 1634612853@qq.com 彭锴，研究方向为InSAR理论及沉降监测应用。E-mail： 1634612853@qq.com 1634612853@qq.com 1 2 1 2 corresp 汪云甲，研究方向为国土环境及灾害监测、室内外无缝定位、资源环境评价及管理。E-mail： wyj4139@163.com 汪云甲，研究方向为国土环境及灾害监测、室内外无缝定位、资源环境评价及管理。E-mail： wyj4139@163.com wyj4139@163.com 1 2 1 2 中国矿业大学 环境与测绘学院, 徐州 221116 中国矿业大学 环境与测绘学院, 徐州 221116 School of Environment Science and Spatial Information, China University of Mining and Technology, Xuzhou 221116, China School of Environment Science and Spatial Information, China University of Mining and Technology, Xuzhou 221116, China 自然资源部国土环境与灾害监测重点实验室, 徐州 221116 自然资源部国土环境与灾害监测重点实验室, 徐州 221116 The Key Laboratory of Land Environment and Disaster Monitoring, MNR, China University of Mining and Technology,Xuzhou 221116, China The Key Laboratory of Land Environment and Disaster Monitoring, MNR, China University of Mining and Technology,Xuzhou 221116, China 高质量监测点空间密度是时序InSAR开展形变监测的重要指标，依靠分布式散射体DS（Distributed Scatters）开展InSAR形变监测可有效解决传统时序InSAR监测点空间密度不够的缺陷，但分布式散射体干涉相位极易受到去相干影响，造成干涉相位失真，因此分布式散射体相位优化是DS-InSAR技术的关键，针对这一情况，本文提出一种新的基于奇异值分解的DS相位优化方法，该方法利用同质像元时序相位重构相位矩阵，对矩阵进行主成分分析得到优化相位。采用模拟数据和33景覆盖郑州东部白沙镇Sentinel-1A数据对提出方法的可靠性与有效性进行验证和分析。采用时序平均相位标准偏差、平均相位梯度、平均残差点数目作为干涉图优化效果评价指标，提出方法优化后干涉图相较于原始干涉图分别下降了15.61%、25.81%、44.84%，与对比DS相位优化方法相比，提出方法对干涉图DS相位优化效果更好，特别是在一些低相干区域仍然可获得较好DS相位优化结果；提出方法在降低DS相位噪声同时可较好地保持地物细节信息。此外，相较于常规PS（Permanent Scatters）技术形变监测结果，本文方法高质量监测点数量由121471个提升至644789个，提高了4.3倍，较对比方法高质量监测点密度提升更显著。模拟与真实数据结果证实了本文提出DS优化方法的有效性，该方法可用于基于DS-InSAR技术的地表形变监测。 The spatial density of high-quality monitoring points is an important indicator for time-series InSAR to carry out deformation monitoring. Relying on Distributed Scatterers (DS) to carry out InSAR deformation monitoring can effectively solve the defect of insufficient spatial density of traditional time-series InSAR monitoring points, but the interferometric phase of distributed scatterers is easily affected by decoherence, causing the interferometric phase distortion and unreliable, so the phase optimization of distributed scatterers is the key to DS-InSAR technology and is particularly significant. Aiming at this situation, this paper proposes a new DS phase optimization method based on singular value decomposition. This method reconstructs the phase matrix by using the time-series phase of homogeneous pixels, which belong to the same substance within the inspection window, and performs principal component analysis on the matrix to obtain the optimized phase. As a very necessary step, simulation data and 33 scene coverages of Sentinel-1A data in Baisha town, eastern Zhengzhou, are used to verify the reliability and validity of the proposed method. Using time-series average phase standard deviation, average phase gradient, and average number of residual points as the evaluation index of interferogram optimization effect, the index of the interferogram optimized by the proposed method is reduced by 15.61%, 25.81%, and 44.84% respectively compared with the original interferogram, which shows that these indicators have significant decreases. The results show that, compared with the contrast DS phase optimization method, the proposed method has a better effect on the interferogram DS phase optimization, especially in some areas with poor coherence and low signal to noise ratio. In addition, the proposed method can better maintain the detailed information of the ground features while reducing the DS phase noise. Besides, compared with the deformation monitoring results of the conventional PS(Permanent Scatterers)-InSAR technology, the number of high-quality monitoring points in this method has increased from 121471 to 644789, an increase of 4.3 times, and the density of high-quality monitoring points has increased more significantly than the comparison method. The experimental results of the simulation and real data confirm the effectiveness of the DS optimization method proposed in this paper, which can be used in DS-InSAR technology for surface deformation monitoring. 遥感 DS-InSAR 分布式散射体 相位优化 奇异值分解 形变监测 remote sensing DS-InSAR distributed scatterers phase optimization singular value decomposition deformation monitoring 1　引言 1 差分合成孔径雷达干涉测量D-InSAR（Differential interferometric synthetic aperture radar）经过几十年的发展，已成为一种能大范围精确测量地表形变的对地观测技术，具有高精度、高分辨率、全天候等优点（ 朱建军 等，2017 朱建军 等，2017 ）。传统的D-InSAR技术很容易受到地形、时空失相干、大气、轨道等因素影响导致监测精度较低（ 朱建军等，2019 朱建军等，2019 ）。为克服上述因素限制，研究者们提出了时序InSAR技术。其中，最具代表性的永久散射体PS（Permanent Scatters）技术通过使用在时间维度上保持相位稳定的像元（PS目标）进行形变监测，能够较好地解决时、空去相干等问题（ Ferretti等，2001 Ferretti等，2001 ），已被广泛应用于各类形变监测中。 但永久散射体技术在建筑物较少的非城市区域，存在监测点密度低的缺点（ Cao等，2015 Cao等，2015 ）。在这些区域，分布式散射体DS（Distributed Scatters）目标较多，通常对应植被，水体、裸地等地物。因此，提取DS目标信息用于形变监测，以增加监测点密度成为近些年时序InSAR技术研究热点。然而，DS极容易受时、空去相干影响导致其干涉相位噪声较大。一般的，在和PS目标联合进行形变解算前需要对DS目标进行优化处理，处理步骤主要包括：同质像元点识别和相位优化（ 蒋弥 等，2018 蒋弥 等，2018 ）。 Ferretti等（2011） Ferretti等（2011） 提出SqueeSAR技术，采用自适应空间滤波与最大似然法获取时空优化后PS与DS的干涉相位。 Fornaro（2015） Fornaro（2015） 采用协方差矩阵进行特征值分解得到优化后像元相位。2016年， Cao等（2016） Cao等（2016） 对Fornaro方法进行改进，对相干矩阵进行特征值分解，可避免因SAR成像过程中后向散射不均衡对SAR像元优化所带来的不利影响。 上述DS优化算法大多数是通过像元协方差矩阵分解来对参考像元（待优化DS）进行相位优化，而直接基于DS同质像元时序相位构建同质像元时序相位矩阵，并对该矩阵进行分析以达到DS相位优化的研究还较少。本文基于已有的同质像元识别算法，提出一种对同质像元时序相位矩阵分解来优化DS像元相位的方法。为评价本文提出方法的有效性，分别使用模拟数据和真实数据对该方法进行测试和验证。 2　方法流程 1 本节首先简要介绍本文所采取的FaSHPS（Fast Statistically Homogeneous Pixel Selection）同质像元识别方法（ Jiang等，2015 Jiang等，2015 ）；之后，对提出的DS优化方法进行详细介绍。 2.1　DS同质像元点识别 2 采用FaSHPS同质像元点选取方法，其基本思想是将三维数据平均，将平均后的参考像素值作为真值，邻域像素值作为待估值，在给定的置信水平下，当待估值落入由真值构筑的置信区间，即视为与参考像素同质（ Jiang等，2015 Jiang等，2015 ）。根据中心极限定理，振幅样本均值 <math id="M1"><mover accent="true"><mi>A</mi><mo>¯</mo></mover><mfenced separators="|"><mrow><mi>p</mi></mrow></mfenced><mo>=</mo><mfrac><mrow><mn mathvariant="normal">1</mn></mrow><mrow><mi>N</mi></mrow></mfrac><mstyle displaystyle="true"><munderover><mo largeop="true">∑</mo><mrow><mi>i</mi><mo>=</mo><mn mathvariant="normal">1</mn></mrow><mrow><mi>N</mi></mrow></munderover></mstyle><mrow><msub><mrow><mi>A</mi></mrow><mrow><mi>i</mi></mrow></msub><mfenced separators="|"><mrow><mi>p</mi></mrow></mfenced></mrow></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557336&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557352&type= 28.10933495 8.38199997 ，在样本数量 N 增加的过程中逐渐服从高斯分布，即 <math id="M2"><mover accent="true"><mi>A</mi><mo>¯</mo></mover><mfenced separators="|"><mrow><mi>p</mi></mrow></mfenced><mo>~</mo><mi>N</mi><mfenced close="" separators="|"><mrow><mover><mrow><munder><mrow><mi>u</mi></mrow><mrow/></munder></mrow><mrow/></mover><mfenced separators="|"><mrow><mi>p</mi></mrow></mfenced></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557355&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557353&type= 22.94466591 7.19666624 ， <math id="M3"><mfenced open="" separators="|"><mrow><mrow><mrow><mi>V</mi><mi>a</mi><mi>r</mi><mfenced separators="|"><mrow><mi>A</mi><mfenced separators="|"><mrow><mi>p</mi></mrow></mfenced></mrow></mfenced></mrow><mo>/</mo><mrow><mi>N</mi></mrow></mrow></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557356&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557369&type= 22.77533340 7.02733374 ，从而获得区间解析式： <math display="block" id="M4"><mtable columnalign="left"><mtr><mtd><mi>P</mi><mfenced open="{" close="" separators="|"><mrow><msub><mrow><mi>μ</mi></mrow><mrow><mfenced separators="|"><mrow><mi>P</mi></mrow></mfenced></mrow></msub><mo>-</mo><msub><mrow><mi>z</mi></mrow><mrow><mn mathvariant="normal">1</mn><mo>-</mo><mrow><mrow><mi>α</mi></mrow><mo>/</mo><mrow><mn mathvariant="normal">2</mn></mrow></mrow></mrow></msub><mo>⋅</mo><mroot><mrow><mfrac><mrow><mi>V</mi><mi>a</mi><msub><mrow><mi>r</mi></mrow><mrow><mi>A</mi><mfenced separators="|"><mrow><mi>P</mi></mrow></mfenced></mrow></msub></mrow><mrow><mi>N</mi></mrow></mfrac></mrow><mrow/></mroot><mo><</mo><mover accent="true"><mi>A</mi><mo>¯</mo></mover><mfenced separators="|"><mrow><mi>p</mi></mrow></mfenced><mo><</mo><msub><mrow><mi>μ</mi></mrow><mrow><mfenced separators="|"><mrow><mi>P</mi></mrow></mfenced></mrow></msub><mo>+</mo></mrow></mfenced></mtd></mtr><mtr><mtd><mfenced open="" close="}" separators="|"><mrow><msub><mrow><mi>z</mi></mrow><mrow><mn mathvariant="normal">1</mn><mo>-</mo><mrow><mrow><mi>α</mi></mrow><mo>/</mo><mrow><mn mathvariant="normal">2</mn></mrow></mrow></mrow></msub><mo>⋅</mo><mroot><mrow><mfrac><mrow><mi>V</mi><mi>a</mi><msub><mrow><mi>r</mi></mrow><mrow><mi>A</mi><mfenced separators="|"><mrow><mi>P</mi></mrow></mfenced></mrow></msub></mrow><mrow><mi mathvariant="bold-italic">N</mi></mrow></mfrac></mrow><mrow/></mroot></mrow></mfenced><mo>=</mo><mn mathvariant="normal">1</mn><mo>-</mo><mi>α</mi></mtd></mtr></mtable></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557374&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557372&type= 64.17733002 22.18266678 <math id="M5"><mi mathvariant="normal">式中</mi><mo>,</mo><msub><mrow><mi>z</mi></mrow><mrow><mn mathvariant="normal">1</mn><mo>-</mo><mrow><mrow><mi>α</mi></mrow><mo>/</mo><mrow><mn mathvariant="normal">2</mn></mrow></mrow></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557358&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557375&type= 18.79599953 4.06400013 为标准正态概率密度函数的1-α/2分位数，根据SAR图像的统计理论，单视复数图的振幅在同质区域服从瑞利分布，瑞利分布的变异系数是常数，得到： <math display="block" id="M6"><mi mathvariant="normal">C</mi><mi mathvariant="normal">V</mi><mo>=</mo><mroot><mrow><mfrac><mrow><mn mathvariant="normal">4</mn></mrow><mrow><mi mathvariant="normal">π</mi></mrow></mfrac><mo>-</mo><mn mathvariant="normal">1</mn></mrow><mrow/></mroot><mo>≈</mo><mn mathvariant="normal">0.52</mn></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557378&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557359&type= 33.18933105 8.97466660 <math display="block" id="M7"><mi>V</mi><mi>a</mi><msub><mrow><mi>r</mi></mrow><mrow><mi>A</mi><mfenced separators="|"><mrow><mi>p</mi></mrow></mfenced></mrow></msub><mo>=</mo><msup><mrow><mfenced separators="|"><mrow><mi mathvariant="normal">C</mi><mi mathvariant="normal">V</mi><mo>⋅</mo><msub><mrow><mi>μ</mi></mrow><mrow><mfenced separators="|"><mrow><mi>p</mi></mrow></mfenced></mrow></msub></mrow></mfenced></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup><mo>≈</mo><msup><mrow><mfenced separators="|"><mrow><mfrac><mrow><mn mathvariant="normal">0.52</mn></mrow><mrow><mroot><mrow><mi>L</mi></mrow><mrow/></mroot></mrow></mfrac><mo>⋅</mo><msub><mrow><mi>μ</mi></mrow><mrow><mfenced separators="|"><mrow><mi>p</mi></mrow></mfenced></mrow></msub></mrow></mfenced></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557382&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557381&type= 55.88000107 11.59933376 然后得到完整的置信区间解析式： <math display="block" id="M8"><mi>P</mi><mfenced open="{" close="}" separators="|"><mrow><msub><mrow><mi>μ</mi></mrow><mrow><mi>p</mi></mrow></msub><mo>-</mo><msub><mrow><mi>z</mi></mrow><mrow><mn mathvariant="normal">1</mn><mo>-</mo><mrow><mrow><mi>α</mi></mrow><mo>/</mo><mrow><mn mathvariant="normal">2</mn></mrow></mrow></mrow></msub><mo>⋅</mo><mn mathvariant="normal">0.52</mn><mo>⋅</mo><msub><mrow><mi>μ</mi></mrow><mrow><mi>p</mi></mrow></msub><mo>∕</mo><mroot><mrow><mi>N</mi><mo>⋅</mo><mi>L</mi></mrow><mrow/></mroot><mo><</mo><mover accent="true"><mi>A</mi><mo>¯</mo></mover><mfenced separators="|"><mrow><mi>p</mi></mrow></mfenced><mo><</mo><msub><mrow><mi>μ</mi></mrow><mrow><mfenced separators="|"><mrow><mi>p</mi></mrow></mfenced></mrow></msub><mo>+</mo><msub><mrow><mi>z</mi></mrow><mrow><mn mathvariant="normal">1</mn><mo>-</mo><mrow><mrow><mi>α</mi></mrow><mo>/</mo><mrow><mn mathvariant="normal">2</mn></mrow></mrow></mrow></msub><mo>⋅</mo><mn mathvariant="normal">0.52</mn><mo>⋅</mo><mrow><mrow><msub><mrow><mi>μ</mi></mrow><mrow><mi>p</mi></mrow></msub></mrow><mo>/</mo><mrow><mroot><mrow><mi>N</mi><mo>⋅</mo><mi>L</mi></mrow><mrow/></mroot></mrow></mrow></mrow></mfenced><mo>=</mo><mn mathvariant="normal">1</mn><mo>-</mo><mi>α</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557386&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557384&type= 129.96333313 6.43466663 式中， L 为视数。 2.2　DS相位优化 2 在同质像元识别准确情况下，识别窗口内同质像元差分干涉相位基本一致（已去除平地和地形相位）。基于此，本文提出运用奇异值分解SVD（Singular Value Decomposition）主成分分析（ Kalman，1996 Kalman，1996 ）思想对DS像元相位基于其同质像元进行优化。 具体步骤如 图1 图1 所示。首先，在遍历过程中的每一个窗口内，将同质像元时序相位按列排的顺序依次（参考像元序号为 n ，如 图1 图1 中为 n =3）横向展开到二维平面构成矩阵 D ，即矩阵 D 的每一行为所识别一个同质像元的时序相位，则 D 矩阵大小为 M <math id="M9"><mo>×</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557913&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557912&type= 1.77800000 2.62466669 N ，其中 M 为参考像元（待优化DS像元）的同质像元数目， N 为时序干涉图幅数。然后，对 D 矩阵进行奇异值分解： <math display="block" id="M10"><mi mathvariant="bold-italic">D</mi><mo>=</mo><mstyle displaystyle="true"><munderover><mo largeop="true">∑</mo><mrow><mi>i</mi><mo>=</mo><mn mathvariant="normal">1</mn></mrow><mrow><mi>r</mi></mrow></munderover></mstyle><mrow><msub><mrow><mi mathvariant="bold-italic">s</mi></mrow><mrow><mi mathvariant="bold-italic">i</mi></mrow></msub><msub><mrow><mi mathvariant="bold-italic">u</mi></mrow><mrow><mi mathvariant="bold-italic">i</mi></mrow></msub><msubsup><mrow><mi mathvariant="bold-italic">v</mi></mrow><mrow><mi mathvariant="bold-italic">i</mi></mrow><mrow><mi mathvariant="normal">T</mi></mrow></msubsup></mrow></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557393&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557366&type= 19.72733307 8.29733276 式中， r 为矩阵 D 的非零奇异值个数，且 r ≤min（ M , N ）， <math id="M11"><msub><mrow><mi mathvariant="bold-italic">u</mi></mrow><mrow><mi>i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557396&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557400&type= 2.70933342 3.72533321 和 <math id="M12"><msub><mrow><mi mathvariant="bold-italic">v</mi></mrow><mrow><mi>i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557406&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557397&type= 2.37066650 3.72533321 分别是 D 的奇异值分解中左、右正交矩阵的列向量。利用最大奇异值 <math id="M13"><msub><mrow><mi>s</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557419&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557416&type= 2.45533323 3.72533321 及其所对应的左右奇异向量 <math id="M14"><msub><mrow><mi>u</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557410&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557420&type= 3.04800010 3.72533321 和 <math id="M15"><msub><mrow><mi>v</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557428&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557424&type= 2.70933342 3.72533321 构成降噪后的 D 矩阵 <math id="M16"><msup><mrow><mi mathvariant="bold-italic">D</mi></mrow><mrow><mi mathvariant="normal">'</mi></mrow></msup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557432&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557431&type= 3.38666677 2.96333337 ，如 式（6） 式（6） 所示： <math display="block" id="M17"><msup><mrow><mi mathvariant="bold-italic">D</mi></mrow><mrow><mi mathvariant="bold-italic">'</mi></mrow></msup><mo>=</mo><msub><mrow><mi mathvariant="bold-italic">s</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><msub><mrow><mi mathvariant="bold-italic">u</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><msubsup><mrow><mi mathvariant="bold-italic">v</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow><mrow><mi mathvariant="normal">T</mi></mrow></msubsup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557437&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557434&type= 16.93333435 4.06400013 最后，对遍历过程中的每一个矩阵 <math id="M18"><msup><mrow><mi mathvariant="bold-italic">D</mi></mrow><mrow><mi>'</mi></mrow></msup></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557444&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557440&type= 3.38666677 2.96333337 中提取出第 n 行作为新的参考像元时序相位值构成新的优化干涉图。 对分解后的干涉图进行相干性估计（ Jiang等，2014 Jiang等，2014 ），并利用相干性对DS同质像元进行加权平均得到最终的DS优化干涉相位如 式（7） 式（7） 所示： <math display="block" id="M19"><msub><mrow><mi>φ</mi></mrow><mrow><mi>r</mi><mi>e</mi><mi>f</mi></mrow></msub><mo>=</mo><mfrac><mrow><mstyle displaystyle="true"><munderover><mo largeop="true">∑</mo><mrow><mi>i</mi><mo>=</mo><mn mathvariant="normal">1</mn></mrow><mrow><mi>n</mi></mrow></munderover></mstyle><mrow><msub><mrow><mi>φ</mi></mrow><mrow><mi>i</mi></mrow></msub><mo>·</mo><mi>C</mi><mi>o</mi><msub><mrow><mi>h</mi></mrow><mrow><mi>i</mi></mrow></msub></mrow></mrow><mrow><mstyle displaystyle="true"><munderover><mo largeop="true">∑</mo><mrow><mi>i</mi><mo>=</mo><mn mathvariant="normal">1</mn></mrow><mrow><mi>n</mi></mrow></munderover></mstyle><mrow><mi>C</mi><mi>o</mi><msub><mrow><mi>h</mi></mrow><mrow><mi>i</mi></mrow></msub></mrow></mrow></mfrac></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557451&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557463&type= 24.97666740 17.86466599 式中， <math id="M20"><msub><mrow><mi>φ</mi></mrow><mrow><mi>r</mi><mi>e</mi><mi>f</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557470&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557467&type= 4.40266657 4.40266657 表示参考像元的相位值， <math id="M21"><msub><mrow><mi>φ</mi></mrow><mrow><mi>i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557453&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557475&type= 2.96333337 4.40266657 代表周围同质像元点的相位值， <math id="M22"><mi>C</mi><mi>o</mi><msub><mrow><mi>h</mi></mrow><mrow><mi>i</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557480&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557478&type= 6.51933336 3.72533321 为估计得到的相干性。 2.3　地表形变反演数据处理流程 2 首先，采用单幅主影像配准，结合SRTM DEM（30 m）生成差分干涉图并进行Goldstein滤波（ Goldstein和Werner，1998 Goldstein和Werner，1998 ）处理。其次，采用FaSHPS进行同质像元识别和PS与DS分类（ Ferretti等，2011 Ferretti等，2011 ； Jiang等，2015 Jiang等，2015 ），对DS像元，基于其同质像元构建矩阵 D 并进行奇异值分解。随后，基于分解处理后的同质像元相位进行加权平均滤波，得到优化的DS像元相位。最后采用TPC（Temporal Phase Coherence）方法设置一定阈值选取高质量像元（ Zhao和Mallorqui，2019 Zhao和Mallorqui，2019 ），利用StaMPS方法基于优化后干涉图得到地表形变监测结果（ Hooper等，2004 Hooper等，2004 ， 2007 2007 ； Hooper和Zebker，2007 Hooper和Zebker，2007 ），整体流程图如 图2 图2 所示。 3　模拟数据实验 1 首先，详细介绍时序干涉图模拟步骤；随后，分别采用同质像元相干性加权平均法和本文提出方法对模拟数据进行相位优化并作对比分析。为叙述方便，下文对“同质像元相干性加权平均法”统一叙述为“加权平均法”。 3.1　模拟步骤 2 模拟数据只模拟DS点目标，所使用的原始图片从 http://www.perso.telecom-paristech.fr/~deledall/nlinsar.php http://www.perso.telecom-paristech.fr/~deledall/nlinsar.php ［2020-08-15］获取。 基于原始图片可模拟真实的强度 <math id="M23"><mi>R</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557576&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557583&type= 2.37066650 2.62466669 ，相位 <math id="M24"><mi>β</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557611&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557608&type= 2.28600001 3.38666677 以及相干性 <math id="M25"><mi>γ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557617&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557614&type= 2.03200006 3.38666677 （无噪声）。然后根据模拟的真实的强度 <math id="M26"><mi>R</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557576&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557583&type= 2.37066650 2.62466669 ，相位 <math id="M27"><mi>β</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557611&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557608&type= 2.28600001 3.38666677 以及相干性 <math id="M28"><mi>γ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557617&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557614&type= 2.03200006 3.38666677 ，模拟对应的SLC图像 <math id="M29"><msub><mrow><mi>z</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557537&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557536&type= 2.53999996 3.72533321 和 <math id="M30"><msub><mrow><mi>z</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557642&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557629&type= 2.53999996 3.72533321 （加噪声）。具体模拟的过程如下（ Deledalle等，2011 Deledalle等，2011 ； Lee等，2003 Lee等，2003 ）： （1）首先计算矩阵 <math id="M31"><mi mathvariant="bold-italic">L</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557548&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557527&type= 2.03200006 2.62466669 使得 <math id="M32"><mi mathvariant="bold">C</mi><mi mathvariant="bold">V</mi><mo>=</mo><mi mathvariant="bold-italic">L</mi><mi>*</mi><mi mathvariant="bold-italic">L</mi><mi>'</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557530&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557552&type= 16.51000023 2.70933342 ， <math id="M33"><mi mathvariant="bold">C</mi><mi mathvariant="bold">V</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557565&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557563&type= 4.99533367 2.70933342 是模拟干涉图协方差矩阵，如 式（5） 式（5） 所示： <math display="block" id="M34"><mtable><mtr><mtd><mrow><mi mathvariant="bold">C</mi><mi mathvariant="bold">V</mi><mo>=</mo><mfenced open="[" close="]" separators="|"><mrow><mtable><mtr><mtd><mi>R</mi></mtd><mtd><mroot><mrow><mi>R</mi><msup><mrow><mi>R</mi></mrow><mrow><mi>'</mi></mrow></msup></mrow><mrow/></mroot><mi>γ</mi><msup><mrow><mi>e</mi></mrow><mrow><mi>i</mi><mi>β</mi></mrow></msup></mtd></mtr><mtr><mtd><mroot><mrow><mi>R</mi><msup><mrow><mi>R</mi></mrow><mrow><mi>'</mi></mrow></msup></mrow><mrow/></mroot><mi>γ</mi><msup><mrow><mi>e</mi></mrow><mrow><mo>-</mo><mi>i</mi><mi>β</mi></mrow></msup></mtd><mtd><msup><mrow><mi>R</mi></mrow><mrow><mi>'</mi></mrow></msup></mtd></mtr></mtable></mrow></mfenced></mrow></mtd></mtr></mtable></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557571&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557569&type= 45.97399902 12.27666664 通过Cholesky分解得到的 L 如下： <math display="block" id="M35"><mtable><mtr><mtd><mrow><mi mathvariant="bold-italic">L</mi><mo>=</mo><mfenced open="[" close="]" separators="|"><mrow><mtable><mtr><mtd><mroot><mrow><mi>R</mi></mrow><mrow/></mroot></mtd><mtd><mn mathvariant="normal">0</mn></mtd></mtr><mtr><mtd><mi>γ</mi><mroot><mrow><msup><mrow><mi>R</mi></mrow><mrow><mi>'</mi></mrow></msup></mrow><mrow/></mroot><msup><mrow><mi>e</mi></mrow><mrow><mo>-</mo><mi>i</mi><mi>β</mi></mrow></msup></mtd><mtd><mroot><mrow><mfenced separators="|"><mrow><mn mathvariant="normal">1</mn><mo>-</mo><msup><mrow><mi>γ</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msup></mrow></mfenced><msup><mrow><mi>R</mi></mrow><mrow><mi>'</mi></mrow></msup></mrow><mrow/></mroot></mtd></mtr></mtable></mrow></mfenced></mrow></mtd></mtr></mtable></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557581&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557579&type= 46.22800064 13.37733269 式中， <math id="M36"><mi>R</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557576&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557583&type= 2.37066650 2.62466669 和 <math id="M37"><mi>R</mi><mi>'</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557596&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557589&type= 3.55599999 2.70933342 代表主、辅影像的强度， <math id="M38"><mi>β</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557611&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557608&type= 2.28600001 3.38666677 是干涉相位， <math id="M39"><mi>γ</mi></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557617&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557614&type= 2.03200006 3.38666677 是相干性。 （2）生成两个服从零均值复圆高斯分布的独立随机复数变量 <math id="M40"><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557621&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557601&type= 2.87866688 3.72533321 和 <math id="M41"><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557623&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557605&type= 2.87866688 3.72533321 ，然后通过 式（10） 式（10） 即可得到模拟干涉图所对应的单视复数影像 <math id="M42"><msub><mrow><mi>z</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557641&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557625&type= 2.53999996 3.72533321 和 <math id="M43"><msub><mrow><mi>z</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557642&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557629&type= 2.53999996 3.72533321 ： <math display="block" id="M44"><mfenced open="[" close="]" separators="|"><mrow><mtable><mtr><mtd><msub><mrow><mi>z</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi>z</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></mtd></mtr></mtable></mrow></mfenced><mo>=</mo><mi>L</mi><mfenced open="[" close="]" separators="|"><mrow><mtable><mtr><mtd><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></mtd></mtr><mtr><mtd><msub><mrow><mi>x</mi></mrow><mrow><mn mathvariant="normal">2</mn></mrow></msub></mtd></mtr></mtable></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557635&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557643&type= 18.45733452 8.21266651 （3）重复上述两个步骤，生成时序干涉图集和对应的时序SLC影像。所得到干涉图的噪声水平与给定相干性相关，为尽量接近真实情况，本文模拟的干涉图相干性随着影像时间间隔的增加而减少。 （4）考虑到需要依据SAR强度做同质像元识别，本文模拟的强度随时间变化。不失一般性，强度图模拟为周期性变化： <math display="block" id="M45"><msub><mrow><mi>R</mi></mrow><mrow><mi>n</mi></mrow></msub><mo>=</mo><msub><mrow><mi>R</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mo>+</mo><mi mathvariant="normal">s</mi><mi mathvariant="normal">i</mi><mi mathvariant="normal">n</mi><mtext> </mtext><mfenced separators="|"><mrow><mfrac><mrow><mn mathvariant="normal">2</mn><mi mathvariant="normal">π</mi><mi mathvariant="normal">n</mi></mrow><mrow><mi>N</mi></mrow></mfrac></mrow></mfenced><mo>×</mo><msub><mrow><mi>R</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub><mfenced separators="|"><mrow><mn mathvariant="normal">1</mn><mo>-</mo><mi>γ</mi></mrow></mfenced></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557637&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557636&type= 50.20733261 8.46666718 式中， <math id="M46"><msub><mrow><mi>R</mi></mrow><mrow><mn mathvariant="normal">1</mn></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557657&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557654&type= 3.47133350 3.72533321 和 <math id="M47"><msub><mrow><mi>R</mi></mrow><mrow><mi>n</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557660&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557672&type= 3.47133350 3.72533321 表示第一幅和第 n 幅SLC图像的强度图， N 代表模拟影像的总数。 3.2　模拟结果 2 经过上述过程即可得到模拟的时序干涉图，将 式（11） 式（11） 中的 N 值设为33，使之与真实数据SLC数目相同。 图3 图3 为第一幅模拟干涉图所对应的真实相位、相干性及强度。需要注意的是，在模拟的过程中，由于实验中最大相干性设置为0.8，最小设置为0.4，并从主影像开始随时间基线增加逐渐衰减，故第一幅模拟干涉图真实相干性整体最低。为模拟更多不同干涉相位情况，在后3组模拟实验时只改变干涉相位，强度和相干性均与第1组模拟相同。 1,2,3;4,5,6; 对4组模拟干涉图采用加权平均（ Jiang等，2014 Jiang等，2014 ）和本文提出的方法处理，得到优化的相位结果如 图4 图4 所示。由于篇幅原因，这里仅展示各组第一幅模拟干涉图处理结果。从图中可看出，在相干性较高的区域，加权平均法能得到较好的结果，但在相干性较低区域却得不到理想的优化相位；相比之下，本文提出的方法在相干性高或者较低区域均能得到更接近真实值的DS相位。 1,2,3,4;5; 为进一步整体地评价各方法DS点相位优化效果，计算了各组像元优化相位中误差平均值（ 表1 表1 ）；并且绘制了第250行、250列的相位剖面图和第300行、300列的点时序相位图。由于篇幅原因，这里的仅展示第一幅干涉图的相位剖面图，时序中误差图与点时序相位图对应第一组模拟的结果，如 图5 图5 所示。 第1组 第2组 第3组 第4组 同质像元相干性加权平均 0.7726 0.7421 0.7466 0.7383 本文方法 0.2401 0.2068 0.2083 0.2068 1,2;3,4; 与 图4 图4 结果相符，在相干性较高的区域，两种优化方法中误差值都较低，但在相干性低的区域，本文方法表现更优，且中误差整体保持较低水平（ 图5（a） 图5（a） ）。由 表1 表1 可知，提出方法中误差平均值均低于加权平均方法所得到的结果，说明本文方法的优化相位与真实值更为接近，DS相位优化效果更好。 点时序相位图中（ 图5（b） 图5（b） ），加权平均方法结果虽然与真实相位趋势一致，但在部分干涉图中，优化值偏离真值较远；而本文方法结果几乎与真实值重合。在所选取的行、列相位剖面图中（ 图5（c） 图5（c） 、（d）），加权平均方法优化相位的噪声并没有得到很好地抑制，与模拟加噪相位一同散乱分布；本文提出方法优化相位不仅与真实相位趋势一致，且都以较小的波动幅度分布在真实值周围，异常值数量更少。综上所述，模拟实验结果证明了本文方法DS像元相位优化结果的准确性和可靠性。 4　真实数据实验 1 除使用模拟数据验证本文提出DS相位优化方法外，同时基于33景覆盖郑州的VV极化Sentinel-1A数据，利用提出方法对该地区进行形变监测以验证方法有效性。选择郑州城区东部白沙镇作为研究区域，影像大小4400×1500，分辨率为5 m×20 m。以2018年6月13日获取的SAR数据为主影像，生成32景干涉图，具体时空基线信息如 图6 图6 所示。 FaSHPS同质像元识别采用15 <math id="M48"><mo>×</mo></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557913&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557912&type= 1.77800000 2.62466669 15窗口，结果如 图7 图7 右侧所示。 图7 图7 左侧为Google Earth研究区影像，为了与SAR影像保持一致，对其做上下颠倒处理。同质像元识别结果地物边界清晰，同质区域与实际情况相符。 4.1　差分干涉图优化结果分析 2 选取原始干涉图相位及两种优化方法得到的优化相位做分析， 图8 图8 分别为时间基线最长的第一幅和最短的第16幅优化前后的干涉相位图。从左至右依次为原始干涉图相位、加权平均结果、本文方法结果。与模拟数据结果相似，在噪声水平较低的区域，加权平均能较好地优化干涉相位，但在低相干区域则表现不佳；相比之下，本文提出方法即使在部分相干性较低的区域也能得到较为理想的相位优化结果。例如，在第一幅优化干涉图右下方的道路条纹（ 图8 图8 黑色矩形所示）更加清晰，整体相位条纹更加平滑，质量更好。 1;2;3; 对所有干涉图集计算时序平均相位标准偏差、平均相位梯度、平均残差点数目作为干涉图优化效果评价指标（ Goldstein和Werner，1998 Goldstein和Werner，1998 ； Li等，2004 Li等，2004 ），结果如 表2 表2 所示。值得说明的是，本文使用的这3种指标值越小，对应干涉图集相位噪声水平越低。可看到，在各项指标上，加权平均方法结果相较于原始干涉图分别下降了4.3%、6.66%、12.08%，本文方法结果则分别下降了15.61%、25.81%、44.84%，较加权平均方法都有大幅度的提升，证明了本文方法的有效性。 相位梯度均值（弧度） 相位偏差均值（弧度） 残差点均值（个） 原始干涉图 1.6644 1.7966 1347700 同质点相干性加权平均 1.5929（4.3%↓） 1.6770（6.66%↓） 1184884（12.08%↓） 本文方法 1.4046（15.61%↓） 1.3315（25.81%↓） 743440（44.84%↓） 4.2　时序沉降优化结果分析 2 需要说明的是，与模拟数据处理不同，为尽量保留PS点相位信息，不对同质像元数量小于 <math id="M49"><msub><mrow><mi>T</mi></mrow><mrow><mi mathvariant="normal">S</mi><mi mathvariant="normal">H</mi><mi mathvariant="normal">P</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557991&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557998&type= 5.75733376 3.72533321 的像元进行处理。在本次实验中，阈值 <math id="M50"><msub><mrow><mi>T</mi></mrow><mrow><mi mathvariant="normal">S</mi><mi mathvariant="normal">H</mi><mi mathvariant="normal">P</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557991&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42557998&type= 5.75733376 3.72533321 为20（ Ferretti等，2011 Ferretti等，2011 ）。 基于原始干涉图和两种方法优化后的干涉图，统一采用TPC筛选高质量点目标。阈值 <math id="M51"><msub><mrow><mi>T</mi></mrow><mrow><mi mathvariant="normal">t</mi><mi mathvariant="normal">h</mi><mi mathvariant="normal">r</mi><mi mathvariant="normal">e</mi></mrow></msub></math> http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558003&type= http://html.publish.founderss.cn/rc-pub/api/common/picture?pictureId=42558000&type= 5.33400011 3.72533321 定为0.9，对应相位标准偏差约为25°（ Zhao和Mallorqui，2019 Zhao和Mallorqui，2019 ）。随后，分别基于原始和优化后的干涉图及选取的高质量像元，使用StaMPS做时序分析，得到结果如 图9 图9 所示。 1,2;3; 图9（a） 图9（a） —（c）分别为原始干涉图、加权平均以及本文方法优化干涉图所对应的沉降结果，其最大年平均沉降速率分别为-32.4 mm/a，-33.3 mm/a，-31.8 mm/a，结果接近，并且实验区域整体沉降趋势一致，可从侧面说明监测结果的可靠性。但各方法获得的监测点数量相差较大，其中利用原始干涉图得到的形变监测点共121471个；利用加权平均优化干涉图得到的测点共165743个，较原始干涉图结果仅提升36%，点密度都相对较低，直接通过目视解译很难对研究区域整体沉降范围及趋势做出准确的判断；而采用提出方法得到形变监测点共644789个，较原始干涉图结果提升430%，整体趋势更加明朗，可更好反演研究区地表形变细节，且对准确判断具体沉降范围带来了极大的方便。 图9 图9 中白色框为研究区主要沉降区，右下角编号与主图相同的为对应放大图像。可以看到，本文方法在主要沉降区的监测点较多，较其他方法结果细节信息更为丰富，对该区域做沉降面积统计，分别约为1.25 km 2 ，1.69 km 2 和6.20 km 2 ，新方法监测的沉降范围更广。并且在本文方法沉降结果图中的右上方和左下方，可以看到道路和河流堤坝上的监测点提高较多，表明提出方法在道路、桥梁、堤坝、机场跑道等线形地物形变监测上效果更佳。 5　结论 1 针对传统PS-InSAR在非城市区域形变监测点密度较低问题，为提取出SAR影像中DS目标信息用于地表形变监测，本文提出一种新的DS相位优化方法。该方法主要思想是基于同质像元相位一致性，将同质像元时序相位展开到二维平面构成矩阵 D ，利用奇异值分解对 D 矩阵做主成分分析提取出参考像元（待优化DS像元）优化时序相位，最后再对分解后的干涉图做同质像元空间加权平均得到最终DS相位优化结果。将提出方法分别应用于模拟与真实SAR数据，结果表明：对于模拟数据和真实数据DS像元相位优化，本文方法较对比DS优化方法（加权平均法）具有更好表现。对于地表形变监测，3种形变监测方法得到的结果形变整体趋势一致，说明本文方法形变反演结果可靠；但在监测点数量上本文方法较基于原始干涉图的PS方法提升了430%，且比对比方法（36%）提升幅度更大，能够更好地反演地表形变。 参考文献（References） Cao N , Lee H and Jung H C . 2015 . Mathematical framework for phase-triangulation algorithms in distributed-scatterer interferometry . IEEE Geoscience and Remote Sensing Letters , 12 ( 9 ): 1838 - 1842 [ DOI: 10.1109/LGRS.2015.2430752 http://dx.doi.org/10.1109/LGRS.2015.2430752 ] Mathematical framework for phase-triangulation algorithms in distributed-scatterer interferometry Cao N , Lee H and Jung H C . 2016 . A phase-decomposition-based PSInSAR processing method . IEEE Transactions on Geoscience and Remote Sensing , 54 ( 2 ): 1074 - 1090 [ DOI: 10.1109/TGRS.2015.2473818 http://dx.doi.org/10.1109/TGRS.2015.2473818 ] A phase-decomposition-based PSInSAR processing method Deledalle C A , Denis L and Tupin F . 2011 . NL-InSAR: nonlocal interferogram estimation . IEEE Transactions on Geoscience and Remote Sensing , 49 ( 4 ): 1441 - 1452 [ DOI: 10.1109/TGRS.2010.2076376 http://dx.doi.org/10.1109/TGRS.2010.2076376 ] NL-InSAR: nonlocal interferogram estimation Ferretti A , Fumagalli A , Novali F , Prati C , Rocca F and Rucci A . 2011 . A new algorithm for processing interferometric Data-Stacks: SqueeSAR . IEEE Transactions on Geoscience and Remote Sensing , 49 ( 9 ): 3460 - 3470 [ DOI: 10.1109/TGRS.2011.2124465 http://dx.doi.org/10.1109/TGRS.2011.2124465 ] A new algorithm for processing interferometric Data-Stacks: SqueeSAR Ferretti A , Prati C and Rocca F . 2001 . Permanent scatterers in SAR interferometry . IEEE Transactions on Geoscience and Remote Sensing , 39 ( 1 ): 8 - 20 [ DOI: 10.1109/36.898661 http://dx.doi.org/10.1109/36.898661 ] Permanent scatterers in SAR interferometry Fornaro G , Verde S , Reale D and Pauciullo A . 2015 . CAESAR: an approach based on covariance matrix decomposition to improve multibaseline-multitemporal interferometric SAR processing . IEEE Transactions on Geoscience and Remote Sensing , 53 ( 4 ): 2050 - 2065 [ DOI: 10.1109/TGRS.2014.2352853 http://dx.doi.org/10.1109/TGRS.2014.2352853 ] CAESAR: an approach based on covariance matrix decomposition to improve multibaseline-multitemporal interferometric SAR processing Goldstein R M and Werner C L . 1998 . Radar interferogram filtering for geophysical applications . Geophysical Research Letters , 25 ( 21 ): 4035 - 4038 [ DOI: 10.1029/1998GL900033 http://dx.doi.org/10.1029/1998GL900033 ] Radar interferogram filtering for geophysical applications Hooper A , Segall P and Zebker H . 2007 . Persistent scatterer interferometric synthetic aperture radar for crustal deformation analysis, with application to Volcán Alcedo, Galápagos . Journal of Geophysical Research , 112 ( B7 ): B 07407 [ DOI: 10.1029/2006jb004763 http://dx.doi.org/10.1029/2006jb004763 ] Persistent scatterer interferometric synthetic aperture radar for crustal deformation analysis, with application to Volcán Alcedo, Galápagos Hooper A , Zebker H , Segall P and Kampes B . 2004 . A new method for measuring deformation on volcanoes and other natural terrains using InSAR persistent scatterers . Geophysical Research Letters , 31 ( 23 ): L 23611 [ DOI: 10.1029/2004GL021737 http://dx.doi.org/10.1029/2004GL021737 ] A new method for measuring deformation on volcanoes and other natural terrains using InSAR persistent scatterers Hooper A and Zebker H A . 2007 . Phase unwrapping in three dimensions with application to InSAR time series . Journal of the Optical Society of America A , 24 ( 9 ): 2737 - 2747 [ DOI: 10.1364/JOSAA.24.002737 http://dx.doi.org/10.1364/JOSAA.24.002737 ] Phase unwrapping in three dimensions with application to InSAR time series Jiang M , Ding X L , Hanssen R F , Malhotra R and Chang L . 2015 . Fast statistically homogeneous pixel selection for covariance matrix estimation for multitemporal InSAR . IEEE Transactions on Geoscience and Remote Sensing , 53 ( 3 ): 1213 - 1224 [ DOI: 10.1109/TGRS.2014.2336237 http://dx.doi.org/10.1109/TGRS.2014.2336237 ] Fast statistically homogeneous pixel selection for covariance matrix estimation for multitemporal InSAR Jiang M , Ding X L and Li Z W . 2014 . Hybrid approach for unbiased coherence estimation for multitemporal InSAR . IEEE Transactions on Geoscience and Remote Sensing , 52 ( 5 ): 2459 - 2473 [ DOI: 10.1109/TGRS.2013.2261996 http://dx.doi.org/10.1109/TGRS.2013.2261996 ] Hybrid approach for unbiased coherence estimation for multitemporal InSAR Kalman D . 1996 . A singularly valuable decomposition: the SVD of a matrix . The College Mathematics Journal , 27 ( 1 ): 1 - 23 A singularly valuable decomposition: the SVD of a matrix Lee J S , Cloude S R , Papathanassiou K P , Grunes M R and Woodhouse I H . 2003 . Speckle filtering and coherence estimation of polarimetric SAR interferometry data for forest applications . IEEE Transactions on Geoscience and Remote Sensing , 41 ( 10 ): 2254 - 2263 [ DOI: 10.1109/TGRS.2003.817196 http://dx.doi.org/10.1109/TGRS.2003.817196 ] Speckle filtering and coherence estimation of polarimetric SAR interferometry data for forest applications Li Z L , Zou W B , Ding X L , Chen Y Q and Liu G X . 2004 . A quantitative measure for the quality of INSAR interferograms based on phase differences . Photogrammetric Engineering and Remote Sensing , 70 ( 10 ): 1131 - 1137 [ DOI: 10.14358/PERS.70.10.1131 http://dx.doi.org/10.14358/PERS.70.10.1131 ] A quantitative measure for the quality of INSAR interferograms based on phase differences Zhao F and Mallorqui J J . 2019 . A temporal phase coherence estimation algorithm and its application on DInSAR pixel selection . IEEE Transactions on Geoscience and Remote Sensing , 57 ( 11 ): 8350 - 8361 [ DOI: 10.1109/TGRS.2019.2920536 http://dx.doi.org/10.1109/TGRS.2019.2920536 ] A temporal phase coherence estimation algorithm and its application on DInSAR pixel selection Jiang M , Ding X L and Li Z W . 2018 . Homogeneous pixel selection algorithm for multitemporal InSAR . Chinese Journal of Geophysics , 61 ( 12 ): 4767 - 4776 Homogeneous pixel selection algorithm for multitemporal InSAR 蒋弥 , 丁晓利 , 李志伟 . 2018 . 时序InSAR同质样本选取算法研究 . 地球物理学报 , 61 ( 12 ): 4767 - 4776 [ DOI: 10.6038/cjg2018L0490 http://dx.doi.org/10.6038/cjg2018L0490 ] 时序InSAR同质样本选取算法研究 Zhu J J , Li Z W and Hu J . 2017 . Research progress and methods of InSAR for deformation monitoring . Acta Geodaetica et Cartographica Sinica , 46 ( 10 ): 1717 - 1733 Research progress and methods of InSAR for deformation monitoring 朱建军 , 李志伟 , 胡俊 . 2017 . InSAR变形监测方法与研究进展 . 测绘学报 , 46 ( 10 ): 1717 - 1733 [ DOI: 10.11947/j.AGCS.2017.20170350 http://dx.doi.org/10.11947/j.AGCS.2017.20170350 ] InSAR变形监测方法与研究进展 Zhu J J , Yang Z F and Li Z W . 2019 . Recent progress in retrieving and predicting mining-induced 3D displace-ments using InSAR . Acta Geodaetica et Cartographica Sinica , 48 ( 2 ): 135 - 144 Recent progress in retrieving and predicting mining-induced 3D displace-ments using InSAR 朱建军 , 杨泽发 , 李志伟 . 2019 . InSAR矿区地表三维形变监测与预计研究进展 . 测绘学报 , 48 ( 2 ): 135 - 144 [ DOI: 10.11947/j.AGCS.2019.20180188 http://dx.doi.org/10.11947/j.AGCS.2019.20180188 ] InSAR矿区地表三维形变监测与预计研究进展 致谢:感谢欧洲空间局提供的免费Sentinel-1A数据以及蒋弥教授提供的同质像元选取开源工具箱。 10.11834/jrs.20210454 P237 1007-4619（2023）02-0533-10 彭锴，赵峰，汪云甲，闫世勇，冯瀚.2023.基于奇异值分解的DS-InSAR相位优化.遥感学报，27（2）： 533-542 彭锴，赵峰，汪云甲，闫世勇，冯瀚.2023.基于奇异值分解的DS-InSAR相位优化.遥感学报，27（2）： 533-542 Peng K，Zhao F，Wang Y J，Yan S Y and Feng H. 2023. DS-InSAR phase optimization based on singular value decomposition. National Remote Sensing Bulletin， 27（2）：533-542 Peng K，Zhao F，Wang Y J，Yan S Y and Feng H. 2023. DS-InSAR phase optimization based on singular value decomposition. National Remote Sensing Bulletin， 27（2）：533-542 国家自然科学基金（编号：41874044，42004011）；国家博士后面上基金（编号：2020M671646）；中央高校基本业务费（编号：2020QN27）；江苏省优势学科：测绘科学与技术 Supported by National Natural Science Foundation of China (No.41874044,42004011); China Postdoctoral Science Foundation (No. 2020M671646); Fundamental Research Funds for the Central Universities of China (No. 2020QN27); Priority Academic Program Development (PAPD) of Jiangsu Higher Education Institutions (Science and Technology of Surveying and Mapping) 2020-10-27 2023-02-07 2023-04-20 2021-03-11 遥感学报 2 27