基于改进低秩张量分解的小数据集PSInSAR地表形变估计方法

Small dataset PSInSAR surface deformation estimation method based on improved low-rank tensor decomposition

摘要

针对PSInSAR技术在小数据集下地表形变估计精度不理想的问题，提出双加权低秩张量分解算法对时序干涉相位的观测数据去噪，使其更适用于较大城市场景，增强了张量分解在小数据集PSInSAR技术上的实用性。以天津地区的29景TerraSAR图像为数据源，抽取其中6～11景作为小数据集对所提算法进行验证，采用双加权张量分解与PSInSAR结合的技术估计地面沉降。实验结果显示，所提算法对小数据集PSInSAR技术估计的地表形变质量有明显提升效果。在11景图像条件下的细节形变区域，与原低秩张量分解时序算法相比，形变速率估计的误差降低约82%。

Abstract

In addressing the issue of inadequate precision in surface deformation estimation using the PSInSAR technique with small datasets， this study proposes a dual-weighted low-rank tensor decomposition algorithm for denoising the observed data of temporal interferometric phases. This algorithm enhances the applicability of tensor decomposition in the context of small dataset PSInSAR technology， particularly in larger urban areas. This article utilizes a dataset comprising 29 TerraSAR images from the Tianjin region and extracts a subset of 6 to 11 images as a small dataset for validating the proposed algorithm. The combined approach of dual-weighted tensor decomposition and PSInSAR is employed for ground subsidence estimation. Experimental results demonstrate a significant improvement in the quality of surface deformation estimation using the proposed algorithm for small dataset PSInSAR technology. In the region exhibiting detailed deformation under the condition of 11 images， the error in deformation rate estimation is reduced by approximately 82% compared to the results obtained using the original low-rank tensor decomposition algorithm with the same number of SAR images.

Graphical abstract

关键词

PSInSAR / 小数据集 / 低秩张量分解 / 地面沉降

Key words

PSInSAR / small datasets / low rank tensor decomposition / ground settlement

近十几年来，由于我国城市化过程发展十分迅速，出现了地下资源开发过度及软土地压实等问题，导致城市地面沉降灾害频繁发生且规模越来越大^［1］，给居民生活、交通、环境等很多方面带来不便^［2］。为了能更好地控制地面沉降问题，尽量降低其危害，形变监测是必不可少的手段。目前为止，主要的传统沉降监测技术有水准测量和布设全球定位系统等，尽管这些传统的测量技术能满足精度的要求，却存在测量点密度低、成本高等劣势问题，若要满足大范围的连续形变测量需要大量的人力物力^［3］。

合成孔径雷达干涉（interferometric synthetic aperture radar，InSAR）测量技术成功克服了以上传统技术存在的缺陷，能够满足大范围内的全天时、全天候的形变监测需求^［4］。该技术从早期的差分InSAR^［5］发展到后来的时序InSAR，其中包括永久散射体干涉（permanent scatterers InSAR， PSInSAR）技术^［6］、小基线集干涉技术^［7］等，在很大程度上解决了时空失相干和大气延迟的影响，精度上也已达到mm级别^［8］，被广泛应用于城市地表沉降等的形变测量中。到现在为止，国内外对于时序InSAR的研究已经处于一个相对成熟的阶段，然而传统的PSInSAR技术仍然存在一定的难点^［9］，其中之一就是对图像数量的需求，一般而言至少也需要25景合成孔径雷达（synthetic aperture radar，SAR）的单视复影像（single look complex，SLC）才能反演出高精度的形变结果。

