CN108051810B - 一种InSAR分布式散射体相位优化方法 - Google Patents

Abstract

本发明属于合成孔径雷达干涉测量数据处理技术领域，公开了一种InSAR分布式散射体相位优化方法，所述方法包括：获取同一地区的SAR影像堆栈，选取主影像，将剩余影像依次配准到主影像几何空间；探测影像内像元的振幅序列时间异变值，将其剔除；构造置信区间，判断影像内像元为中心的固定窗口内各像元的振幅序列均值是否落入置信区间，若是，则认为该点为中心像元的统计同质点SHPs；计算所述像元的样本相干系数估值，并进行偏差纠正；根据阈值，识别影像内的分布式散射体；构建所述分布式散射体相位序列的最大似然估计函数；求解所述最大似然函数，得到最优解，完成相位优化。与现有技术相比，本发明具有精度高、适用范围广、运算效率高等优势。

Description

一种InSAR分布式散射体相位优化方法

技术领域

本发明属于合成孔径雷达干涉测量数据处理技术领域，具体涉及一种InSAR分布式散射体相位优化方法。

背景技术

合成孔径雷达干涉测量技术(Interferometry Synthetic Aperture Radar，InSAR)因其覆盖范围广、测量精度高等优势，已逐渐成为大地测量领域的新宠。为消除时间、空间失相干与大气延迟效应对InSAR技术精度的影响，基于时序分析的高精度差分干涉测量技术应运而生。现有时序分析方法中，应用最为广泛的为PSI(Persistent ScattererInterferometry)技术，现有研究表明，该技术在人工地表可获得高精度高分辨率的地表形变信息，但是在非人工地表精度较低，并不适用于自然地表。究其原因在于永久散射体(Persistent Scatterer，PS)多分布于人工地表，极少分布于自然地表。

分布式散射体(distributed Scatterers，DS)分布均匀且广泛，为提高多时相InSAR技术观测点空间密度与形变解算精度，近年来研究人员将研究重点转移至基于分布式散射体的雷达干涉测量技术，该技术包含两个主要步骤：一是同质点探测，二是DS相位优化。通过同质点探测识别DS，增加观测点的空间密度；通过DS相位优化将DS与PS规划到同一解算框架，提高形变解算精度。目前为止，该领域最卓越的成果为SqueeSAR^TM技术。

但是该技术存在一些问题。首先，该技术利用K-S检验进行同质点探测，K-S检验探测精度易受样本分布影响，且并不普适，其次，该技术构建DS相位序列最大似然估计函数时未考虑样本相干性矩阵存在偏差的问题，且利用BFGS进行DS相位优化的解算效率很低。

发明内容

为解决现有技术的不足，本发明的目的在于提供一种InSAR分布式散射体相位优化方法，能够提高同质点探测与分布式散射体相位优化的精度与效率。

为了实现上述目标，本发明采用如下的技术方案：一种InSAR分布式散射体相位优化方法，其特征在于，包括以下步骤：

步骤一，利用SAR卫星对同一地区周期性成像获取所述地区的SAR影像堆栈，选取主影像，将剩余影像依次配准到主影像几何空间；

步骤二，探测SAR影像内所有像元的振幅序列时间异变值，并将其剔除；

步骤三，利用剔除时间异变后的振幅期望点估计

与单视时域平均振幅影像变差系数

构造置信区间，逐个判断SAR影像内每个像元为中心的固定窗口内各像元的振幅序列均值是否落入置信区间，若是，则认为该点为中心像元的统计同质点SHPs；

步骤四，根据像元的所述统计同质点SHPs计算所述像元的样本相干系数估值，并对所述样本相干系数估值进行偏差纠正；

步骤五，根据阈值，在偏差纠正后的无偏相干系数分布图中识别SAR影像内的分布式散射体；

步骤六，利用包含相位序列的投影矩阵与样本相干性矩阵，构建所述分布式散射体相位序列的最大似然估计函数；

步骤七，求解所述分布式散射体相位序列的最大似然函数，得到最优解，完成相位优化。

进一步地，步骤一中所述选取主影像的方法为：

依次计算SAR影像堆栈中各影像的总体相干系数ρ_total，计算公式如下：

其中：

式中，T、T_⊥、T_DC分别为时间基线、空间基线与多普勒相位中心基线，T^C、

分别为时间基线、空间基线与多普勒相位中心基线的极限值；

