CN116859391A - 基于高频纹理互相关的InSAR对流层延迟矫正方法 - Google Patents

Abstract

本发明公开一种基于高频纹理互相关的InSAR对流层延迟矫正方法，涉及雷达测量领域，通过对干涉相位和地形高程进行高通滤波得到的高频纹理图之间的相关系数实现。具体地，本发明由低分辨率分层延迟估计模块以及高分辨率分层延迟矫正模块实现。低分辨率分层延迟估计模块在空间域对每个干涉图逐窗口执行，以估计出低分辨率的相位‑高程斜率图和截距图。其中，斜率的估计用到高频纹理相关性。高分辨率分层延迟矫正模块为解决窗口内的对流层延迟异质性而设计，其在时间域和空间域迭代执行以对残余对流层延迟进行矫正。本发明可以有效去除干涉相位中的对流层延迟，恢复出正确的形变信息，提升InSAR形变监测的精度。

Description

基于高频纹理互相关的InSAR对流层延迟矫正方法

技术领域

本发明属于雷达测量领域，具体涉及一种基于高频纹理互相关的InSAR（Interferometric Synthetic Aperture Radar, 干涉合成孔径雷达）对流层延迟矫正方法。

背景技术

得益于高分辨率、广域覆盖、全天时全天候测量，InSAR成为地表形变监测的一种强有力手段。然而InSAR基于相位信息的测量方式使其精度受限于SAR信号传播路径上各种误差的干扰。其中，对流层延迟误差是长期以来限制InSAR信息正确提取和解译的主要干扰来源之一。对流层延迟的复杂性在于其关联因素众多（如，温度、湿度、压强、风速等），且每种关联因素自身在时间域和空间域上特性复杂。

幸运的是，在统计学上按照系统性和随机性，对流层延迟在时间域和空间域分别可以被划分为不同的分量，可以针对各个分量的不同表现形式进行针对性地处理。在空间域，对流层延迟可以分为垂直分层分量（系统）以及紊流分量（随机）。在时间域，对流层延迟可以分为季节性分量（系统）和非季节性分量（随机）。其中，垂直分层分量与地形起伏直接相关，且表现出季节相关性，而紊流分量则表现出时间上的随机性。研究表明仅一天时间间隔，紊流分量即可失去相关性。通常，随机分量在一定程度上可以通过滤波进行缓解，因此矫正系统性的垂直分层分量成为对流层延迟矫正的首要任务。

矫正对流层分层延迟的方法主要有两大类，一类是基于外部数据进行矫正，另一类是基于InSAR相位自身信息进行矫正。尽管基于外部数据的矫正方法由于融合外部独立测量数据显得更为可靠，但当前所能实现的外部数据产品的时间、空间分辨率对于应用到现代高分辨率InSAR信息提取上还略有不足。而基于InSAR相位自身信息进行矫正的方法，一般利用相位和高程的关系对对流层分层延迟进行估计并去除。但该方法存在两大挑战，一是对流层延迟具有空间异质性，二是地形相关的形变信号会对估计造成影响。异质性要求在小尺度上进行相位-高程关系的估计，然而在小尺度上，地形相关的形变信号对估计的影响将会变得显著。

针对对流层分层延迟异质性与地形相关的形变相互制约这一问题，以往的研究要么将窗口设计的较大通过牺牲估计的精度（由忽略部分对流层延迟异质性导致）以削弱地形相关形变的影响，要么根据经验关系对分层延迟以及形变信号等相位成分进行建模，联合求解各个模型参数。然而，根据经验关系建模的方法对于复杂形变情况可能并不适用，且联合求解的计算量和计算复杂度将大大增加。因此亟待一种能有效解耦对流层分层延迟的异质性与地形相关的形变的相互制约关系，又不牺牲估计精度以及计算简便性的方法。

发明内容

为解决上述技术问题，本发明提出一种基于高频纹理互相关的InSAR对流层延迟矫正方法。该方法所依赖的基本原理为：在空间高频带，形变信号、紊流等相位分量与地形相关性变弱，而对流层分层延迟由于直接与地形起伏相关，因此依然保持与地形相对较高的相关性。这样，通过减小干涉相位和地形的高频纹理相关性，即可实现不受形变影响的对流层分层延迟估计。