尽管欧洲航天局发射的哨兵1号卫星可提供免费开放、可持续的地球观测数据，但其分辨率较低，无法满足所有的研究需求。而现有的高分辨率SAR图像数据，例如TerraSAR和COSMO-SkyMed，数据资源价格昂贵，不利于国内时序InSAR的科研发展。另外大量图像的累积是以时间为代价的，以TerraSAR为例，积累25景数据至少需要约3个季度的时间。因此开始有学者对小数据集时序InSAR进行相关的研究：张奎等^［10］对所有多视像元的相位逐一萃取以增强干涉对的相干性，提高了小数据集多基线InSAR技术监测密度，但由于图像数据有限，易受到低质量像素的干扰；Kang等^［11］提出时序干涉相位序列在三维上的低秩性，利用低秩张量分解算法实现基于目标的小数据集时序干涉处理，且可对三维时序干涉相位序列进行整体处理，直接剔除低质量像素，很大程度上提高了小数据集下的形变反演精度，但算法是基于单一目标的，即针对某一建筑物的形变进行研究，易受到研究区域范围大小及形变复杂度的限制，对大范围观测区域的细节形变估计精度不理想。

针对以上问题，本文深入研究了低秩张量分解算法在PSInSAR技术中的应用，并进行改进，对分解过程中的部分较小奇异值赋予更大的权值，增强其惩罚力度，防止细节形变信息的丢失，并对其改进前后的去噪效果进行实测验证。

1 算法理论

1.1 时序干涉相位理论与建模

对于重复轨道模式下的InSAR测量技术，每个干涉对的单个像素相位由多个分量组合构成^［12］

\phi ={\phi }_{\mathrm{f}\mathrm{l}\mathrm{a}\mathrm{t}}+{\phi }_{\mathrm{t}\mathrm{o}\mathrm{p}\mathrm{o}}+{\phi }_{\mathrm{d}\mathrm{e}\mathrm{f}}+{\phi }_{\mathrm{a}\mathrm{t}\mathrm{m}\mathrm{o}}+{\phi }_{\mathrm{n}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{e}}，（1）

其中，{\phi }_{\mathrm{f}\mathrm{l}\mathrm{a}\mathrm{t}}为参考椭球面引起的平地相位；{\phi }_{\mathrm{t}\mathrm{o}\mathrm{p}\mathrm{o}}为地面高程引起的地形相位；{\phi }_{\mathrm{d}\mathrm{e}\mathrm{f}}为雷达观测期间地表位移引起的形变相位；{\phi }_{\mathrm{a}\mathrm{t}\mathrm{m}\mathrm{o}}为大气状态不同而引起的大气相位；{\phi }_{\mathrm{n}\mathrm{o}\mathrm{i}\mathrm{s}\mathrm{e}}为噪声相位。

平地相位{\phi }_{\mathrm{f}\mathrm{l}\mathrm{a}\mathrm{t}}可以借助卫星精密轨道状态矢量消除^［5］，在不考虑大气和噪声相位的情况下，时序InSAR的理想去平地相位序列应由像素点的高程和雷达观测期间的形变决定。因此对于上述前提下的单像素点而言，其相位序列可以建模为

\widehat{\mathit{g}}(s,p)=\mathit{A}\mathrm{e}\mathrm{x}\mathrm{p}\left\{-\mathrm{j}\left(\frac{\mathrm{4}\mathrm{\pi }}{\lambda r}s\times \mathit{b}+\frac{\mathrm{4}\mathrm{\pi }}{\lambda }p\times \mathit{\tau }\right)\right\} .（2）

其中， A 是像素点的振幅建模，\lambda是雷达发射信号的波长，r是雷达与观测像素点之间的斜距， b 是时序干涉对之间的空间基线向量，\mathit{\tau }是时间基线向量，s和p则代表待求解的高程和形变速率参数。

而对于某一雷达观测区域，式（2）的干涉相位序列可以表示为张量的形式\mathit{𝒢}\in {\mathit{𝒜}}^{{\mathrm{I}}_{\mathrm{1}}\times {\mathrm{I}}_{\mathrm{2}}\times {\mathrm{I}}_{\mathrm{3}}}^［13］，张量的本质是高维的数据存储形式，其中I₁与I₂分别表示距离向和方位向的空间维度，I₃则表示干涉图的数量，代表时间维度。对观测区域的干涉相位序列进行张量建模^［11］公式如下