取最大的总体相干系数ρ_total所对应的影像作为主影像。

进一步地，步骤二中所述探测SAR影像内所有像元的振幅序列时间异变值的方法为：对于SAR影像内的任意像元，以其振幅序列为样本，记Q₁为样本下四分位点，Q₂为样本中位数，Q₃为样本上四分位点，Q₃-Q₁为样本四分位距IQR；记Q₃+h_u(MC)IQR为样本上内限，Q₁-h_l(MC)IQR为样本下内限，位于所述上内限和下内限之外的样本值即为振幅序列时间异变值；

h(MC)的计算公式为：

其中，

式中，x_i为样本中小于Q₂的任意值，x_j为样本中大于Q₂的任意值；

核函数g为：

进一步地，步骤三中所述置信区间为SAR影像内的任意像元p振幅期望点估计的一个置信水平为1-α的置信区间：

进一步地，步骤三中所述固定窗口的大小为15像素×15像素。

进一步地，步骤四中所述像元的样本相干系数估值

的计算公式为：

式中，K为窗口内的SHPs，*表示共轭相乘，s₁(t)、s₂(t)分别表示像元在两幅影像上的复数值。

进一步地，步骤四中所述对样本相干系数估值进行偏差纠正的方法为：

获取样本相干系数估值

与真实值|γ|的偏差B_R，B_R的计算公式为：

式中，

为抽样样本的相干系数，R为抽样样本个数；

计算修正后的样本相干系数

计算公式为：

进一步地，所述抽样样本个数R为200。

进一步地，步骤五中所述阈值为0.6。

进一步地，步骤六中所述投影矩阵为：

P_A(θ)＝A(θ)(A(θ)^HA(θ))^-1A(θ)^H

式中，A(θ)为空间转向序列，θ为相位序列，H表示艾尔米特转置；

所述样本相干性矩阵为：

式中，Ω表示分布式散射体D的同质点合集，||·||为归一化算子，p为分布式散射体D的同质点，x(P)为复数序列。

本发明的有益之处在于：

与现有技术相比，本发明所提供的一种InSAR分布式散射体相位优化方法具有精度高、适用范围广、运算效率高等优势，具体表现为：其一，本发明提出的探测同质点方法，不受样本分布的影响，适用性更高，通过剔除振幅序列时间异变值提高了探测精度，通过构建置信区间提高了探测效率；其二，本发明考虑样本相干系数估计的有偏性，并采用合适的方法纠正了偏差，从而提高分布式散射体识别与相位序列最大似然估计函数构造的精度；其三，综合交替投影与影像总体相干系数求解相位序列最大似然估计函数，可以快速有效的逼近全局最优解。

附图说明

图1是本发明一个实施例的InSAR分布式散射体相位优化方法的流程示意图。

具体实施方式

以下结合附图和具体实施例对本发明作具体的介绍。

参照图1所示，本发明一种InSAR分布式散射体相位优化方法，其特征在于，包括以下步骤：

S10，利用SAR卫星对同一地区周期性成像获取所述地区的SAR影像堆栈，选取主影像，将剩余影像依次配准到主影像几何空间。其中，所述主影像的选取方法为：

依次计算SAR影像堆栈中各影像的总体相干系数ρ_total，计算公式为：

其中：

式中，T、T_⊥、T_DC为时间基线、空间基线与多普勒相位中心基线，T^C、

为时间基线、空间基线与多普勒相位中心基线的极限值，常规取值为：T^C＝5years、

取最大ρ_total所对应的影像作为主影像。

S20，探测SAR影像内所有像元的振幅序列时间异变值，并将其剔除。

所述振幅序列时间异变值探测方法为：

对于SAR影像内的任意像元，以其振幅序列为样本，记Q₁为样本下四分位点，Q₂为样本中位数，Q₃为样本上四分位点，Q₃-Q₁为样本四分位距IQR。上内限为Q₃+h_u(MC)IQR为样本上内限，Q₁-h_l(MC)IQR为样本下内限，位于内限之外的样本值即为时间异变值。h(MC)的计算公式为：

其中MC是偏态分布中的一个稳健估计值，其定义为：

式中，x_i为样本中小于Q₂的任意值，x_j为样本中大于Q₂的任意值，核函数g为：

S30，利用剔除时间异变后的振幅期望点估计

与单视时域平均振幅影像变差系数

构造置信区间，逐个判断SAR影像内每个像元为中心的15像素×15像素窗口内各像元的振幅序列均值是否落入置信区间，若是，则认为该点为中心像元的统计同质点SHPs。

