CN111856459B - 一种改进的DEM最大似然约束多基线InSAR相位解缠方法 - Google Patents

Abstract

本发明涉及一种改进的DEM最大似然约束多基线InSAR相位解缠方法包括：S1、获取同一区域的多张干涉图，针对每张干涉图分别构建每个同名像素点所对应的概率密度函数；S2、将同名像素点对应的概率密度函数按照组合关系进行组合得到联合概率密度函数，求取一定区间内函数极大值所对应的相位值，形成相位集合；S3、计算得到相位集合的极差，与噪声判定阈值进行比较，在大于噪声判定阈值的像素点处加入外部DEM进行约束；S4、求得不同像素点对应的联合概率密度函数的最大值解，得到最终的解缠相位。与现有技术相比，本发明能够有效地抑制噪声影响，并减少对外部DEM的依赖性，得到精度较高的反演结果。

Description

一种改进的DEM最大似然约束多基线InSAR相位解缠方法

技术领域

本发明涉及多基线InSAR相位解缠技术，尤其是涉及一种改进的DEM最大似然约束多基线InSAR相位解缠方法。

背景技术

合成孔径雷达干涉测量(Interferometry Synthetic Aperture Radar，InSAR)技术在近年来发展较为迅猛，但是传统的单基线InSAR技术还是会受到大气效应、几何畸变等影响。基线长度也会影响到相位解缠的复杂度以及精度。当空间基线长度较短时，其对应的相位干涉条纹较为稀疏，相位解缠较为容易，但是其对应的数字高程模型精度较低。反之，虽然长基线的相位干涉条纹比较密集，对应的相位解缠复杂，但是其对应的数字高程模型精度较高。此外，由于Itoh条件的约束，对于陡峭复杂的地形或断续的场景，如悬崖山谷、高层建筑等，传统的单基线InSAR相位解缠方法反演DEM精度会显著降低甚至完全失效，导致最终没有办法实现高精度的高程测量。基于传统单基线InSAR的优缺点，多基线InSAR系统通过结合不同基线所获取的多个通道的数据，从而反演得到精度较高的高程信息。多基线InSAR扬长避短，具有极大的应用研究价值。例如，多基线InSAR系统充分结合了短基线对应干涉图相位解缠容易以及长基线对应干涉相位图反演DEM精度高的优势，有效减小了相位解缠的复杂度，提升了DEM的反演精度，对复杂地形三维测绘问题的解决起到了很大作用；而且，由于数据来源于多个通道，这能够降低或抑制噪声、大气效应以及几何畸变等因素的影响，进而得到高精度的三维地形信息。

国外有关多基线干涉处理的研究大多是采用多幅InSAR数据来提高干涉图质量和相位解缠精度或者直接利用多幅InSAR数据来估计侧视角以获得高程信息。而国内在多基线InSAR数据处理的研究上虽已有一些科研成果，但现今大多实验研究水平尚且还在起步期，现有多基线InSAR技术多是基于仿真实验进行研究，对实测真实数据的处理技术还有待完善。例如基于差分相位的多基线InSAR算法解缠精度受最短基线干涉图解缠精度影响较大，且存在误差累计问题；基于多基线最小二乘算法对噪声比较敏感，存在高程压缩现象；而基于多基线极大似然估计法存在信噪比低时解缠误差大，反演结果毛刺现象严重的弊端。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供一种改进的DEM最大似然约束多基线InSAR相位解缠方法。

本发明的目的可以通过以下技术方案来实现：

一种改进的DEM最大似然约束多基线InSAR相位解缠方法，包括：

S1、获取同一区域的多张干涉图，针对每张干涉图分别构建每个同名像素点所对应的概率密度函数；

S2、将同名像素点对应的概率密度函数按照组合关系进行组合得到联合概率密度函数，求取一定区间内函数极大值所对应的相位值，形成相位集合；

S3、计算得到相位集合的极差，与噪声判定阈值进行比较，在大于噪声判定阈值的像素点处加入外部DEM进行约束；