\overline{\mathit{𝒢}}(\mathit{S},\mathit{P})=\mathit{A}\mathrm{e}\mathrm{x}\mathrm{p}\left\{-\mathrm{j}\left(\frac{\mathrm{4}\mathrm{\pi }}{\lambda r}\mathit{S}\otimes \mathit{b}+\frac{\mathrm{4}\mathrm{\pi }}{\lambda }\mathit{P}\otimes \mathit{\tau }\right)\right\}，（3）

其中，\overline{\mathit{𝒢}}是观测区域的理想干涉相位张量； S 和 P 分别表示该区域内所有像素点构成的高程和形变速率矩阵，\otimes表示外积，用于维数扩展。

对于城区而言，研究区域的 S 和 P 在空间和时间上通常是具有相关性的，因此其理想去平地干涉相位张量通常具有固有的三维低秩特性^［11］。为便于理解，图1对上述张量模型进行了示例展示。

1.2 低秩张量分解算法

低秩张量分解算法的核心思想与主成分分析法^［14］一致，它相当于矩阵分解的高阶泛化，通过提取数据的主要特征分量实现对原数据的降维处理，保留其低秩主成分。式（3）所示的城区理想去平地干涉相位张量在三维结构上具有低秩性，它存在于包含噪声值的真实观测干涉相位数据中，因此可以通过张量分解算法从观测数据中提取出该相位张量的低秩成分^［15］，将其视为观测区域理想相位张量的低秩近似。此过程直接在三维结构上对时序干涉信号进行去噪处理，获取尽可能理想的干涉相位，最后通过PSInSAR技术估计形变速率值。

算法基本思想就是通过张量分解将观测相位张量分解为低秩张量部分和稀疏张量部分^［16］，求解秩的最小化