为达到上述目的，本发明的技术方案为：

一种基于高频纹理互相关的InSAR对流层延迟矫正方法，包括如下步骤：

步骤1、对原始干涉相位时间序列进行估计窗口划分；

步骤2、进行基于高频纹理互相关的相位-高程斜率估计，获得空间离散分布的斜率估计值；

步骤3、对步骤2所得的空间离散分布的斜率估计值进行低通滤波；

步骤4、对步骤3所得的低通滤波后的空间离散分布的斜率估计值进行空间插值；

步骤5、对去除由步骤4所得的斜率估计值乘以地形高程后的相位进行低通滤波以估计截距；

步骤6、基于步骤4所得的斜率估计值与步骤5估计的截距对原始干涉相位进行低分辨率分层延迟去除；

步骤7、对步骤6所得的低分辨率分层延迟去除后的干涉相位进行高分辨率分层延迟矫正。

本发明提供的基于高频纹理互相关的InSAR对流层延迟矫正系统，包括低分辨率分层延迟估计模块和高分辨率分层延迟矫正模块，其中，所述低分辨率分层延迟估计模块实现步骤1-步骤6，所述高分辨率分层延迟矫正模块实现步骤7。

有益效果：

本发明所提方法可以有效解耦对流层分层延迟异质性与地形相关的形变相互制约关系，并实现高精度的对流层分层延迟矫正，从而恢复出正确的形变信息，进而提升合成孔径雷达干涉形变监测的精度，为地表形变机制解译提供更准确的证据。此外，方法计算复杂度低，稳健性、普适性强，且由于不依赖于外部大气测量数据，可以进行实时处理。

附图说明

图1为本发明的一种基于高频纹理互相关的InSAR对流层延迟矫正方法的流程图；

图2为对流层延迟矫正前后的干涉相位（主影像获取日期为2020-02-02，参考影像获取日期为2015-02-23）的相位-高程散点图；其中

（a）为未矫正的干涉相位的相位-高程散点图；

（b）为利用外部数据ERA5矫正对流层延迟后的干涉相位的相位-高程散点图；

（c）为采用全局线性拟合矫正对流层延迟后的干涉相位的相位-高程散点图；

（d）为采用局部线性拟合矫正对流延迟后的干涉相位的相位-高程散点图；

（e）为本发明所提方法矫正对流层延迟后的干涉相位的相位-高程散点图。

图3为残余干涉相位的空间域标准偏差结果图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整的描述，显然，所描述的实施例仅为本发明的一部分实施例，而不是全部的实施例，基于本发明中的实施例，本领域的普通技术人员在不付出创造性劳动的前提下所获得的所有其他实施例，都属于本发明的保护范围。

结合本发明的一个实施例，提出一种基于高频纹理互相关的InSAR对流层延迟矫正方法，如图1所示，包括如下步骤：

步骤101：对原始干涉相位时间序列进行估计窗口划分：

对原始干涉相位时间序列中每一个干涉相位（后简称为原始干涉相位）进行估计窗口划分。估计窗口的大小（以方形窗为例，记边长为）的选取主要在两个方面进行折中。一方面，对流层延迟的空间异质性要求估计窗口尽可能小。另一方面，估计所需要的丰富纹理信息要求估计窗口尽可能大。估计步长/>（假设步长在距离向方位向等距），即表征斜率估计值的分布密度的参数，其选取也需要在两个方面进行折中。一方面，低通滤波（步骤105）为避免破坏形变信号要求斜率估计值在二维空间中的分布尽可能稀疏，即估计步长要尽可能大。而另一方面，对斜率进行空间插值（步骤104）时为避免数据过度外推又约束估计步长不能太大。综合起见，本发明优选估计步长/>设置为最大形变范围。确定完估计步长/>后，估计窗口的边长/>的设置满足窗口重叠率在30%~50%之间即可，这样的设置能在窗口大小的两个要求之间达到一个折中。

步骤102：进行基于高频纹理互相关的相位-高程斜率估计：

估计窗口划分完毕之后，接下来对原始干涉相位在每一个估计窗口内，基于干涉相位和地形高程的高频纹理相关性估计斜率。以第（/>，/>为时间序列的个数）个时间序列的第/>（/>，/>为估计窗口的行数，/>为估计窗口的列数）个估计窗口为例，为了更快的搜索到解，首先需要根据估计窗口内的相位-高程线性关系计算出斜率的初值，即：

