CN116338607A

CN116338607A - 时间域和空间域两步式InSAR对流层延迟矫正方法

Info

Publication number: CN116338607A
Application number: CN202310565150.6A
Authority: CN
Inventors: 杨晴月; 李俊峰; 陆萍萍; 蔡永华; 王宇
Original assignee: Aerospace Information Research Institute of CAS
Current assignee: Aerospace Information Research Institute of CAS
Priority date: 2023-05-19
Filing date: 2023-05-19
Publication date: 2023-06-27
Anticipated expiration: 2043-05-19
Also published as: CN116338607B

Abstract

本发明公开一种时间域和空间域两步式InSAR对流层延迟矫正方法，涉及雷达测量领域。该方法包括时间域矫正和空间域矫正两大主要步骤。时间域矫正联合求解形变速率和季节调制的大气对流层延迟相位模型参数。根据时间域矫正后的形变速率图识别形变区域，进而对每个形变区域进行针对性地空间域矫正。空间域矫正基于相位‑高程线性关系矫正，估计窗口的大小根据相关性峰值自适应确定。经时间、空间域矫正后即可得到每一形变区域的形变速率和形变时间序列。本发明可以有效提高广域微弱形变识别的准确度以及形变时间序列监测的精度，且由于不依赖外部大气数据，可以应用于实时时间序列干涉处理。

Description

时间域和空间域两步式InSAR对流层延迟矫正方法

技术领域

本发明属于雷达测量领域，具体涉及一种时间域和空间域两步式InSAR（Interferometric Synthetic Aperture Radar）对流层延迟矫正方法，可应用于时间序列干涉处理中以提高形变速率和形变时间序列的测量精度。

背景技术

大气对流层延迟（以下简称为大气延迟）相位由雷达信号经过大气对流层引起。易变的大气条件（温度、压强和水汽等）使得雷达信号发生不同程度的相位延迟。一方面，使雷达信号在传播时间上发生延迟，另一方面，使雷达信号的空间传播路径发生弯曲，从而导致波的实际传播路径比理论路径长。其中信号传播路径弯曲的影响可以忽略不计。就重轨干涉测量来说，由于两次SAR成像时间及视角不同，电磁波传播所经历的大气状态也不同，这种影响在干涉处理过程中不能被对消，因而在干涉相位中引入了相位误差。这种额外的相位成分会给干涉图的正确解译和信息提取带来困难，是干涉测量精度的主要限制因素之一。

干涉相位中的大气延迟在空间上由三个分量组成：短空间尺度分量（几公里），由对流层的紊流或相干动力过程引入；长空间尺度分量（几十公里），由压强、温度及湿度的横向变化引入；垂直分层分量，由压强、温度及相对湿度随海拔高程的变化引入。在时间上，大气延迟表现为受季节调制和与季节变化无关的两个部分。其中，受季节调制的部分主要与空间上垂直分层分量相关。鉴于此，时间序列InSAR算法通常所采用的空时滤波去除大气延迟影响的方法，由于其所基于的大气延迟服从高斯分布的随机性假设不能被很好的满足，必然会引入形变估计误差，进而导致形变区域的误判和漏判。

为了准确识别形变区域，需要去除具有时间相关性的垂直分层延迟的影响。理论上，干涉相位中的垂直分层延迟可以沿高程对折射率进行积分得到。然而，要获得干涉相位中每个像素的垂直折射率剖线是不现实的。因而在实际应用中，垂直分层延迟一般基于外部数据进行矫正或基于高程和相位的经验关系建模而得。

基于外部数据的大气对流层延迟矫正方法根据所使用的外部数据的不同分为四类。第一类是基于地面气象观测信息建模的校正方法，利用地面气象观测数据（温度、湿度、气压等）估计对流层分层延迟；第二类是基于GPS数据的校正方法，反演SAR数据覆盖范围内各GPS站点的天顶对流层延迟，并对GPS网络所覆盖区域的天顶对流层延迟进行插值，进而获得完整的InSAR对流层延迟；第三类是基于空间辐射计测量的水汽数据的矫正方法，利用MODIS、MERIS等近红外水汽产品计算出的延迟量进行对流层湿延迟矫正；第四类是基于数值大气模型的对流层延迟校正法，利用数值大气模型计算出SAR数据获取时刻的大气参量，对干涉相位进行对流层延迟矫正。