\{\widehat{\mathit{𝒳}},\widehat{\mathit{ℰ}}\}=\underset{\mathit{𝒳},\mathit{ℰ}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\mathrm{r}\mathrm{a}\mathrm{n}\mathrm{k}(\mathit{𝒳})+\gamma \Vert \mathit{ℰ}{\Vert }_{\mathrm{0}}\hspace{1em}\mathrm{s}.\mathrm{t}.\mathit{𝒳}+\mathit{ℰ}=\mathit{𝒢}，（4）

其中，\mathit{𝒢}代表去平地干涉相位张量的观测数据；\mathit{𝒳}为分解出的低秩部分，代表理想去平地干涉相位张量的低秩近似；\mathit{ℰ}为分解出的稀疏部分，代表观测数据中的噪声。

然而式（4）存在多项式复杂度的非确定性问题。为解决该问题，有学者陆续提出式（4）的替代算法，其中Kang等^［11］处理小数据集时序InSAR时采用基于目标的低秩张量分解稳健多基线干涉合成孔径雷达技术（robust multi-baseline InSAR technique via object-based low-rank tensor decomposition，RoMIO），既解决了函数的复杂度非确定性问题，又在解算过程中增强了\mathit{𝒳}的低秩性和\mathit{ℰ}的稀疏性^［17-18］，公式如下

\begin{array}{c}\{\widehat{\mathit{𝒳}},\widehat{\mathit{ℰ}}\}=\underset{\mathit{𝒳},\mathit{ℰ}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\displaystyle \sum _{n=\mathrm{1}}^{\mathrm{3}}}{‖{\mathit{w}}_{\mathrm{ℒ},n}\odot \mathit{\sigma }\left({\mathit{X}}_{(n)}\right)‖}_{\mathrm{1}}+\gamma {‖{\mathit{𝒲}}_{\mathit{ℰ}}\odot \mathit{ℰ}‖}_{\mathrm{1}},\end{array}\mathrm{s}.\mathrm{t}.\mathit{𝒳}+\mathit{ℰ}=\mathit{𝒢}.（5）

其中，{\mathit{X}}_{(n)}为张量\mathit{𝒳}的模-n展开矩阵，{\mathit{w}}_{\mathrm{ℒ},n}为矩阵{\mathit{X}}_{(n)}的奇异值权重向量，\mathit{ℰ}为稀疏向量，{\mathit{𝒲}}_{\mathit{ℰ}}是张量\mathit{E}的权重张量。优化过程中对L₁范数进行加权处理能够更好地逼近张量的秩，增强优化结果\mathit{𝒳}的低秩性。

RoMIO算法在基于目标的小数据集时序InSAR中可以取得高精度的形变反演结果，然而当研究区域范围比较大且形变比较复杂时，一些细节形变区的低秩性不够，式（5）的算法在优化过程中容易出现过度平滑的现象，导致参数反演精度较低^［11］。为尽量摆脱研究区域大小和形变复杂度的限制，增强张量分解在小数据集时序InSAR中的实用性，本文对RoMIO算法进行了双加权操作。

双加权算法针对低秩恢复张量\mathit{𝒳}进行额外约束以达到分离有效信息与噪声信息的目的。对于空间范围较大的研究区域，细节形变的相位信息在整个干涉相位图像信息中占比很小，也就是说，在上述函数的求解过程中，其对应的空间维度上的2个展开矩阵{\mathit{X}}_{(\mathrm{1})}和{\mathit{X}}_{(\mathrm{2})}的奇异值很小。二次权重向量目的在于增强此部分奇异值在优化过程中的惩罚力度及该部分的形变相位信息量，使其在迭代过程中更好地与噪声信息区分开，从而解决过度平滑问题，以防止细节损失。由奇异值分解原理可知，细节信息所占比重随奇异值成正比下降，然而对于形变估计而言细节区重要性不变，因此本文对该部分对应奇异值的惩罚力度进行正比增强，目标函数如下

\begin{array}{c}\begin{array}{l}\{\widehat{\mathit{𝒳}},\widehat{\mathit{ℰ}}\}=\underset{\mathit{𝒳},\mathit{ℰ}}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}{\displaystyle \sum _{n=\mathrm{1}}^{\mathrm{3}}}{‖{\mathit{w}}_{d,n}{\mathit{w}}_{\mathrm{ℒ},n}\odot \mathit{\sigma }\left({\mathit{X}}_{(n)}\right)‖}_{\mathrm{1}}+\\ \gamma {‖{\mathit{𝒲}}_{\mathit{ℰ}}\odot \mathit{ℰ}‖}_{\mathrm{1}},\end{array}\end{array}

\mathrm{s}.\mathrm{t}.\mathit{𝒳}+\mathit{ℰ}=\mathit{𝒢}，（6）

{\mathit{w}}_{d,\mathrm{1}}(r_{\mathrm{1}}:p_{\mathrm{1}})=\frac{{\sigma }_{r_{\mathrm{1}}}({\mathit{X}}_{(\mathrm{1})})}{\mathbf{\sigma }({\mathit{X}}_{(\mathrm{1})})(r_{\mathrm{1}}:p_{\mathrm{1}})}，

{\mathit{w}}_{d,\mathrm{2}}(r_{\mathrm{2}}:p_{\mathrm{2}})=\frac{{\sigma }_{r_{\mathrm{2}}}({\mathit{X}}_{(\mathrm{2})})}{\mathit{\sigma }({\mathit{X}}_{(\mathrm{2})})(r_{\mathrm{2}}:p_{\mathrm{2}})}，

{\mathit{w}}_{d,\mathrm{3}}=[\mathrm{1,1},\dots ,\mathrm{1}].