(1)

其中，表示地形高程，/>表示原始干涉相位时间序列，右下标“/>”表示第/>个时间序列，/>为原始干涉相位。右上标“/>”表示由第/>个估计窗口限制的子矩阵，上标“/>”表示将矩阵转换为列向量的操作，“/>”表示求矩阵伪逆的操作，“/>”表示求矩阵转置的操作，/>表示元素全为1的矩阵，矩阵大小与/>相同，/>和/>分别为在第/>个时间序列，第/>个窗口内根据相位-高程线性关系估计的斜率和截距的初值，其中，截距初值在下文不被利用。

得到斜率的初值后，在/>的基础上进行微调，直到/>的高频分量与地形高程/>的高频分量的相关性达到最小时，输出此时的/>作为第/>个时间序列，第个窗口内估计得到的斜率，用公式描述为：

(2)

其中，表示高通滤波操作，作用于矩阵的符号/>表示计算矩阵转换为列向量后的模长。记第/>个窗口的中心坐标为/>，且记第/>个时间序列估计的斜率为/>，则把/>中坐标/>为/>的元素赋值为/>，即，/>。

对每一个时间序列的原始干涉相位逐窗口执行步骤102，运行完毕后将得到个维空间离散分布的斜率估计值矩阵。

步骤103：对步骤102所得的空间离散分布的斜率估计值进行低通滤波：

对步骤102所得的维空间离散分布的斜率估计值进行空间低通滤波。低通滤波采用高斯滤波，一方面，低通滤波可以为低分辨率分层延迟估计模块所得到的斜率图与截距图分辨率相匹配做准备，另一方面，低通滤波可以有效避免一些局部斜率估计误差（表现为空间高频分量）造成的影响。

步骤104：对步骤103所得的空间离散分布的斜率估计值进行空间插值：

对步骤103所得维空间离散分布的斜率估计值进行二维空间插值。插值可以选用任一常用的二维空间插值方法，但需要注意的是，插值后得到的完整的估计的斜率/>的数据尺寸应与原始干涉相位/>的数据尺寸大小相同。本发明推荐使用三次样条插值，其可以实现在数据边缘进行适当的数据外推。

步骤105：对去除由步骤4所得的完整的估计的斜率乘以地形高程后的相位进行低通滤波以估计截距；

用步骤104得到的完整的估计的斜率乘以地形高程/>并从原始干涉相位中减去，此时得到的相位/>还需要进一步解决截距的影响。这里，本发明直接采用低通滤波来估计截距。记估计的截距为/>，则步骤105所表示的操作可以用数学语言描述为：

(3)

其中，表示低通滤波操作。需要注意的是，低通滤波的参数需要选取的比较大以避免破坏形变信号。具体参数的选择需要视研究区域的形变分布范围而定，根据对研究区域的形变类型的先验知识即可进行判断。

步骤106：对原始干涉相位进行低分辨率垂直分层延迟去除：

对于第个时间序列，通过步骤101-步骤105已得到完整的估计的斜率/>和估计的截距/>，则估计的对流层分层延迟为/>。将其从原始干涉相位/>中减去，得到初步矫正后的干涉相位。此时大部分对流层分层延迟已经去除，相位以形变相位为主导，因此记作/>，则有：

(4)

其中，为形变相位为主导的相位。

步骤107：进行高分辨率垂直分层延迟矫正：

在空间上任意一点从时间域来看，原始相位时间序列为，其中下标/>表示第/>个时间序列，N为时间序列的个数；且对于/>有/>，即/>为/>中空间点/>处的值。类似地有定义：空间点/>处真实和估计的矫正对流层分层延迟相位后的相位时间序列分别为和/>，且，/>；空间点/>处真实的斜率时间序列和估计的斜率时间序列分别为/>和，且/>，/>；空间点/>处真实的截距时间序列分别和估计的截距时间序列分别为和/>，且，/>；空间点/>处的地形高程为/>。那么有：

(5)

(6)

由于截距是根据空间低通滤波获取的，因此截距误差与斜率误差直接相关，将斜率误差记作/>，因此有：