然而，目前基于外部数据的方法的应用效果通常受到外部数据的空间、时间分辨率以及数据的精度限制。尽管地基的测量结果相当精确，地面气象站和GPS站点一般布设比较稀疏，单纯利用地面气象信息建模或GPS数据进行对流层延迟校正的方法并没有太大的优势。MERIS和MODIS数据分辨率较高，但只能在白天无云的条件下工作，对于夜间采集的SAR数据，无法利用MERIS和MODIS的水汽产品进行对流层湿延迟校正，对于有云区则只能得到云上部分的水汽含量。数值大气模型虽不受云和日光的限制，但其空间分辨率较低，且模型计算的对流层延迟精度严重依赖模型的初值和边界条件，稳定性相对较差。另外外部数据的获取都有不同程度的时间延迟，不能做到实时处理。

基于相位的分层延迟矫正方法能有效克服外部数据空间、时间分辨率低和精度低的缺点。相位基的方法可靠性主要受形变区域以及大气对流层延迟的空间异质性影响。一方面，大气分层延迟受季节调制呈现时间相关性，干扰形变区域的准确识别，而形变区域又会影响大气分层延迟的估计。另一方面，当大气存在异质性的情况下，单一的经验模型可能不适用于整个区域。此外，现有模型都是经验模型（线性模型、指数模型或幂律模型），在特定地区和情况下能提供好的结果，但不一定适用于所有地区。

综合以上分析，基于外部数据的大气对流层延迟方法由于引入了外部测量数据或专业的数值模型，在一些情况下能达到很好的效果，但在其它情况下依然需要基于相位的对流层延迟矫正方法作为补充。因此亟待提出一种高精度，稳健性、普适性强的基于相位的大气对流层延迟矫正方法。

发明内容

为解决上述技术问题，本发明提出一种时间域和空间域两步式InSAR对流层延迟矫正方法，其为一种基于相位的大气对流层延迟矫正方法，以提高时间序列干涉处理所得形变速率和形变时间序列的精度。本发明的方法精度高，稳健性、普适性强，且由于不依赖于外部大气测量数据，可以进行实时处理。

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

一种时间域和空间域两步式InSAR对流层延迟矫正方法，包括如下步骤：

步骤1、预处理及建立干涉相位信号模型；

步骤2、逐点分析联合求解形变速率和季节性对流层延迟模型参数；

步骤3、季节性对流层延迟模型参数空间一致性修正；

步骤4、重新逐点分析联合求解形变速率和季节性对流层延迟模型参数；

步骤5、形变区域识别；

步骤6、季节性对流层延迟去除；

步骤7、空间域矫正。

进一步地，所述步骤1包括：预处理步骤包括干涉组合，根据高相干面积比筛选干涉对，解缠干涉相位，矫正解缠误差以及DEM误差，去趋势性相位；经过预处理后，干涉相位中剩余相位成分包含大气延迟和形变相位，根据随时间变化情况对二者分别在时间域建模；接下来用建模的信号构建干涉观测方程。

进一步地，所述步骤2包括：对干涉图进行分块，每个干涉图块内各自设置参考点，在每一干涉图块内，逐点求解干涉观测方程得到形变速率和季节性对流层延迟模型参数。

进一步地，所述步骤3包括：根据参数比众数确定每一干涉图块内季节性对流层延迟模型参数空间一致的参数比关系。

进一步地，所述步骤4包括：将空间一致性的参数比关系带入干涉观测方程，重新在每一干涉图块内逐点分析求解观测方程得到形变速率和季节性对流层延迟模型参数。

进一步地，所述步骤5包括：根据步骤2和步骤4得到的形变速率结果以及拟合方差分别识别具有线性形变模式的形变区域以及具有非线性形变模式的形变区域。

进一步地，所述步骤6包括：对具有非线性形变模式的形变区域的季节性对流层延迟模型参数依据空间相关性信息进行插值修正，进而得到建模的季节性对流层延迟，并从反演的相位时间序列中去除。

进一步地，所述步骤7包括：在每个形变区域根据高程与相位的相关性峰值自适应确定窗口大小，并根据窗口内的高程-相位线性模型在逐个相位时间序列中去除残余对流层延迟相位。

有益效果：

