一种基于小波滤波的小基线集干涉测量方法

1.本发明属于地表形变探测技术领域，尤其涉及一种基于小波滤波的小基线集干涉测量方法。


背景技术：

2.小基线集干涉测量技术(small baseline subset interferometric synthetic aperture radar，sbas-insar)是一种将短时、空基线和相对小的多普勒中心频率差异的干涉对进行组合，识别长时间间隔内地表稳定高相干点的相位信息，并建立时序相位模型估算地表位移量的测量技术，被广泛用于长时间地表形变探测。如何准确识别地表稳定高相干点，解算相干点位移信息，同时有效去除误差影响，是保证其形变探测精度的关键因素。
3.传统sbas-insar技术在进行地表稳定点识别时，采用“一刀切”的硬阈值方法，造成数据损失，从而导致识别出的地表稳定点数量和密度有限；在利用最小费用流方法进行相位解缠时，未进行迭代和优化，相位不连续区域解缠效果不理想，容易造成解缠错误；在估算dem误差和消大气延迟影响时，直接采用时间域高通滤波和空间域低通滤波的方式，可能引起数据低频部分和非稳定成分(如硬阈值造成的跳变等)的丢失，且没有考虑真实大气条件，导致气象条件复杂区域形变测量误差较大。上述问题容易增大地表位移估算误差，从而降低形变探测精度。
4.针对传统sbas-insar技术由于存在上述缺陷，导致地表形变探测精度低的问题，本发明提出一种基于小波滤波的小基线集干涉测量方法。


技术实现要素：