(7)

其中，表征了/>中与/>相关的信号分量。为了恢复出/>，需要估计出和/>。

整体性的误差一般来自空间异质性，由于本发明的窗口还是不够小，即使在单个窗口内空间异质性也会存在。不过需要强调的是，垂直分层延迟在空间具有一定的相关性，从而/>与/>具有相关性。这样，式（7）又可写作：

(8)

其中，为随着/>到/>的变化而变化的系数。

为简单起见，以线性信号模型近似/>，其中/>和/>分别为线性模型速率和常数项，/>表征SAR数据采集时间，/>表示/>在第/>个时间序列的值，则/>可以通过最小二乘估计出来：

(9)

注意到，这里还有一个潜在的假设，即与/>不相关，不符合这一条件的形变则会受到影响。为此，本发明设计了一个准则：相邻SAR采样间隔的形变速率之差的标准偏差，用以表示速率变化的紊乱性，即：

(10)

其中，

(11)

(12)

其中，和/>为中间参数，/>表示任一一维时间序列，/>表示/>在第/>个时间序列的值。接下来，将/>的改善作为触发条件用来控制是否对/>进行更新：如果/>相比/>能够满足触发条件，即，那么就进行如下更新：

(13)

(14)

其中，矩阵记录了每个像素对触发条件的满足性，满足则对应像素位置的元素值为1，否则为0，/>为/>中空间位置为/>处的值。公式（9），公式（13）以及公式（14）所代表的操作可以迭代进行以提高估计的精度。迭代的终止条件设置为：当有90%以上的像素满足/>，则迭代终止，/>表示所设置的阈值。

下面结合具体实施例对本发明的技术方案作进一步详细说明。

实施例1

数据使用2015-02-09至2017-09-18期间采集的28个L波段ALOS-2降轨（轨道号为23）InSAR时间序列，并设置2015年2月23日获取的SAR影像为参考影像。该数据集覆盖了约的区域，空间分辨率为0.000185185°。

DEM误差引起的地形地形相位误差以及可能的电离层延迟误差在常规的InSAR预处理过程中去除。并将此时得到的相位时间序列称为原始相位时间序列，在此基础上应用本发明所提出的方法进行对流层延迟矫正。为了更好的体现本发明的矫正效果，基于外部独立数据ERA5、基于全局线性拟合以及基于局部线性拟合的方法用来作为对比。在实施案例中，对于本发明所提出基于高频纹理互相关的对流层延迟矫正方法，我们将斜率估计窗口大小设置为141×141个像素，相邻窗口之间的重叠为40%，使用高斯滤波器得到高频纹理，滤波器大小为13像素。对于用于对比的局部线性方法，使用相同的窗口相关参数，去除估计斜率和截距的高频成分并将斜率和截距插值到整个二维空间。

首先以一张干涉图为例（获取日期为2015年2月23日和2017年8月21日），在图2中以相位-高程散点图的形式显示了本发明所提方法和对比方法的矫正结果。其中，基于ERA5的矫正结果在一定程度上可以作为参考。从图2的（a）可以看出，原始干涉相位的相位-高程散点图呈现出很明显的线性趋势，表明此时相位与地形高程存在较强的相关性，即对流层延迟的影响比较严重。在非形变区域计算相位的空间域标准偏差为2.0cm。经ERA5矫正后（图2的（b）），相位和地形高程的相关性大大减弱，表明大部分分层延迟已被去除，但较大的空间域标准偏差（1.3cm）说明干涉相位中还存在局部未校正的对流层延迟。得益于原始干涉相位与高程的关系近乎线性，全局线性拟合方法（图2的（c））也实现了与ERA5接近的效果，但全局拟合不能解决对流层延迟异质性问题，因此相位的标准偏差依然较大（1.2cm）。局部线性拟合方法（图2的（d））可以更好地解决异质性问题，从而提供比全局线性方法更精细的结果，实现了相位的空间域标准偏差进一步减小（1.1cm）。然而局部线性拟合会使地形相关的形变信号被平滑。相比较而言，本发明所提方法（图2的（e））不仅很大程度上去除了残余分层延迟，实现了相位的空间域标准偏差最小（0.8cm），还能有效保留两座活火山的形变信号。