其中，{\mathit{w}}_{d,n}是矩阵{\mathit{X}}_{(n)}的重新加权的向量，\mathit{\sigma }({\mathit{X}}_{(\mathrm{1})})(r_{\mathrm{1}}:p_{\mathrm{1}})与\mathit{\sigma }({\mathit{X}}_{(\mathrm{2})})(r_{\mathrm{2}}:p_{\mathrm{2}})分别为细节形变区域在空间维度上对应的奇异值部分。

由于线性形变在时间上是属于高度相关的信息，细节形变相位信息在观测张量\mathit{𝒢}的模-3展开矩阵{\mathit{X}}_{(\mathrm{3})}中依然能够保留足够的低秩性，理论上不会出现过度平滑问题，因此本文设置权值向量{\mathit{w}}_{d,\mathrm{3}}=[\mathrm{1,1},\cdots ,\mathrm{1}]。根据截断核范数的低秩矩阵近似原理，最大的r个非零奇异值对矩阵的秩没有影响，作为秩函数的近似包含的信息非常少，因此为避免前r奇异值惩罚过重导致的过缩问题，选择将权值{\mathit{w}}_{d,\mathrm{1}}(\mathrm{1}:r_{\mathrm{1}})与{\mathit{w}}_{d,\mathrm{2}}(\mathrm{1}:r_{\mathrm{2}})所有元素都置1，r₁与r₂的取值参考文献［19］。而截断点p₁与p₂则没有先验知识，需借助经验进行大概的判断。本文在进行实验之前，选择多个与实验所用研究区域大小一致且细节形变较多的城区沉降干涉相位图，通过研究不同截断奇异值阈值下的低秩近似与原始图像之间的信息损失来分析细节信息所对应的奇异值占比。借此判断，p₁与p₂的值大致在190~210，实验过程中选择p₁=p₂=200，同时将权值向量剩余元素置1以防止对噪声信息的惩罚过重而无法有效滤波。

本文采用交替方向乘子法^［20-21］将目标问题分解成几个子问题进行迭代求解。下式为式（6）的增广拉格朗日函数

（7）

其中，张量\mathit{𝒴}表示拉格朗日乘数，\mu为引入的惩罚参数，求解流程如表1所示。

2 数据处理流程

观测干涉信号经张量分解去噪后，采用传统的PSInSAR进行后续处理，后文简称为PSI。PSI由Ferretti等^［6］在2001年提出，该技术首先需要识别相位稳定的永久散射体，通过永久散射体的差分相位对形变速率参数进行建模和解算，最终获取每个永久散射体（permanent scatterers，PS）点的形变速率时间序列。如前所述，具有三维低秩特性的干涉相位张量由高程和形变决定，因此在进行张量分解之前需先对观测干涉信号进行去平地处理，即去除所有干涉像对中与平地效应相关的相位成分。图2给出了基于双加权张量分解算法的小数据集时序InSAR的整个处理流程。

由于高程相位信息的低秩性高于高程残差，因此先对去平地干涉相位进行低秩张量分解再进行地形相位的去除。而大气相位在时间上是不相关的，它和噪声相位一样在三维张量结构上是不具备低秩特性的，因此原则上在对去平地观测干涉相位进行分解之后，大气与噪声相位都被提取到稀疏张量部分，高程与形变相位则被提取到低秩张量部分。最后采用单主图像的PSI技术对低秩张量部分进行时序反演操作，将外部数字高程模型（digital elevation model，DEM）投影到雷达坐标系并利用差分干涉的思想去除高程相位，进而得到形变相位分量。

3 实验分析

3.1 研究区域和数据

为验证所提算法在小数据集时序InSAR技术中的可行性，以天津市静海县（现静海区）为研究区域，将双加权RoMIO算法运用在实测数据中，并与原算法的结果进行对比。静海县位于天津西南部，属于海河冲积平原，地势平坦，地面沉降主要是由地下水的过度使用造成。对于本文的时序InSAR技术分析，文章使用时间跨度为2009年3月—2010年3月的TerraSAR-X的29景SLC，重访周期为11 d，数据空间分辨率为3 m，极化方式为HH极化，影像参数见表2。同时选择由美国“奋进号”航天飞机搭载雷达地形测绘系统获取到的分辨率为30 m的DEM作为参考高程去除地形相位。