所述置信区间构建方法为：对于SAR影像内的任意像元p，随着影像幅数的增加，该点振幅期望的点估计

近似服从期望为μ(p)，方差为Var(A(p))/N的高斯分布，根据标准正态分布的上分位点定义与未知参数置信水平为1-α的置信区间公式可得：

式中，z表示上分位数，将上式(6)展开：

式中，Var(A(p))＝(0.52·μ(p))²。假设SAR影像像元空间不相干，单视时域平均振幅影像

的变差系数可表示为：

将式(8)带入展开式(7)，像元p振幅期望点估计的一个置信水平为1-α的置信区间为：

上述公式(9)所述的置信水平均为95％。

S40，根据像元的所述统计同质点SHPs计算所述像元的样本相干系数估值，并利用bootstrap法对所述样本相干系数估值进行偏差纠正。

由于SAR影像任意像元的相干系数|γ|很难直接求得，通常采用固定窗口内的样本相干系数估计

代替|γ|：

式中，K为窗口内的SHPs，*表示共轭相乘，s₁(t)，s₂(t)分别表示像元在两幅影像上的复数值。

所述样本相干系数估计的偏差纠正方法为：

对于样本数为K，样本相干系数为

的任意像元p，随机且有放回的对此样本进行抽样，依次生成R个样本数为K的抽样样本x^*r，r＝1，2...，R，分别计算样本相干系数

则样本相干系数

与真实值|γ|的偏差B_R可以近似为：

其中，R取200可兼顾计算精度与效率。

偏差改正后的样本相干系数为：

S50，根据阈值，在偏差纠正后的无偏相干系数分布图中识别SAR影像内的分布式散射体。所述阈值优选为0.6。

S60，利用包含相位序列的投影矩阵与样本相干性矩阵构，构建所述分布式散射体相位序列的最大似然估计函数。

对于识别到的任意分布式散射体D，其同质点具有相同相位序列θ，以D的同质点为随机变量的联合密度函数可表示为：

式(13)中，Ω表示D的同质点合集，p为D的同质点，x(户)为复数序列，A(θ)为空间转向序列，θ为相位序列，s(户)为后向反射强度量级序列，σ²I为随机噪声的方差分量。对式(13)两端取对数并忽略常数项，则相位序列对数似然函数可表示为：

对式(14)求σ²的偏导，将σ²的最大似然估计值带回上式，忽略常数项，利用最小二乘消除θ与s之间的线性关系，则相位序列的最大似然估计函数可近似为：

式(15)中，tr|·|为矩阵对角线元素求和算子，P_A(θ)＝A(θ)(A(θ)^HA(θ))^-1A(θ)^H为将s投影到A(θ)列向量空间的投影矩阵，H表示艾尔米特转置，

为D的样本相干性矩阵。

对复数序列x(P)进行归一化处理，则样本相干性矩阵

可近似为：

式(16)中，Ω表示D的同质点合集，||·||为归一化算子。

S70，求解所述分布式散射体相位序列的最大似然函数，得到最优解，完成相位优化。

利用交替投影算法循环求解相位序列的最大似然函数，若循环终止，则相位序列收敛到全局最优解。具体方法如下：

根据影像总体相干系数ρ_total对相位序列重新排序，对于相位序列为θ＝[θ₁，θ₂，...，θ_N]^T的任意分布式散射体户，在第k轮(k＞2)迭代中，将

作为已知值，则θ_i，i＝1，2，...，N的最大似然估计值可近似为：

式(17)中，

表示第k-1轮迭代中计算的相位序列估计值：

重复上述过程直至θ中所有元素被估计，则第k轮迭代后相位序列估计值为

相邻两次迭代相位序列估计值的差值小于给定阈值10^-6时，迭代终止。迭代终止时的相位序列估值即为相位序列最大似然估计函数的最优解。

此外，为降低运算负担，根据投影矩阵分解公式，将

分解为

结合矩阵迹的性质，则θ_i，i＝1，2，...，N的最大似然估计值可改写为：

式(19)中，b为单位向量：

式(20)中，||·||为归一化算子：

迭代初值的获取步骤为：

首先，对相位序列进行重新排序，主影像相位排首位，其余相位按影像总体相干系数ρ_total的大小依次排序，重新排序后的相位序列为