S4、求得不同像素点对应的联合概率密度函数的最大值解，得到最终的解缠相位。

优选的，所述每个同名像素点所对应的概率密度函数为：

其中，φ₀为真实相位的估计值，φ_n，norm为第n张干涉影像的相位值，|γ_n|为第n张相干系数图值的绝对值，φ_n为第n张干涉影像的真实相位值，ξ_n为相位比例因子。

优选的，所述ξ_n为：

其中，λ为雷达波长，B_n⊥、B_l⊥为不同干涉图对应的垂直基线长度，R为雷达到目标的距离，φ_l为第l张干涉影像的真实相位值，h为高程值，θ为雷达入射角，m为发射-接收因子，当雷达天线“一发多收”时，m＝1，当雷达天线“一发一收”时，m＝2。

优选的，所述联合概率密度函数为：

F_MB(φ|φ₀)＝Πpd_fnorm(φ_n，norm|φ₀)

其中，φ对应为真实干涉图相位。

优选的，所述噪声判定阈值为：

其中，λ为雷达对应波长，R(i，j)、θ(i，j)分别为第(i，j)像素点所对应的斜距及入射角，h_threshold为高程阈值，B_⊥为垂直基线长度。

优选的，所述在大于噪声判定阈值的像素点处加入外部DEM进行约束包括：在大于噪声判定阈值的像素点处加入外部DEM的高程信息作为先验信息，再通过高程反演得到DEM。

优选的，所述大于噪声判定阈值的像素点实际所对应的相位值的先验概率分布函数为：

其中，N为滑动窗口的大小，T为滑动窗口内像素个数，φ_DEM，i为滑动窗口内各像素对应的DEM转化后的相位值，σ_local为滑动窗口内所有DEM转化后相位值的标准差，σ_h为外部DEM转化后相位值的已知标准差。

优选的，所述外部DEM采用各干涉图得到的DEM并基于相干系数加权进行融合后得到。

优选的，所述S1中同一区域的多张干涉图通过将得到的影像依次经过影像配准过程、平地相位去除步骤、滤波处理过程得到。

优选的，所述影像配准过程中从得到的影像中选取一张主影像，对于其余影像分别进行：

通过卫星轨道信息进行粗配准得到近似的待匹配点；基于复数值的相似性测度，选取相干系数作为参考，采用窗口移动的方法找出两张影像中相干系数最大的点对作为同名像素；通过迭代拟合，剔除残差值较大的点，进而求取主辅影像最佳配准多项式模型和系数，根据最终计算出的每个像素的偏移量得到重采样之后的辅影像。

与现有技术相比，本发明首先通过组合方式判断每个像素点受噪声影像的程度，然后采用外部DEM在噪声影响严重的像素点处进行约束，缩短模糊周期，提高了高程重建的精度，有利于保持无噪声影响的像素点的准确性，又能够抑制噪声对相位解缠的影响；通过TerraSAR InSAR数据实验，对实验结果定性及定量的分析可以看出，本发明的方法通过在部分像素点处加入外部参考DEM作为约束，能够有效地抑制噪声影响，并减少对外部DEM的依赖性，有利于高精度DEM的生成。

附图说明

图1为相位集合中最大值与最小值差值示意图；

图2为实施例中研究区域示意图；

图3为实施例中基于外部DEM的局部约束的多基线最大似然法流程图；

图4为实施例中多基线最大似然反演DEM结果图；

图5为实施例中三组实验高程差分图；

图6为实施例中差异数量统计分布图。

具体实施方式

下面结合附图和具体实施例对本发明进行详细说明。本实施例以本发明技术方案为前提进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

多基线干涉SAR技术结合了长短基线的优点，从而可获得较高的高程测量精度。多基线最大似然相位解缠算法通过联合多组不同干涉相位的概率密度函数，并求取联合后概率密度函数的极大值来估计相位，其应用颇为广泛。但是在相位噪声存在的情况下，该算法不能够完全得到像素点缠绕相位的准确值，采用基于外部DEM对每个像素进行约束的多基线最大似然算法虽然能够成功的抑制噪声影响，但是最终得到的DEM对外部DEM的依赖性较大，当外部DEM精度较低时可能会使得反演精度甚至低于未加约束前的结果。