3.2 实验结果与分析

按照PSI技术的常规流程对前述获取的29景SLC进行时序分析，首先选择2009年9月19日获取的图像作为主图像，对其余28景图像进行配准并重采样到主图像中。实验中采用振幅离差指数与相干系数结合的方式识别PS点，在小数据集下也能获取较为准确的结果。为确保实验对比结果的一致性，在每个数据集情况下应采用相同的PS分布结果，而实际上小数据集下很难获取29景大数据集下可识别的PS点，因此全文采用11景图像条件下获取的PS分布结果，如图3（a）所示。再由主SLC图像与28景从SLC图像组成时序干涉对，主从图像逐像素共轭相乘导出干涉图，并导入精准轨道数据与外部DEM数据去除平地与地形相位，从而生成差分干涉图并解算研究区域的形变相位分量。图3（b）为所得形变速率图，并以此作为后续小数据集形变速率估计的参考值。

从图（3）可见，在研究区域的中间大部分位置的形变趋势比较平滑且形变值相对较小。为展示双加权张量分解算法在复杂形变区的反演效果加成，本文主要观察图3红色方框圈出的边缘细节区的形变情况，该区域形变值较大且相对复杂。

原则上，PSI技术在29景图像的基础上足以获取较高精度的形变估计结果，而本文的研究目的是使小数据集时序InSAR的估计结果尽量达到大数据集下的精度，因此后续直接以上述形变速率值作为参考标准值，以此来衡量所提算法在小数据集下的反演精度。为使估计结果更具可比性，从所有图像数据中选取时空基线分布与大数据集大致相似的11景小数据集，图4显示了2组数据集的时空基线分布。

按照图2中的技术流程进行小数据集时序形变反演，并对图3方框内的细节区在传统PSI、RoMIO+PSI及本文所提PSI+双加权RoMIO这3种算法下的形变速率估计结果进行展示，如图5所示。

其中，图5（a）为29景影像下获取的参考值，图5（b）~5（d）分别为使用传统PSI、PSI+RoMIO、PSI+双加权RoMIO算法的11景小数据集形变反演。图中左侧红色箭头为剖面位置，从图中可以看出，对于本文所选的研究区域，在小数据集情况下，原PSI+RoMIO算法对比传统PSI并未起到明显的作用，两者在边缘的复杂形变区细节损失严重。观察红色方框标出的细节部分，传统PSI的形变速率反演结果噪点较多且相对参考值低很多；PSI+RoMIO算法的结果虽然相对光滑噪点较少，但是张量分解过程中的过度平滑使得形变速率反演值整体上出现偏移，导致估计值的低精度；而本文所提的双加权张量分解法的估计结果则明显与传统PSI的结果更加接近，噪点几乎很少，对边缘的细节形变区实现了更稳健的估计。

为更客观地验证上述结论，对以上几个估计结果在红色箭头所指的横向剖面进行了绘制，如图6所示。

很明显可以看出所提算法的形变速率剖面曲线与参考值的拟合度更高，传统PSI和原RoMIO的结果与参考值之间的误差都比较大。另外计算以上各结果细节区形变速率估计值与参考值的平均绝对误差（mean absolute error，MAE），如表3所示。从图6也可以看出，PSI+RoMIO算法由于在低秩分解过程中过度平滑的原因导致结果出现整体偏移问题，因此尽管传统PSI的结果偏离参考值的噪声散点更多，但它与参考值之间的MAE反而较小。以上数据都证明了双加权RoMIO算法确实可以提高形变反演精度，增强张量分解在时序InSAR技术上的实用性。