5.本发明的目的在于提供一种基于小波滤波的小基线集干涉测量方法，以解决上述背景技术中提出的传统sbas-insar技术存在的识别出的地表稳定点数量和密度有限，气象条件复杂区域形变测量误差较大，相位不连续区域解缠效果不理想等问题。
6.为实现上述目的，本发明采用以下技术方案实现：
7.一种基于小波滤波的小基线集干涉测量方法，包括如下步骤：
8.s1、获取sar图像进行处理，得到干涉图集合；
9.s2、采用软阈值规则进行稳定点的识别选取；计算干涉相位小波系数方差，并以此为标准识别干涉图中稳定像元；
10.s3、进行稳定点构网与相位解缠；利用s2中识别出的稳定像元构建德劳内三角网络，并利用最小费用流解缠算法对稀疏稳定点网络进行解缠；
11.s4、进行相位滤波与迭代修正；使用勒让德小波对相位解缠结果进行滤波估算高频干扰，修正解缠后的相位错误并去除空间失相干噪声，并重复s3，直到估算的高频干扰小于给定阈值；
12.s5、对s4步骤处理后的干涉图进行dem误差去除；将dem误差相位进行勒让德展开，估算并消除dem误差相位；
13.s6、生成地表形变时间序列；采用加权迭代最小二乘方法，把观测值重赋权作为前一步所得残差的函数(如双平方)，并使用求解广义逆矩阵的方法进行方程的求解，估算地表形变时间序列；
14.s7、对地表形变时间序列进行大气延迟成分去除，得到最终的地表形变时间序列；
15.s8、利用最小二乘拟合的方法，估算地表形变的线性形变速率。
16.优选地，所述s1中获取sar图像进行处理，具体如下：
17.获取n+1幅相同覆盖范围的sar图像，获取时间序列为t0，t1，
…
，tn，n为自然数；
18.通过sar图像具有相似的观测几何，进行图像配准、小基线组对、干涉处理、轨道误差估算、地形相位去除的处理步骤，得到k幅具有适当空间基线的干涉图集合。
19.优选地，所述s2中识别干涉图中稳定像元，具体如下：
20.s21、在复数域通过小波分析方法估算相位噪声；
21.s22、从原始相位加以去除，进行降噪；在对干涉相位进行降噪时，采用软阈值规则进行阈值控制；
22.s23、为地表分布的目标散射体建立高斯散射模型，得到干涉相位密度分布函数。
23.优选地，所述s21中小波分析方法基于小波系数方差获取噪声成分，具体如下：
24.对于像元x在第i幅干涉图中估算的干涉相位噪声相位成分为：
[0025][0026]
其中和分别表示噪声相位的实部和虚部；
[0027]
则其相关的时域方差为：
[0028]
其中c和s分别表示实部和虚部；
[0029]
干涉相位和复相位观测结果之间关系满足：
[0030]
其中，r和i分别为包含相位噪声和的初始相位观测值的实部和虚部；
[0031]
则干涉相位标准差为：
[0032][0033]
优选地，所述s22中采用软阈值规则进行阈值控制，具体如下：
[0034]
将σ作为干涉相位的参考方差，根据方差的计算公式，则
[0035][0036]
其中χ2为卡方概率密度函数；
[0037]
σ通过基于insar视线向形变的期望方差进行估算：
[0038][0039]
其中，δr为雷达侧视斜距差，σ
l
为其方差；
[0040]
对于给定显著水平α，置信区间为：
[0041][0042]
通过该检验标准的像元为作为稳定像元，否则为非稳定像元。
[0043]
优选地，所述s5中将dem误差相位进行勒让德展开，具体如下：
[0044]
对解缠后干涉图中的任意像元x(r,c)，r、c分别为该点在方位—斜距坐标系中的坐标，dem误差的相位贡献表示为：
[0045][0046]
其中，λ为雷达波长，r代表卫星与目标点的斜距，b
⊥
表示卫星两次观测的垂直基线距离，θ为雷达侧视角；
[0047]
勒让德多项式展开式形式为：
[0048][0049]
其中，pm为m阶勒让德多项式展开式；
[0050]
dem误差相位的勒让德展开形式为：
[0051][0052]
上式中前四项成分(m＝0,1,2,3)包含85％以上的dem误差，在干涉图中进行去除。
[0053]
优选地，所述s6中生成地表形变时间序列，具体如下：
[0054]
定义ψ＝[ψ1,k,ψn]
t
为n个未知形变相位组成的向量，δφ＝[1δφ,k,qδφ]
t
为k个已知解缠后相位观测值；
[0055]
干涉图中相干点相位与地表形变相位关系表示为：
[0056][0057]
其中，a为设计矩阵，为观测误差，为时间序列位移量估计值；
[0058]
上式的最小二乘解为：其中，p为权重矩阵；
[0059]
通过迭代处理计算新权重并更新系数，直到满足预设阈值为止，如δ为预设阈值。
[0060]
优选地，所述s7中对地表形变时间序列进行大气延迟成分去除，具体如下：
[0061]
利用era-5估算雷达视线方向的大气延迟相位，高程z处的los单路对流层延迟是表面高程z与参考平面高程z
ref
之间的空气折射率的积分，表示为：
[0062][0063]
其中，θ是入射角；
[0064]
rd＝287.05j
·
kg-1
·
k-1
，rv＝461.495j
·
kg-1
·
k-1
分别是干燥空气和水蒸气比气体常数；
[0065]gm
是加权平均值z和z
ref
之间的重力加速度；
[0066]
p是以pa为单位的干燥空气分压；
[0067]
e是以pa为单位的水蒸气分压；
[0068]
t是以k为单位的温度；
[0069]
常数为k1＝0.776k
·
pa-1
，k2＝0.716k
·
pa-1
，k3＝3.75
·
103k2·
pa-1
；
[0070]zref
设置为10000m；
[0071]
上式中第一项对应于延迟路径的干燥空气分量，第二项与空气湿度有关；通过上式估算并消除时间序列中的大气延迟成分，得到最终的地表形变时间序列。
[0072]
与现有技术相比，本发明的有益效果是：
[0073]
(1)、本发明提出一种基于小波滤波的小基线集干涉测量方法，该方法提升了传统sbas-insar方法识别的目标点密度和形变探测精度，实现了一种更优的小基线集干涉测量方法。
[0074]
(2)、本发明中该方法在识别地表稳定点时采用一种尺度相关的软阈值规则进行控制，避免了传统“一刀切”的硬阈值方法造成的数据损失，提升了传统sbas-insar方法识别的目标点密度。
[0075]
(3)、本发明中该方法利用勒让德小波对相位解缠结果进行滤波估算高频干扰，通过迭代修正解缠后的相位错误并去除空间失相干噪声，降低解缠误差，提高解缠精度。
[0076]
(4)、本发明中该方法采用勒让德展开式估算并去除dem误差，从而提高地表形变探测的精度。
[0077]
(5)、本发明中该方法利用era-5大气模型代替时间域高通滤波和空间域低通滤波的方式估算并去除大气延迟的影响，从而提高地表形变探测的精度。
附图说明
[0078]
图1为本发明中的一种基于小波滤波的小基线集干涉测量方法的流程图；
[0079]
图2为实施例2中采用本发明识别出稳定点与传统方法识别出稳定点的对比图；
[0080]
图3为实施例2中时间序列的解缠结果图；
[0081]
图4为实施例2中地表形变时间序列估算图；
[0082]
图5为实施例2中结果的精度验证图。
具体实施方式
[0083]
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整的描述，显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
[0084]
实施例1：
[0085]
一种基于小波滤波的小基线集干涉测量方法，包括如下步骤：
[0086]
步骤一：假设共有n+1幅相同覆盖范围的sar图像，获取时间序列为t0，t1，
…
，tn，n
为自然数。这些sar图像具有相似的观测几何，通过图像配准、小基线组对、干涉处理、轨道误差估算、地形相位去除等步骤，得到k幅具有适当空间基线的干涉图集合。
[0087]
步骤二：计算干涉相位小波系数方差，并以此为标准识别稳定像元。首先在复数域通过小波分析方法估算相位噪声，并从原始相位加以去除(即降噪)；之后，为地表分布的目标散射体建立高斯散射模型(gaussian scattering model)，得到干涉相位密度分布函数。
[0088]
在对干涉相位进行降噪时，为了避免硬阈值函数“一刀切”效应在小波域产生跳变或者抖动，采用一种尺度相关软阈值规则(scale-dependent soft thresholding scheme)进行阈值控制。该方法基于小波系数方差获取噪声成分，因此不需要先验信息，并且避免了数据损失。
[0089]
假设对于像元x在第i幅干涉图中估算的干涉相位噪声相位成分为其中和分别表示噪声相位的实部和虚部；则其相关的时域方差为其中c和s分别表示实部和虚部。干涉相位和复相位观测结果之间关系满足：其中，r和i分别为包含相位噪声和的初始相位观测值的实部和虚部。
[0090]
则干涉相位标准差为：
[0091][0092]
假设σ为干涉相位的参考方差，根据方差的计算公式，则其中χ2为卡方概率密度函数。对于给定显著水平α，置信区间为：
[0093][0094]
通过该检验标准的像元为作为稳定像元，否则为非稳定像元。需要说明的是，σ可以基于insar视线向形变的期望方差进行估算：
[0095][0096]
其中，δr为雷达侧视斜距差，σ
l
为其方差。
[0097]
研究表明，对sar数据来说，相位噪声可以表示为视数(number of looks)的函数，如(4)式所示，而sar图像的干涉相干性与视数无关，如(5)式所示。因此，基于时序干涉相位噪声统计特征的方法可以在全分辨率下鉴别稳定像元，优于传统的小基线集技术。
[0098]
[0099][0100]
其中，nr和n
az
分别表示距离向和方位向视数，snr为雷达图像信噪比，μ
db
表示相位噪声。干涉图的相干性可以表示为：
[0101][0102]
其中，e[x]表示x的期望，|u1|和|u1|分别表示形成干涉的两幅sar图像的振幅，γ表示其干涉相干性。
[0103]
步骤三：利用步骤二中识别出的稳定像元构建德劳内三角网络(delaunay triangulation network)，并采用最小费用流解缠算法对稀疏稳定点网络进行解缠。
[0104]
步骤四：使用勒让德小波对相位解缠结果进行滤波估算高频干扰，并修正相应目标点相位，然后重复步骤三，直到估算高频干扰小于给定阈值。经该步骤处理后，相位噪声中主要包含dem误差和大气延迟成分，将在后续步骤中进行去除。
[0105]
步骤五：将dem误差相位进行勒让德展开，估算并消除dem误差相位。
[0106]
对解缠后干涉图中的任意像元x(r,c)，r、c分别为该点在方位—斜距坐标系中的坐标，用δh(x)表示dem与真实地形之间的误差，那么dem误差的相位贡献可以表示为：
[0107][0108]
其中，λ为雷达波长，r代表卫星与目标点的斜距，b
⊥
表示卫星两次观测的垂直基线距离，θ为雷达侧视角。
[0109]
勒让德多项式展开式形式为：
[0110][0111]
其中，pm为m阶勒让德多项式展开式。dem误差相位的勒让德展开形式为：
[0112][0113]
上式中前四项成分(m＝0,1,2,3)包含了85％以上的dem误差，在干涉图中进行去除。经过该步骤处理后，相位噪声中主要包含大气延迟成分，将在后续步骤中进行去除。
[0114]
步骤六：采用加权迭代最小二乘(iteratively reweighted least squares，irls)方法，把观测值重赋权作为前一步所得残差的函数(如双平方)，并使用求解广义逆矩阵的方法进行方程的求解，估算地表形变时间序列。
[0115]
雷达卫星在分别时刻a和时刻b两次成像，经干涉处理生成的相位干涉图中，其相
干点的干涉相位值与地表形变相位的关系可以表示如下：
[0116][0117]
其中，ψ＝[ψ1,k,ψn]
t
为n个未知形变相位组成的向量，δφ＝[1δφ,k,qδφ]
t
为k个已知解缠后相位观测值。(10)式的矩阵形式表示为：关系可以表示为：
[0118][0119]
其中，a为设计矩阵，为观测误差，为时间序列位移量估计值。为了求解(11)式，同时也为了降低异常值(如相位解缠误差)的影响，采用加权迭代最小二乘(iteratively reweighted least squares，irls)方法，把观测值重赋权作为前一步所得残差的函数(如双平方)。为了使用irls方法，考虑(11)式最小二乘解为：
[0120][0121]
其中，p为权重矩阵。由于观测结果是相互独立、且在时间上有重叠的子集，因此在进行求逆运算时，必须使用求解广义逆矩阵的方法进行方程的求解。如此便可更为真实地评估观测噪声，并产生更为稳健的结果。irls的迭代处理是从估算观测值的新权重开始，迭代过程如(13)式所示：
[0122][0123]
其中，t为谐调因子，求解方法见holland和wetsch的文章，df为自由度(干涉图数目减去sar图像数目再加1，即k－n+1)，计算新权重和更新系数的过程会一直迭代，直到满足预设阈值(如δ为预设阈值)为止。经过n次迭代后，同时可以得到地表形变时间序列的完成方差协方差矩阵其中，和qn分别为最终方差因子和权重矩阵。
[0124]
步骤七：利用era-5(european centre for medium-range weather forecasts reanalysis v5，第五代ecmwf大气再分析数据集)估算雷达视线方向(light of sight，los)的大气延迟相位。高程z处的los单路对流层延迟是表面高程z与参考平面高程z
ref
之间的空气折射率的积分，表示为：
[0125][0126]
其中，θ是入射角，rd＝287.05j
·
kg-1
·
k-1
，rv＝461.495j
·
kg-1
·
k-1
分别是干燥空气和水蒸气比气体常数，gm是加权平均值z和z
ref
之间的重力加速度，p是以pa为单位的干燥空气分压，e是以pa为单位的水蒸气分压，t是以k为单位的温度。常数为k1＝0.776k
·
pa-1
，k2＝0.716k
·
pa-1
，k3＝3.75
·
103k2·
pa-1
，z
ref
通常设置为10000m。(14)式中第一项对应于延迟路径的干燥空气分量，第二项与空气湿度有关。通过(14)式估算并消除时间序列中的大气延迟成分，得到最终的地表形变时间序列。
[0127]
步骤八：最后利用最小二乘法估算线性形变速率。
[0128]
实施例2：
[0129]
采用实施例1中基于小波滤波的小基线集干涉测量方法进行的实施案例。
[0130]
步骤一：收集覆盖北京地区的120幅sentinel-1(c波段)slc产品，并进行小基线组对、配准、干涉处理、轨道误差估算、地形相位去除等步骤，得到235幅具有适当空间基线的干涉图集合。
[0131]
步骤二：采用(2)式软阈值规则进行稳定点的识别选取。经过对比，传统硬阈值方法处理得到的稳定相干点密度为299个/km2，本发明提出的软阈值方法识别到的相干点密度提高到532个/km2。识别得到的稳定相干点密度和数量明显提升，参见图2。图2中左侧图片为本发明所识别出的稳定点，右侧图片为传统方法所识别出的稳定点。
[0132]
步骤三：构建德劳内三角网络，采用最小费用流解缠算法对所有干涉图中稀疏稳定点网络进行解缠，解缠结果参见图3。
[0133]
步骤四：利用勒让德小波对相位解缠结果进行滤波估算高频干扰，根据数据的采样周期，设置高频干扰部分阈值为12，重复步骤三进行迭代，修正解缠后的相位错误并去除空间失相干噪声，降低解缠误差。
[0134]
步骤五：按照(9)式对步骤四中的相位解缠结果进行勒让德展开，计算前四项，估算dem误差，并从相位中去除。
[0135]
步骤六：按照(10)-(13)式迭代，将观测值重赋权作为前一步所得残差的函数(如双平方)，并使用求解广义逆矩阵的方法进行方程的求解，估算地表形变时间序列。参见图4。
[0136]
步骤七：按照(14)式逐像元估算大气延迟误差，并从时间序列中去除。
[0137]
步骤八：利用最小二乘拟合的方法，估算线性形变速率。参见图5，经过水准实测数据对比，二者相关性(r2)达到0.82，误差优于6mm/year。
[0138]
以上所述，仅用于帮助理解本发明的方法及其核心要义，但本发明的保护范围并不局限于此，对于本技术领域的一般技术人员在本发明揭露的技术范围内，根据本发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护范围之内。综上所述，本说明书内容不应理解为对本发明的限制。