本申请提出一种改进的DEM最大似然约束多基线InSAR相位解缠方法，对不同影像对生成的干涉图，先通过组合方式判断干涉图上各像素点受噪声影响的程度，再在噪声影响程度较大的像素点处加入外部参考DEM作为约束进行高程反演实验，有利于保持无噪声影响的像素点的准确性，又能够抑制噪声对相位解缠的影响。

本方法包括：

S1、获取同一区域的多张干涉图，针对每张干涉图分别构建每个同名像素点所对应的概率密度函数；

S2、将同名像素点对应的概率密度函数按照组合关系进行组合得到联合概率密度函数，求取一定区间内函数极大值所对应的相位值，形成相位集合；

S3、计算得到相位集合的极差，与噪声判定阈值进行比较，在大于噪声判定阈值的像素点处加入外部DEM进行约束；

S4、求得不同像素点对应的联合概率密度函数的最大值解，得到最终的解缠相位。

在多基线干涉SAR条件下，若同一区域对应有n张干涉图，则对于同一地物P而言，每个干涉图均对应一个概率密度函数，由于不同的干涉相位与其相应的基线长度近似成比例关系，即有：

其中，λ为雷达波长，B_n⊥、B_l⊥为不同干涉图对应的垂直基线长度，R为雷达到目标的距离，φ_l为第l张干涉影像的真实相位值，h为高程值，θ为雷达入射角，m为发射-接收因子，当雷达天线“一发多收”时，m＝1，当雷达天线“一发一收”时，m＝2。

基于多基线最大似然相位解缠方法利用了干涉图干涉相位的概率分布，其算法的核心原理在于不同的干涉相位之间相互独立，对于单基线情况而言，干涉相位存在如下的概率密度函数式：

式中：|γ|代表相干系数值的绝对值，φ为真实干涉相位，φ₀为真实相位的估计值。从上式看出，由于概率密度函数对相位差作了余弦运算，所以其周期为2π。当解缠后相位φ与缠绕相位φ₀存在φ＝φ₀+2kπ(k为正整数)的关系时，概率密度函数取得最大值。故可以通过检测概率密度函数的最大值来进行相位估计。因为函数存在周期，所以对于单基线而言，可能存在多个最大值解。

对于多基线而言，联立式(1)和式(2)得到每个像素所对应的概率密度函数如式(3)所示，即每个同名像素点所对应的概率密度函数：

其中，φ_n，norm为第n张干涉影像的相位值，|γ_n|为第n张相干系数图值的绝对值，φ_n为第n张干涉影像的真实相位值，ξ_n为相位比例因子。

所以对于多基线而言，若将多张干涉图的概率密度函数相乘即得到联合概率密度函数F_MB为：

F_MB(φ|φ₀)＝∑pdf_norm(φ_n，norm|φ₀) (4)

其中，φ对应为真实干涉图相位。

相位的极大似然估计值为：

在有噪声存在的情况下，多基线极大似然相位解缠算法会受到影响，不能够完全恢复像素点缠绕相位的真实值。解缠结果往往出现离散分布的尖锐毛刺，这是因为部分像素点含有噪声，噪声的存在使得在求取联合概率密度函数的极大值时出现周期跳变，从而表现为该像素点与周围像素点的值差别较大。为了约束相位噪声造成的极值周期性转移的情况，可采用低分辨率DEM作为先验条件的最大后验估计相位解缠方法。本方法在该方法的基础上先随机组合概率密度函数计算出各曲线极大值对应的相位，形成相位集合，并将集合内最大值与最小值的差值

与给定阈值进行比较，如图1所示。当大于阈值则需要加入外部DEM作为约束。

当获取一幅DEM后，首先结合雷达系统参数信息将外部DEM高程值转化为相位：