本发明所提算法可以提高时间序列干涉处理所得形变速率和形变时间序列的精度。方法精度高，稳健性、普适性强，且由于不依赖于外部大气测量数据，可以进行实时处理。

附图说明

图1为本发明的时间域和空间域两步式InSAR对流层延迟矫正方法流程图；

图2为形变速率（绝对值）图；其中，(a)为不对大气延迟进行处理得到的形变速率，(b)为用高程—相位线性模型以及时空滤波处理大气延迟后得到的形变速率图，(c)为对ERA5数据进行精确插值获取大气延迟并矫正得到的形变速率图，(d)为利用本发明提出的方法对大气延迟矫正后得到的形变速率图。

图3为选取的三个典型形变时间序列图；其中，(a)为一个位于线性形变区域的像素点的形变时间序列图，(b)为一个位于非线性形变区域的像素点的形变时间序列图，(c)为一个位于稳定区域的像素点的形变时间序列图。所有图中分别展示了不做处理得到的原始形变时间序列、时间域矫正后的形变时间序列、空间域矫正后的形变时间序列，即本发明所展示方法得到的最终形变时间序列，以及时间域估计的形变速率。

具体实施方式

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

结合本发明的一个实施例，提出一种时间域空间域两步式InSAR对流层延迟矫正方法，如图1所示，包括如下步骤：

步骤101：预处理及建立干涉相位信号模型：

首先进行干涉组合，根据高相干面积比筛选干涉对。解缠干涉相位，矫正必要的解缠误差等处理误差和DEM误差，并去趋势性相位。

接下来对干涉相位信号建模。由于预处理中已经对解缠误差等处理误差、DEM误差、趋势性相位进行了矫正，因此接下来的分析中不再强调这几个相位项。首先对大气对流层延迟中的季节性变化部分建模。与以往的研究不同，对季节性对流层延迟的空间特性做弱约束，允许其充分的空间异质性，而在时间上做强约束。假设季节性对流层延迟时间调制模式在空间上一致，不同的位置仅存在幅度的不同，则有：

（1）

其中，

表示大气对流层延迟中被建模的部分，/>

表示空间坐标，/>

表示时间，单位为天，/>

用以表示季节性大气对流层延迟在空间上的变化，/>

用以表示季节性大气对流层延迟在时间上的变化。考虑到以一年一周期、半年一周期、四个月一周期以及三个月一周期的季节性调制，/>

可以建模为：

（2）

其中，

为一年一周期所对应的频率,/>

为一年包含的天数，可以取值为365.25天。/>

，/>

，/>

为待求解的季节性对流层延迟模型参数。

将对流层延迟的剩余分量表示为

，则其中包含季节性对流层延迟建模不完备的分量以及大气紊流分量。由于/>

相当于季节性对流层延迟的低频分量，剩余高频分量经SAR获取采样后表现为噪声特性，不会对形变的估计造成较大影响。至于大气紊流分量，研究表明，相隔一天时间大气紊流就会失去相关性，可被视作时间随机分量。完整的对流层延迟相位可以表示为：

（3）

接下来对形变信号进行建模。为了更简便的说明，这里仅建模到一次项，即线性形变模型：

（4）

其中，

表示线性形变速率，/>

表示截距，/>

为波长，/>

将形变量转化为相位。当研究区域的形变演变状态更复杂时，可以将形变模型提升到二阶或更高阶项。

类似地，将剩余形变分量表示为

，则其包含模型之外的形变演变状态。则完整的形变相位可以表示为：

（5）

最后，完整的单个SAR获取相位可以表示为：

（6）

其中，

表示时间建模的相位分量，包含/>

和/>

。

表示未进行时间建模的分量，主要包含/>

，/>

以及相位噪声（包括残余轨道误差相位、地形误差相位、相位解缠误差以及失相干噪声等），除了其中的未建模的形变分量，其余分量均表现出时间随机性。

基于以上分析，假设有N个SAR获取相位

，构成M个干涉相位/>

，对于空间上任意一点/>

，则可以构建干涉观测方程：

（7）

其中：

（8）

（9）

（10）

其中，下标

表示第/>

个SAR获取，/>

表示所有/>

贡献的干涉相位分量。

步骤102：逐点分析联合求解形变速率和季节性对流层延迟模型参数：