以上结论都验证了所提算法对细节形变区的保留效果更优，算法的二次加权增强了细节信息的权重，在成功过滤噪声的前提下又能维持细节形变区的低秩性并进行更有效恢复。为测试所提算法在不同小数据集时序InSAR场景下的性能，又分别选取了10、9、8、7、6共5种不同低数量SLC进行处理，并分别计算这5种小数据集的形变速率估计值与参考值的MAE，结果如表4所示。随着SAR图像景数越少，形变估计误差大致上呈现逐步变大的趋势。图像景数越少，时间样本不够，形变相位在时间维度上的相关性不够，导致时序干涉相位的观测张量在时间维度上低秩性会降低，估计精度也就越低。从数值上来看，在影像数量少至8景以下时，所提算法的形变速率误差超过5 mm/a，反演效果变差较为明显。

4 结语

本文通过实测验证了所提算法在小数据集时序InSAR技术中对形变反演精度的提升。对地面沉降严重的天津城区在各个算法下的形变估计结果进行了多方面的对比分析，以传统PSI技术的29景图像反演结果为参考，相比于传统PSI算法和PSI+RoMIO算法，本文所提PSI+双加权RoMIO算法的小数据集形变速率估计值在整体上与完整数据集下参考值的MAE要小很多，再结合与参考值剖面曲线的拟合度，充分验证了该算法估计小数据集形变的可靠性和高精度。

传统的时序InSAR技术需要以大量的SAR图像为前提，本文的实测结果也显示，传统算法的小数据集情况下形变估计精度非常低，很难对地面沉降做出正确的判断和预警。考虑到时序干涉相位序列的三维低秩特性，理论上利用张量分解算法可以从去平地观测干涉相位张量中提取出理想的干涉相位，将除地形和形变外的其他相位直接剔除。然而由于原优化算法的限制，实际上无法实现精准的相位提取，分析实测数据可知，通过对原算法进行双加权处理，使其在算法迭代中更好地分离形变相位信息与噪声相位信息，实现了更有效的低秩分解，在一定范围内可以取得更高精度的形变结果。因此本文所提算法对城区的高分辨率形变反演研究具有一定价值，能够在保证形变估计高精度的同时，节约时间和数据资源，使得PSI技术的实时性更高。

基于张量分解算法对于数据低秩性的要求，在保证高精度的前提下，所提算法能相对提高对大范围研究区域的适用性，使其不再局限于单个建筑物的形变反演。然而，对于更大的场景仍存在过度平滑的问题，且运算效率会相对下降，但在同等场景下，其性能依然优于原算法。因此，今后还需要就此问题继续进行更深入的小数据集时序InSAR的研究。

参考文献

[1] 侯俊东， 肖人彬. 地质灾害防治项目治理的社会网络结构优化及策略［J］. 管理评论， 2015， 27（2）： 13-25. DOI：10.14120/j.cnki.cn11-5057/f.2015.02.002 .

[2] 肖宁宁. 基于时间序列InSAR技术的成都地区地表沉降研究［D］. 成都： 电子科技大学， 2019.

[3] 陈洋， 陶秋香， 刘国林， 等. InSAR与概率积分法联合的矿区地表沉降精细化监测方法［J］. 地球物理学报， 2021， 64（10）： 3554-3566. DOI： 10.6038/ cjg2021O0421 .

[4] 汪钰晴， 唐伶俐， 王新鸿， 等.基于时序InSAR技术的中巴走廊盖孜河谷段地表形变分析［J］.中国科学院大学学报， 2022， 39（3）： 343-351. DOI： 10.7523/j. ucas.2020.0046 .

[5] 朱建军， 杨泽发， 李志伟. InSAR矿区地表三维形变监测与预计研究进展［J］. 测绘学报， 2019， 48（2）： 135-144.