其中，H为雷达天线高度，h_DEM为外部DEM所对应的高程值。

则多基线InSAR的相位最大后验估计为：

式中：pdf(φ₀)是从外部DEM中获取的相位先验分布，pdf_norm(φ_n，norm|φ₀)为每个同名像素点所对应的概率密度函数，K为干涉相位图数量。

设外部DEM中所对应的SAR干涉图中一个像素的相位值为φ_DEM，i，外部DEM对应相位的标准差为σ_DEM，则在该像素点实际所对应的相位值的先验概率分布函数为：

其中，N为滑动窗口的大小，T为滑动窗口内像素个数，φ_DEM，i为滑动窗口内各像素对应的DEM转化后的相位值，σ_local为滑动窗口内所有DEM转化后相位值的标准差，σ_DEM为外部DEM转化后相位值的已知标准差，σ_h为外部DEM转化后相位值的已知标准差。

噪声判定阈值可采用：

其中，λ为雷达对应波长，R(i，j)、θ(i，j)分别为第(i，j)像素点所对应的斜距及入射角，h_threshold为高程阈值，B_⊥为垂直基线长度。

S1中同一区域的多张干涉图通过将得到的影像依次经过影像配准过程、平地相位去除步骤、滤波处理过程得到。

多基线最大似然法通过联合多组不同干涉相位的概率密度函数，并求取联合后概率密度函数的极大值来估计相位，而基于外部DEM的多基线最大似然估计法是一种最大后验估计的方法，该方法是在多基线最大似然估计的基础上在每个像素点处加入目标高程的先验分布，其先验信息取自外部的DEM(后文将该方法称为“全局约束的多基线最大似然法”)。最大后验估计同最大似然估计具有较强的关联性，前者主要将所要估计的高程的先验分布融入到后者，可将前者等价于规则化的最大似然估计。

实施例

本实施例中，试验区位于伊比利亚半岛东北部西班牙巴塞罗那(Barcelona)地区(约1°53′00″～1°55′00″E，41°38′00″～41°40′00″N)，覆盖区域如图2矩形框所示。试验区主要覆盖山地地区，此区域有相当大的部分被植被覆盖，研究区域的高程分布在海拔200m至800m之间的这种多植被覆盖地域。TerraSAR-X影像的长波段雷达影像能够更好地避免空间以及时间失相干性的影响，且可穿透植被测得地表的高程。

实验数据采用由德国宇航中心(German Aerospace Center,DLR)和EADS Astrium及Infoterra公司共同开发的TerraSAR-X卫星摄取的几组L波段(λ＝0.031m)SAR数据。数据是在SM(Stripemap)模式下获取的升轨影像数据，源自联合国环境规划署国际环境资料系统(Infoterra)网站。实验共使用了四景TerraSAR影像数据，其详细信息如表1所示

表1实验数据的基本信息

	DATA 1	DATA 2	DATA 3	DATA 4
					获取时间	2009/01/29	2009/02/09	2009/02/20	2009/03/14
近距端入射角	33.869deg	33.871deg	33.871deg	33.872deg
					远距端入射角	36.687deg	36.686deg	36.685deg	36.688deg
距离向采样间隔	0.909m	0.909m	0.909m	0.909m
					方位向采样间隔	1.888m	1.888m	1.888m	1.888m
轨道编号	9028	9195	9028	9696

理论上上述4景数据可任意两两组合形成6幅干涉图，但是为了避免干涉图配准引起的误差，并结合空间及时间的相干性，本实施例选用了2009/02/20的影像作为三对影像组合的主影像，其余三张影像分别作为辅影像。具体的影像对信息如表2所示。

表2 InSAR测量SAR影像对基本信息

	干涉图1	干涉图2	干涉图3
				主影像获取时间	2009/02/20	2009/02/20	2009/02/20
辅影像获取时间	2009/02/09	2009/01/29	2009/03/14
				时间基线	11days	22days	22days
垂直基线	150.81m	86.72m	-35.61m
				高程模糊度	-36.73m	-63.87m	155.57m