对干涉图进行分块，分块大小选择约12.5km，块与块的重叠区域设置为2.5km。在实际应用中，分块参数根据数据的具体尺寸可以调整。块内参考点结合相干性、所处位置、地物类型、利用基本方法得到的线性形变速率等综合因素分析确定。

在每一干涉图块内逐点分析，求解干涉观测方程（7）中的未知参数，即线性形变速率和季节性对流层延迟模型参数。

步骤103：季节性对流层延迟模型参数空间一致性修正：

求解干涉观测方程（7）后，由于大气异质性、非线性形变、紊流、相位噪声等因素的影响，季节调制参数

在空间上可能并不一致。而通过干涉观测方程（7）只能求出被空间异质性调制的参数/>

，其中

随空间变化，因而不能直接对季节调制参数进行空间一致性修正。为此，通过求解参数的比值来消掉空间变化的/>

项：

（11）

其中：

（12）

根据空间统计分布直方图确定，具体为：分别计算比值

的空间统计分布直方图，其中最收敛到0的统计直方图所对应的/>

的取值即为/>

中的最大值。

进而，通过

的众数来对季节调制参数进行空间一致性修正，优化后的季节调制参数用/>

表示。

步骤104：重新逐点分析联合求解形变速率和季节性对流层延迟模型参数：

在每个干涉图块重新逐点求解干涉观测方程（7），此时式（8）变为：

（13）

其中：

（14）

进而，

即为时间域分析期望求得的季节性对流层延迟。

步骤105：形变区域识别：

可以注意到，在季节调制参数空间一致性优化前后共求了两次形变速率，记为

和/>

以作区分，并且可以得到两次模型拟合方差记作/>

和/>

。则具有线性形变模式的形变区域/>

（简称为线性形变区域）以及具有非线性形变模式的形变区域/>

（简称为非线性形变区域）可以通过以下准则进行识别：

（15）

其中，

为满足条件/>

的区域，/>

表示空间高通滤波，/>

表示取绝对值操作；/>

为阈值，根据/>

的全局精度来定，可以取值为/>

，其中/>

为/>

的标准偏差。

步骤106：季节性对流层延迟去除：

由于时间域矫正采用的是逐点分析法，所以非线性形变的存在可能会使得估计的季节调制参数存在误差，且误差的分布范围与非线性形变区域一致，应对其赋予掩膜以标识为需要进一步处理的区域。此外，部分像素在部分SAR获取中的相干性较低导致时间域参数估计不满足奈奎斯特采样定理，也需要对其赋予掩膜。接下来根据空间相关性信息对掩膜区域的参数进行克里金插值得到完整区域的参数。最后，根据关系

得到建模的大气对流层延迟相位。通过最小二乘反演相位时间序列并减去/>

得到时间域矫正后的相位时间序列。

步骤107：空间域矫正：

识别出潜在的形变区域并得到时间域矫正后的相位时间序列后，在形变区域形成掩膜，接下来对残余对流层延迟

进行估计并去除。

残余对流层延迟相位的去除在逐个相位时间序列上进行。首先根据相位与高程的相关性确定估计窗口的大小：改变窗口的大小，当相关性达到最大的时候，所得窗口大小即为估计残余对流层延迟的窗口大小。由于在以上步骤已经去除了大部分对流层延迟，残余的对流层延迟的幅度会大大减小，这样就可以采用以往研究中通常采用的线性模型进行估计：

（16）

其中，

是地形高程，/>

是截距，/>

是待求解系数。

值得注意的是，由于已经识别出了潜在的形变区域，因此只需要对形变区域或感兴趣的区域进行残余对流层延迟去除以及形变时间序列反演即可。这出于两点考虑：一是误差随着原理参考点而累积，二是减小运算量。

得到形变时间序列后，可以对其进行时间上的平滑以达到去除对流层延迟中紊流分量影响的目的。

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

实施例1

本发明采用C波段哨兵一号A星（Sentinel-1A）升轨数据进行实验验证。数据时间覆盖范围为2019年10月17日到2020年11月16日，采样间隔为12天。研究区域覆盖四川鲜水河和庆大河两个个流域。该区域由于地势起伏较大，气候湿润，季节变化显著，因此大气垂直分层现象明显。此外，地质条件、气候条件以及人类活动综合作用，使得该区域滑坡频发。