技术特征：
1.一种基于小波滤波的小基线集干涉测量方法，其特征在于，包括如下步骤：s1、获取sar图像进行处理，得到干涉图集合；s2、采用软阈值规则进行稳定点的识别选取；计算干涉相位小波系数方差，并以此为标准识别干涉图中稳定像元；s3、进行稳定点构网与相位解缠；利用s2中识别出的稳定像元构建德劳内三角网络，并利用最小费用流解缠算法对稀疏稳定点网络进行解缠；s4、进行相位滤波与迭代修正；使用勒让德小波对相位解缠结果进行滤波估算高频干扰，修正解缠后的相位错误并去除空间失相干噪声，并重复s3，直到估算的高频干扰小于给定阈值；s5、对s4步骤处理后的干涉图进行dem误差去除；将dem误差相位进行勒让德展开，估算并消除dem误差相位；s6、生成地表形变时间序列；采用加权迭代最小二乘方法，把观测值重赋权作为前一步所得残差的函数，并使用求解广义逆矩阵的方法进行方程的求解，估算地表形变时间序列；s7、对地表形变时间序列进行大气延迟成分去除，得到最终的地表形变时间序列；s8、利用最小二乘拟合的方法，估算地表形变的线性形变速率。2.根据权利要求1所述的一种基于小波滤波的小基线集干涉测量方法，其特征在于，所述s1中获取sar图像进行处理，具体如下：获取n+1幅相同覆盖范围的sar图像，获取时间序列为t0，t1，
…
，t
n
，n为自然数；通过sar图像具有相似的观测几何，进行图像配准、小基线组对、干涉处理、轨道误差估算、地形相位去除的处理步骤，得到k幅具有适当空间基线的干涉图集合。3.根据权利要求1或2所述的一种基于小波滤波的小基线集干涉测量方法，其特征在于，所述s2中识别干涉图中稳定像元，具体如下：s21、在复数域通过小波分析方法估算相位噪声；s22、从原始相位加以去除，进行降噪；在对干涉相位进行降噪时，采用软阈值规则进行阈值控制。4.根据权利要求3所述的一种基于小波滤波的小基线集干涉测量方法，其特征在于，所述s21中小波分析方法基于小波系数方差获取噪声成分，具体如下：对于像元x在第i幅干涉图中估算的干涉相位噪声相位成分为：其中和分别表示噪声相位的实部和虚部；则其相关的时域方差为：其中w＝c,s，c和s分别表示实部和虚部；干涉相位和复相位观测结果之间关系满足：其中，r和i分别为包含相位噪声和的初始相位观测值的实部和虚部；则干涉相位标准差为：
5.根据权利要求4所述的一种基于小波滤波的小基线集干涉测量方法，其特征在于，所述s22中采用软阈值规则进行阈值控制，具体如下：将σ作为干涉相位的参考方差，根据方差的计算公式，则其中χ2为卡方概率密度函数；σ通过基于insar视线向形变的期望方差进行估算：其中，δr为雷达侧视斜距差，σ
l
为其方差；对于给定显著水平α，置信区间为：通过该检验标准的像元为作为稳定像元，否则为非稳定像元。6.根据权利要求1所述的一种基于小波滤波的小基线集干涉测量方法，其特征在于，所述s5中将dem误差相位进行勒让德展开，具体如下：对解缠后干涉图中的任意像元x(r,c)，r、c分别为该点在方位—斜距坐标系中的坐标，dem误差的相位贡献表示为：其中，λ为雷达波长，r代表卫星与目标点的斜距，b
⊥
表示卫星两次观测的垂直基线距离，θ为雷达侧视角；勒让德多项式展开式形式为：其中，p
m
为m阶勒让德多项式展开式；dem误差相位的勒让德展开形式为：上式中前四项成分(m＝0,1,2,3)包含85％以上的dem误差，在干涉图中进行去除。7.根据权利要求1所述的一种基于小波滤波的小基线集干涉测量方法，其特征在于，所述s6中生成地表形变时间序列，具体如下：定义ψ＝[ψ1,k,ψ
n
]
t
为n个未知形变相位组成的向量，δφ＝[1δφ,k,
q
δφ]
t
为k个已知解缠后相位观测值；干涉图中相干点相位与地表形变相位关系表示为：
其中，a为设计矩阵，为观测误差，ψ)为时间序列位移量估计值；上式的最小二乘解为：其中，p为权重矩阵；通过迭代处理计算新权重并更新系数，直到满足预设阈值为止，如δ为预设阈值。8.根据权利要求1所述的一种基于小波滤波的小基线集干涉测量方法，其特征在于，所述s7中对地表形变时间序列进行大气延迟成分去除，具体如下：利用era-5估算雷达视线方向的大气延迟相位，高程z处的los单路对流层延迟是表面高程z与参考平面高程z
ref
之间的空气折射率的积分，表示为：其中，θ是入射角；r
d
＝287.05j
·
kg-1
·
k-1
，r
v
＝461.495j
·
kg-1
·
k-1
分别是干燥空气和水蒸气比气体常数；g
m
是加权平均值z和z
ref
之间的重力加速度；p是以pa为单位的干燥空气分压；e是以pa为单位的水蒸气分压；t是以k为单位的温度；常数为k1＝0.776k
·
pa-1
，k2＝0.716k
·
pa-1
，k3＝3.75
·
103k2·
pa-1
；z
ref
设置为10000m；上式中第一项对应于延迟路径的干燥空气分量，第二项与空气湿度有关；通过上式估算并消除时间序列中的大气延迟成分，得到最终的地表形变时间序列。

技术总结
本发明公开了一种基于小波滤波的小基线集干涉测量方法，涉及地表形变探测技术领域；该方法包括如下步骤：获取SAR图像进行处理；采用软阈值规则进行稳定点的识别选取，计算干涉相位小波系数方差，并以此为标准识别干涉图中稳定像元；进行稳定点构网与相位解缠；进行相位滤波与迭代修正；对干涉图进行DEM误差去除，将DEM误差相位进行勒让德展开，估算并消除DEM误差相位；生成地表形变时间序列；对地表形变时间序列进行大气延迟成分去除，得到最终的地表形变时间序列；利用最小二乘拟合的方法，估算地表形变的线性形变速率。该方法提升了传统SBAS-InSAR方法识别的目标点密度和形变探测精度，实现了一种更优的小基线集干涉测量方法。法。法。


技术研发人员：高明亮 宫辉力 陈蓓蓓 李小娟 周超凡 冷靖 李翔 叶楹桐
受保护的技术使用者：首都师范大学
技术研发日：2023.05.05
技术公布日：2023/8/28