平均相干系数	0.55	0.56	0.43

为了验证单基线InSAR DEM的精度，本实施例还使用了实验区域对应的ASTERGDEM作为参考数据。

对生成的三组干涉图进行联合，共进行了三组多基线反演DEM高程实验。第一组实验采用多基线最大似然估计算法对高程进行反演(无约束的多基线最大似然法)，第二组实验通过在每个像素点处加入外部DEM高程值作为先验信息，再进行最大似然高程估计反演DEM(全局约束的多基线最大似然法)，第三组实验首先通过判断噪声对于每个像素的影响程度，然后在噪声影响较大的像素点处加入外部DEM的高程信息作为先验信息，最终反演得到DEM(局部约束的多基线最大似然法，即本申请提出的)。具体流程如图3所示。

图3关键步骤中对应的方法如下所示：

(1)在影像配准过程中，本实施例通过卫星轨道信息先进行粗配准得到近似的待匹配点，然后基于复数值的相似性测度，选取相干系数作为参考，采用64*64窗口移动的方法找出两张影像中相干系数最大的点对作为同名相对；通过迭代拟合，剔除残差值较大的点，进而求取主辅影像最佳配准多项式模型和系数，然后根据最终计算出的每个像素的偏移量即可得到重采样之后的辅影像。

(2)在平地相位去除步骤中，本实施例基于分辨率为90m的SRTM，对研究区域平地相位进行计算，然后再从干涉图中将其去除。

(3)滤波处理过程中，本实施例基于经典Goldstein的优化方法，即Zhao滤波(ZHAOChaoying et al.，2012)对干涉相位图进行去噪处理，该方法通过对经典Goldstein的部分参数进行精化处理，滤波窗口大小设置为32*32，迭代次数为6次，达到了较好的滤波效果。

(4)本实施例外部DEM采用的是三组干涉图按照常规流程得到的DEM并基于相干系数加权进行融合后的结果。

(5)在实施例数据处理中，判断阈值φ_threshold可表示为：

其中，λ为雷达对应波长，R(i，j)、θ(i，j)分别为第(i，j)像素点所对应的斜距及入射角。本实施例中，h_threshold根据对数据反复处理最终设置为11m。

本实施例采用多基线最大似然相位解缠方法(Unconstrained MLE)、全局约束的多基线最大似然相位解缠方法(Global Constrained MLE)与本申请优化的局部约束的多基线最大似然相位解缠方法(Partial Constrained MLE)分别对TerraSAR数据进行处理，得到对应地区的DEM如图4所示。

为了对三种方法生成的多基线InSAR DEM精度进行评价，采用ASTER GDEM进行后续分析。首先将结果DEM降采样至与ASTER GDEM同等分辨率，并将对应像素高程值做差(结果DEM-ASTER GDEM)，得到三幅DEM高程值差异图如图5所示。

对DEM高程差值图进行统计分析，得到三者的统计信息如表3所示。根据统计信息可以看出无论是哪种方法所获取的DEM的均方根误差值均较小，但极大值与极小值都较为异常。通过对三种方法每一个高程差值样本的概率密度进行估计，可以发现这种差异值较大的像素点个数相对较少，三种方法对应的差异值数量统计如图6所示。

表3试验区差异图统计信息(与ASTER GDEM差异值)

通过对图4进行目视解译可以看出，采用无约束的多基线最大似然估计算法直接对多张干涉图联合进行相位解缠能够反演出地形的整体轮廓，但是因为地形相对复杂，会使得反演结果存在大量毛刺现象。而加入外部低DEM的全局约束的多基线最大似然估计算法反演出来的DEM地形清晰且几乎不存在毛刺，但是在局部区域出现了与外部DEM相似的高程跳变现象。而基于外部DEM的局部约束的多基线最大似然估计算法生成的DEM没有毛刺现象，且无论在高山或者低洼地区，地形均比较连续，不再存在高程的阶跃性跳变(如图4的矩形框内所示)。通过对三者整体进行定量评估可以发现，局部约束的最大似然估计方法反演处来的DEM所对应的高程差分值的均方根相对前面两种方法有所提高，同时其高程差异在±30m、±20m、±15m三个区间内的像元比例均较高，这充分表明了本申请优化后的最大似然算法能够有效抑制噪声对反演精度的影像，减少对外部DEM的依赖性，能够提高DEM的反演精度。