为了体现本发明的效果，将其它方法也应用于所选数据进行对比，包括：不对大气延迟进行处理，采用高程—相位线性模型以及时空滤波处理，根据ERA5数据进行精确插值处理，这三种对比方法。本发明所提出的方法及对比方法用于大气对流层延迟去除后的形变速率结果如图2所示。

当不对大气对流层延迟进行处理的时候（图2中的(a)），对流层分层延迟干扰形变速率，使得形变速率图中存在明显的地形相关性，这样的结果与地表真实形变状态严重不符。当用高程—相位线性模型以及时空滤波处理后（图2中的(b)），大气延迟的干扰大大削弱，但形变速率图种依然存在由残余大气延迟导致的形变速率估计误差，使得形变区域仍然难以准确识别。当用ERA-5数据对大气延迟进行矫正后（图2中的(c)），同样存在明显的形变速率估计误差。当用本发明提出的方法进行处理后（图2中的(d)），大气延迟导致的形变速率估计误差大大降低，呈现出清晰的形变区域。

为了进一步说明本发明的效果，在图3中展示了分别在形变区域、非线性形变区域以及稳定区域选取的像素点的形变时间序列，分别为图3中的(a)-(c)。所有图中分别展示了不做处理得到的原始形变时间序列、时间域矫正后的形变时间序列、空间域矫正后的形变时间序列，即本发明所展示方法得到的最终形变时间序列，以及时间域估计的形变速率。在每种情况下都能观察到，时间域矫正去除了原始形变时间序列中随季节波动的成分，空间域矫正后将残留的对流层延迟进一步削弱从而得到比较符合实际情况的形变时间序列。而且时间域处理得到的形变速率与矫正后的形变时间序列的速率较为吻合。

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

Claims

1.一种时间域和空间域两步式InSAR对流层延迟矫正方法，其特征在于，包括如下步骤：

步骤1、预处理及建立干涉相位信号模型；

步骤3、季节性对流层延迟模型参数空间一致性修正；

步骤5、形变区域识别；

步骤6、季节性对流层延迟去除；

步骤7、空间域矫正。

2.根据权利要求1所述的一种时间域和空间域两步式InSAR对流层延迟矫正方法，其特征在于：所述步骤1中，所述预处理包括干涉组合，根据高相干面积比筛选干涉对，解缠干涉相位，矫正解缠误差以及DEM误差，去趋势性相位；经过预处理后，干涉相位中剩余相位成分包含大气延迟和形变相位，根据随时间变化情况对二者分别在时间域建模；接下来用建模的信号构建干涉观测方程。

3.根据权利要求2所述的一种时间域和空间域两步式InSAR对流层延迟矫正方法，其特征在于：所述步骤2包括：对干涉图进行分块，每个干涉图块内各自设置参考点，在每一干涉图块内，逐点求解干涉观测方程得到形变速率和季节性对流层延迟模型参数。

4.根据权利要求3所述的一种时间域和空间域两步式InSAR对流层延迟矫正方法，其特征在于：所述步骤3包括：根据参数比众数确定每一干涉图块内季节性对流层延迟模型参数空间一致的参数比关系。

5.根据权利要求4所述的一种时间域和空间域两步式InSAR对流层延迟矫正方法，其特征在于：所述步骤4包括：将空间一致性的参数比关系带入干涉观测方程，重新在每一干涉图块内逐点分析求解观测方程得到形变速率和季节性对流层延迟模型参数。

6.根据权利要求5所述的一种时间域和空间域两步式InSAR对流层延迟矫正方法，其特征在于：所述步骤5包括：根据步骤2和步骤4得到的形变速率的结果以及拟合方差分别识别具有线性形变模式的形变区域以及具有非线性形变模式的形变区域。

7.根据权利要求6所述的一种时间域和空间域两步式InSAR对流层延迟矫正方法，其特征在于：所述步骤6包括：对具有非线性形变模式的形变区域的季节性对流层延迟模型参数依据空间相关性信息进行插值修正，进而得到建模的季节性对流层延迟，并从反演的相位时间序列中去除。

8.根据权利要求7所述的一种时间域和空间域两步式InSAR对流层延迟矫正方法，其特征在于：所述步骤7包括：在每个形变区域根据高程与相位的相关性峰值自适应确定窗口大小，并根据窗口内的高程-相位线性模型在逐个相位时间序列中去除残余对流层延迟相位。