其次，由于干涉相位对应影像相位序列的差值，为不失一般性，将主影像相位设为零；

最后，

初值估计可表示为：

式(21)中，

重复上述过程直至相位序列中所以元素的初值被估计。

以上显示和描述了本发明的基本原理、主要特征和优点。本行业的技术人员应该了解，上述实施例不以任何形式限制本发明，凡采用等同替换或等效变换的方式所获得的技术方案，均落在本发明的保护范围内。

Claims (9)

1.一种InSAR分布式散射体相位优化方法，其特征在于，包括以下步骤：

步骤一，利用SAR卫星对同一地区周期性成像获取所述地区的SAR影像堆栈，选取主影像，将剩余影像依次配准到主影像几何空间；

步骤二，对于SAR影像内的任意像元，以其振幅序列为样本，记Q₁为样本下四分位点，Q₂为样本中位数，Q₃为样本上四分位点，Q₃-Q₁为样本四分位距IQR；记Q₃+h_u(MC)IQR为样本上内限，Q₁-h₁(MC)IQR为样本下内限，位于所述上内限和下内限之外的样本值即为振幅序列时间异变值；

h(MC)的计算公式为：

其中，

式中，x_i为样本中小于Q₂的任意值，x_j为样本中大于Q₂的任意值；

核函数g为：

并将时间异变值剔除；

步骤三，利用剔除时间异变后的振幅期望点估计与单视时域平均振幅影像变差系数构造置信区间，逐个判断SAR影像内每个像元为中心的固定窗口内各像元的振幅序列均值是否落入置信区间，若是，则认为该点为中心像元的统计同质点SHPs；

步骤四，根据像元的所述统计同质点SHPs计算所述像元的样本相干系数估值，并对所述样本相干系数估值进行偏差纠正；

步骤五，根据阈值，在偏差纠正后的无偏相干系数分布图中识别SAR影像内的分布式散射体；

步骤六，利用包含相位序列的投影矩阵与样本相干性矩阵，构建所述分布式散射体相位序列的最大似然估计函数；

步骤七，求解所述分布式散射体相位序列的最大似然函数，得到最优解，完成相位优化。

2.根据权利要求1所述的一种InSAR分布式散射体相位优化方法，其特征在于，步骤一中所述选取主影像的方法为：

依次计算SAR影像堆栈中各影像的总体相干系数ρ_total，计算公式如下：

其中：

式中，T、T_⊥、T_DC分别为时间基线、空间基线与多普勒相位中心基线，T^C、

分别为时间基线、空间基线与多普勒相位中心基线的极限值；

取最大的总体相干系数ρ_total所对应的影像作为主影像。

3.根据权利要求1所述的一种InSAR分布式散射体相位优化方法，其特征在于，步骤三中所述置信区间为SAR影像内的任意像元p振幅期望点估计的一个置信水平为1-α的置信区间：

4.根据权利要求1所述的一种InSAR分布式散射体相位优化方法，其特征在于，步骤三中所述固定窗口的大小为15像素×15像素。

5.根据权利要求1所述的一种InSAR分布式散射体相位优化方法，其特征在于，步骤四中所述像元的样本相干系数估值

的计算公式为：

式中，K为窗口内的SHPs，*表示共轭相乘，s1(t)、s2(t)分别表示像元在两幅影像上的复数值。

6.根据权利要求5所述的一种InSAR分布式散射体相位优化方法，其特征在于，步骤四中对所述样本相干系数估值进行偏差纠正的方法为：

获取样本相干系数估值

与真实值

的偏差B_R，B_R的计算公式为：

式中，

为抽样样本的相干系数，R为抽样样本个数；

计算修正后的样本相干系数

计算公式为：

7.根据权利要求6所述的一种InSAR分布式散射体相位优化方法，其特征在于，所述抽样样本个数R为200。

8.根据权利要求1所述的一种InSAR分布式散射体相位优化方法，其特征在于，步骤五中所述阈值为0.6。

9.根据权利要求1所述的一种InSAR分布式散射体相位优化方法，其特征在于，步骤六中所述投影矩阵为：

P_A(θ)＝A(θ)(A(θ)^H A(θ))^-1A(θ)^H,

式中，A(θ)为空间转向序列，θ为相位序列，H表示艾尔米特转置；

所述样本相干性矩阵为：

式中，Ω表示分布式散射体D的同质点合集，||·||为归一化算子，p为分布式散射体D的同质点，x(P)为复数序列。

Images