Claims (9)

1.一种改进的DEM最大似然约束多基线InSAR相位解缠方法，其特征在于，包括：

S1、获取同一区域的多张干涉图，针对每张干涉图分别构建每个同名像素点所对应的概率密度函数；

S2、将同名像素点对应的概率密度函数按照组合关系进行组合得到联合概率密度函数，求取一定区间内函数极大值所对应的相位值，形成相位集合；

S3、计算得到相位集合的极差，与噪声判定阈值进行比较，在大于噪声判定阈值的像素点处加入外部DEM进行约束；

S4、求得不同像素点对应的联合概率密度函数的最大值解，得到最终的解缠相位。

2.根据权利要求1所述的一种改进的DEM最大似然约束多基线InSAR相位解缠方法，其特征在于，所述每个同名像素点所对应的概率密度函数为：

其中，φ₀为真实相位的估计值，φ_n,norm为第n张干涉影像的相位值，|γ_n|为第n张相干系数图值的绝对值，φ_n为第n张干涉影像的真实相位值，ξ_n为相位比例因子。

3.根据权利要求2所述的一种改进的DEM最大似然约束多基线InSAR相位解缠方法，其特征在于，所述ξ_n为：

其中，λ为雷达波长，B_n⊥、B_l⊥为不同干涉图对应的垂直基线长度，R为雷达到目标的距离，φ_l为第l张干涉影像的真实相位值，h为高程值，θ为雷达入射角，m为发射-接收因子，当雷达天线“一发多收”时，m＝1，当雷达天线“一发一收”时，m＝2。

4.根据权利要求2所述的一种改进的DEM最大似然约束多基线InSAR相位解缠方法，其特征在于，所述联合概率密度函数为：

F_MB(φ|φ₀)＝Πpdf_norm(φ_n，norm|φ₀)

其中，φ对应为真实干涉图相位。

5.根据权利要求1所述的一种改进的DEM最大似然约束多基线InSAR相位解缠方法，其特征在于，所述噪声判定阈值为：

其中，λ为雷达对应波长，R(i，j)、θ(i，j)分别为第(i，j)像素点所对应的斜距及入射角，h_threshold为高程阈值，B_⊥为垂直基线长度。

6.根据权利要求1所述的一种改进的DEM最大似然约束多基线InSAR相位解缠方法，其特征在于，所述在大于噪声判定阈值的像素点处加入外部DEM进行约束包括：在大于噪声判定阈值的像素点处加入外部DEM的高程信息作为先验信息，再通过高程反演得到DEM。

7.根据权利要求1所述的一种改进的DEM最大似然约束多基线InSAR相位解缠方法，其特征在于，所述外部DEM采用各干涉图得到的DEM并基于相干系数加权进行融合后得到。

8.根据权利要求1所述的一种改进的DEM最大似然约束多基线InSAR相位解缠方法，其特征在于，所述S1中同一区域的多张干涉图通过将得到的影像依次经过影像配准过程、平地相位去除步骤、滤波处理过程得到。

9.根据权利要求8所述的一种改进的DEM最大似然约束多基线InSAR相位解缠方法，其特征在于，所述影像配准过程中从得到的影像中选取一张主影像，对于其余影像分别进行：

通过卫星轨道信息进行粗配准得到近似的待匹配点；基于复数值的相似性测度，选取相干系数作为参考，采用窗口移动的方法找出两张影像中相干系数最大的点对作为同名像素；通过迭代拟合，剔除残差值较大的点，进而求取主辅影像最佳配准多项式模型和系数，根据最终计算出的每个像素的偏移量得到重采样之后的辅影像。

Images