[6] Ferretti A， Prati C， Rocca F. Permanent scatterers in SAR interferometry［J］. IEEE Transactions on Geoscience and Remote Sensing， 2001， 39（1）： 8-20. DOI： 10.1109/36.898661 .

[7] Berardino P， Fornaro G， Lanari R， et al. A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms［J］. IEEE Transactions on Geoscience and Remote Sensing， 2002， 40（11）： 2375-2383. DOI： 10.1109/TGRS. 2002. 803792 .

[8] 刘星， 吕孝雷. 基于卡尔曼滤波的PS-InSAR地表形变预测方法［J］. 中国科学院大学学报， 2017， 34（6）：743-750. DOI： 10.7523/j.issn.2095-6134.2017. 06.011 .

[9] Ma P F， Wang W X， Zhang B W， et al. Remotely sensing large- and small-scale ground subsidence： a case study of the Guangdong–Hong Kong–Macao Greater Bay Area of China［J］. Remote Sensing of Environment， 2019， 232： 111282. DOI： 10.1016/j.rse. 2019.111282 .

[10] 张奎， 汪慧， 宋瑞庆， 等. 相位萃取技术在InSAR小数据集上的应用［J］. 遥感学报， 2017， 21（4）： 645-652. DOI： 10.11834/jrs.20175337 .

[11] Kang J， Wang Y Y， Schmitt M， et al. Object-based multipass InSAR via robust low rank tensor decomposition［J］. IEEE Transactions on Geoscience and Remote Sensing， 2018， 56（6）： 3062-3077. DOI： 10.1109/TGRS. 2018.2790480 .

[12] Hanssen R F. Radar Interferometry： data interpretation and error analysis［M］. Dordrecht： Springer Netherlands， 2001. DOI： 10.1007/0-306-47633-9 .

[13] Kolda T G， Bader B W. Tensor decompositions and applications［J］. SIAM Review， 2009， 51（3）： 455–500. DOI： 10.1137/07070111X .

[14] Hotelling H. Analysis of a complex of statistical variables into principal components［J］. Journal of Educational Psychology， 1933， 24（7）： 498-520. DOI： 10.1037/h0070888 .

[15] Kang Z， Peng C， Cheng Q. Robust PCA via nonconvex rank approximation［C］//2015 IEEE International Conference on Data Mining. Atlantic City， NJ， USA. IEEE， 2015： 211-220. DOI： 10.1109/ICDM.2015.15 .

[16] Goldfarb D， Qin Z W. Robust low-rank tensor recovery： models and algorithms［J］. SIAM Journal on Matrix Analysis and Applications， 2014， 35（1）： 225-253. DOI： 10.1137/130905010 .

[17] Candès E J， Wakin M B， Boyd S P. Enhancing sparsity by reweighted ℓ1 minimization［J］. Journal of Fourier Analysis and Applications， 2008， 14（5）： 877-905. DOI： 10.1007/s00041-008-9045-x .

[18] Peng Y G， Suo J L， Dai Q H， et al. Reweighted low-rank matrix recovery and its application in image restoration［J］. IEEE Transactions on Cybernetics， 2014， 44（12）： 2418-2430. DOI： 10.1109/TCYB.2014. 2307854 .

[19] Hu Y， Zhang D B， Ye J P， et al. Fast and accurate matrix completion via truncated nuclear norm regularization［J］. IEEE Transactions on Pattern Analysis and Machine Intelligence， 2013， 35（9）： 2117-2130. DOI： 10.1109/TPAMI.2012.271 .

[20] Boyd S. Distributed optimization and statistical learning via the alternating direction method of multipliers［J］. Foundations and Trends® in Machine Learning， 2010， 3（1）： 1-122. DOI： 10.1561/ 2200000016 .

[21] Parikh N. Proximal algorithms［J］. Foundations and Trends® in Optimization， 2014， 1（3）： 127-239. DOI： 10.1561/240000000 .

基金资助

国家自然科学基金(62001451)