接下来，从统计意义上展示本发明所提方法和对比方法的性能。图3 给出了所有残余相位时间序列的空间域标准偏差，其中形变区域被赋予掩膜。经ERA5校正后，空间域标准偏差的中位数从1.31厘米减少到0.97厘米，而全局线性拟合、局部线性拟合和纹理相关的方法分别实现了30%、44%和51%的空间域标准偏差降低，凸出了局部窗口估计的优越性，特别是本发明所提方法，它使空间域标准偏差的中位数下降到0.64厘米。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，且应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (7)

1.一种基于高频纹理互相关的InSAR对流层延迟矫正方法，其特征在于，包括如下步骤：

步骤1、对原始干涉相位时间序列进行估计窗口划分；

步骤2、进行基于高频纹理互相关的相位-高程斜率估计，获得空间离散分布的斜率估计值；

步骤3、对步骤2所得的空间离散分布的斜率估计值进行低通滤波；

步骤4、对步骤3所得的低通滤波后的空间离散分布的斜率估计值进行空间插值；

步骤5、对去除由步骤4所得的斜率估计值乘以地形高程后的相位进行低通滤波以估计截距；

步骤6、基于步骤4所得的斜率估计值与步骤5估计的截距对原始干涉相位进行低分辨率分层延迟去除；

步骤7、对步骤6所得的低分辨率分层延迟去除后的干涉相位进行高分辨率分层延迟矫正。

2.根据权利要求1所述的一种基于高频纹理互相关的InSAR对流层延迟矫正方法，其特征在于：所述步骤1包括：对原始干涉相位时间序列中每一个干涉相位，即原始干涉相位进行估计窗口划分；其中估计窗口的大小同时考虑估计所需的地形纹理信息的丰富程度和对流层延迟的空间异质性；估计步长的选择同时考虑与截距滤波分辨率的匹配和避免插值时数据过度外推。

3.根据权利要求1所述的一种基于高频纹理互相关的InSAR对流层延迟矫正方法，其特征在于：所述步骤2包括：对原始干涉相位逐估计窗口进行基于高频纹理互相关的相位-高程斜率估计；首先对估计窗口内原始干涉相位和地形高程做线性拟合，拟合所得斜率作为初值；然后在初值的基础上逐步调整斜率值并乘以地形高程后从原始干涉相位中减去，得到新的干涉相位，评估新的干涉相位和地形高程的高频纹理相关性，即高频纹理图的相关系数，输出使得高频纹理相关性达到最小的斜率作为该估计窗口的斜率估计值；逐估计窗口执行所述步骤2后，得到空间离散分布的斜率估计值。

4.根据权利要求3所述的一种基于高频纹理互相关的InSAR对流层延迟矫正方法，其特征在于：所述步骤4中，采用二维数据插值方法，将空间离散分布的斜率估计值插值到整个二维空间中，使得斜率的数据大小和原始干涉相位的数据大小相同。

5.根据权利要求4所述的一种基于高频纹理互相关的InSAR对流层延迟矫正方法，其特征在于：所述步骤5包括：用步骤4所得的斜率估计值乘以地形高程后，从原始干涉相位中减去，对剩余的相位做空间低通滤波得到截距的估计值；其中空间低通滤波的参数的选取不能破坏形变信号。

6.根据权利要求5所述的一种基于高频纹理互相关的InSAR对流层延迟矫正方法，其特征在于：所述步骤6包括：估计的低分辨率分层延迟由步骤4得到的斜率估计值乘以地形高程并加上步骤5所得到的估计的截距生成，将其从原始干涉相位中减去得到低分辨率分层延迟去除后的干涉相位。

7.根据权利要求6所述的一种基于高频纹理互相关的InSAR对流层延迟矫正方法，其特征在于：所述步骤7中，所述高分辨率分层延迟矫正在时间域和空间域迭代执行，包括：

步骤（1）在时间域基于分层延迟空间相关的特性，逐像素求解残余斜率估计误差；

步骤（2）在空间域进行低通滤波以估计残余截距；

步骤（3）对满足相邻速率之差的标准偏差减小的像素迭代执行步骤（1）、步骤（2），直到相邻速率之差的标准偏差收敛。