CN113311432B - 一种基于相位导数方差的InSAR长短基线融合相位估计方法 - Google Patents

Abstract

本发明公开了一种基于相位导数方差的InSAR长短基线融合相位估计方法，InSAR系统中基线长度的增加会提高测高精度但同时也会加大相位解缠的难度。因此在本发明中充分发挥了长短基线的优势，利用长短基线获取的同一场景的互补信息提高了相位估计的精度，基于相位导数方差准确标记了长短基线各自的不连续区域像素点，在此基础上利用长短基线各自的绝对相位中每个像素点的连续特性，自适应选取不同的融合方法，实现了长基线绝对融合相位的估计。经过仿真验证，本发明提出的方法可以有效地抑制噪声引起的解缠错误，显著减少了跳变点的个数，与传统多基线最大似然相位解缠算法相比，运算效率显著提升，最终提高了长基线域融合绝对相位的估计精度和鲁棒性。

Description

一种基于相位导数方差的InSAR长短基线融合相位估计方法

技术领域

本发明涉及合成孔径雷达干涉处理领域，具体是指一种基于相位导数方差的InSAR长短基线融合相位估计方法。

背景技术

合成孔径雷达(SAR)系统目前已经发展成为一种非常重要的微波遥感探测技术，该技术通过发射雷达波，实现对地面的主动式观测。由于雷达波是微波波段，具备较强的穿透能力，因此具有全天时、全天候工作能力的独特优势。合成孔径雷达干涉测量技术(InSAR)则是将SAR成像技术和干涉技术结合起来，通过获取两副或多副相对于同一地区的SAR图像相位信息和地形的几何关系，从而高效地获取大面积地形高程信息，是目前地形测高主要手段之一。鉴于以上优点，该技术目前已经迅速发展为一项极具应用价值的微波遥感技术，在地表测绘、形变检测、灾害监测、冰川研究等科学研究及民用领域发挥着越来越重要的作用，具有极高的研究价值。

传统单基线InSAR技术利用单一基线对应的两幅SAR主辅图像配准生成缠绕干涉相位，经过相位解缠和控制点校正后获得绝对相位，最后根据InSAR的成像空间几何关系从绝对相位反演得到观测场景的高程信息。传统的单基线InSAR系统对系统噪声、相位噪声比较敏感，并且容易受到大气效应、空间/时间去相干、数据叠掩阴影、顶底倒置等因素的影响，严重降低了测高精度。为了克服单基线InSAR技术的不足，在此基础上，多基线InSAR技术被提出，其基本原理是利用多个基线获取多个通道的数据，增加了观测信息，从而解决单基线InSAR技术难以处理的问题。

基于InSAR技术的原理分析得到：短基线的干涉条纹较稀疏，相干性高，因此比较容易解缠，但是高程反演精度较低，难以展示地形变化细节；而长基线的干涉条纹较密集，且相干性较低，增加了相位解缠的难度，但是长基线反演高程时可以反映更多的地形变化细节，测高精度较高。综上可知，单基线InSAR技术无法兼顾干涉相干性和反演高程精度的要求，即基线长度与高程细节是一个比较难以克服的矛盾。在单基线的解缠算法中，往往是基于相位连续性假设，即相邻像素之间的相位差绝对值需要小于π，而现实中存在许多复杂陡峭的地形，如高楼建筑、悬崖峭壁等等，并且缠绕相位容易受到噪声的影响，这些因素都将导致相位不连续，从而增大了单基线相位解缠的难度，降低了解缠相位的精度。多基线InSAR技术可以获取不同基线下的多幅干涉图像进行处理，这意味着可以得到更多的观测信息，有效克服单基线InSAR系统的不足。

现有的多基线InSAR相位解缠技术主要有以下几类，第一类是多基线最大似然相位解缠算法，其原理为利用不同基线下的干涉相位概率密度函数，对干涉图中每一个像素点的干涉相位实现最大似然估计，但该算法的相位估计精度受SAR工作频率、相干系数以及干涉图数量的影响，当相干系数较低或者干涉图数量较多时，相位噪声也相应增大，利用多基线最大似然相位解缠将会产生较大误差甚至失效；第二类是多基线最大后验相位解缠算法，其原理为利用多基线干涉图中相邻像素信息，采用最大后验估计算法提高反演精度，但是此类算法的局限性在于多基线下最大后验估计算法的运算时间长，内存压力大。第三类是利用多基线干涉相位与模糊数之间的关系进行相位解缠，主要是中国余数定理与其改进算法，但是这种方法的噪声鲁棒性较差，无法在实际中应用。

发明内容

本发明要解决的技术问题为：为了克服单基线InSAR系统中基线长度与反演高程精度之间的矛盾，利用长短基线获取的互补信息，本发明提供一种基于相位导数方差的InSAR长短基线融合相位估计方法，在提高相位估计鲁棒性的前提下，同时保证了高程测量精度不会因此而下降，充分发挥了长短基线各自优势。

本发明技术解决方案：一种基于相位导数方差的InSAR长短基线融合相位估计方法，包括以下几个步骤：

步骤一：利用长基线和短基线的SAR主辅图像分别进行图像配准，得到同一场景下长短基线各自的SAR主辅图像干涉相位；

步骤二：提出了一种新的移除参考地形相位方法：选择所述同一场景下的方位向上近距端和远距端的两个控制点，分别计算所述两个控制点的坐标和长短基线各自的SAR辅图像对应的天线位置坐标之间的距离，即两个控制点的斜距，再计算出所述两个控制点的斜距差，最终拟合出长短基线各自的参考地形相位；

步骤三：将长短基线各自的SAR主辅图像干涉相位均减去长短基线各自的参考地形相位，剩下的相位即为长短基线各自移除参考地形相位后的残余相位，然后计算长短基线各自移除参考地形相位后残余相位的相干系数，再依次对长短基线各自移除参考地形相位后的残余相位进行干涉图滤波、相位解缠以及控制点校正，最终得到长短基线各自移除参考地形相位后的绝对相位；

步骤四：提出了一种基于相位导数方差、跳变点判别以及图像膨胀原理的不连续区域判别方法：针对所述长短基线各自移除参考地形相位后的绝对相位，分别求出长短基线各自移除参考地形相位后的绝对相位归一化的相位导数方差，将所述相位导数方差作为质量图引导，结合绝对相位上的跳变点和图像膨胀原理，标记出长短基线各自移除参考地形相位后的绝对相位的不连续区域；

步骤五：提出了一种基于长短基线相位连续特性的自适应融合相位估计方法：根据步骤四所述的长短基线各自移除参考地形相位后的绝对相位的不连续区域和步骤三中所述的相干系数，针对长短基线长短基线各自移除参考地形相位后的绝对相位的连续特性，自适应选取融合相位估计方法，最终得到相位精度更高的长基线移除参考相位后的绝对融合相位。

所述步骤二中，一种新的移除参考地形相位方法具体实现过程为：

(1)选择所述同一场景下的方位向上近距端和远距端的两个控制点，分别计算所述两个控制点坐标和长短基线各自的SAR辅图像对应的天线位置坐标之间的距离，即两个控制点的斜距，再计算出所述两个控制点的斜距差；

(2)从InSAR空间几何关系出发，假定在SAR主辅图像对应场景的同一个方位向上，SAR主辅图像的斜距差沿距离向的变化是线性的，那么利用所述两个控制点的斜距差和SAR主辅图像同一方位向上像素点的总个数即可求解出参考地形相位线性关系表达式y＝ax+b中的参数a和b；

(3)最终的参考地形相位可以通过线性关系表达式y＝ax+b计算得到，其中y代表SAR主辅图像同一方位向上每一个像素点的参考地形相位，x代表每一个像素点在同一个方位向上的顺序；

(4)按照以上方法分别求出长短基线各自的参考地形相位，将长短基线各自的SAR主辅图像干涉相位均减去长短基线各自的参考地形相位，实现了长短基线各自的移除参考地形相位。

所述步骤四中，一种基于相位导数方差、跳变点判别以及图像膨胀原理的不连续区域判别方法具体实现过程为：

(1)计算移除参考地形相位后绝对相位的归一化相位导数方差，设置一个0～1之间的相位导数方差阈值，将移除参考地形相位后绝对相位中归一化相位导数方差大于相位导数方差阈值的像素点标记为不连续像素点，不连续像素点代表相位梯度变化剧烈的像素点，将这些不连续像素点的集合记为mark1；

(2)根据跳变点的判别方法，将移除参考地形相位后绝对相位中与相邻像素点的相位梯度差绝对值超过π的像素点判定为跳变点，跳变点可以标记出噪声引起的孤立区域和条状区域的边缘不连续像素点，将这些边缘不连续像素点的集合记为mark2；

(3)取集合mark1和集合mark2的并集之后，利用图像膨胀原理，设置一个圆形结构元素，以并集里的所有像素点为圆心进行膨胀操作，最终实现将孤立区域的内部不连续像素点和条状区域附近的不连续像素点都标记出来，将标记出来的所有不连续像素点形成的不连续区域记为mark；

(4)按照以上方法分别求出长短基线各自的移除参考地形相位后绝对相位的不连续区域，分别记为mark_long和mark_short。

所述步骤五中，一种基于长短基线相位连续特性的自适应融合相位估计方法具体实现过程为：

(1)利用步骤四所述的长短基线各自的移除参考地形相位后绝对相位的不连续区域mark_long和mark_short，可以将长短基线各自的移除参考地形相位后绝对相位中的每个像素点标记为连续像素点或不连续像素点。

(2)因为长短基线各自的移除参考地形相位后绝对相位中像素点的个数是相同的，因此可以根据相同坐标位置的像素点在长基线和短基线中的连续特性分成四种情况，分别是：第一种情况，相同坐标位置的像素点在长基线中是不连续像素点，在短基线中是连续像素点；第二种情况，相同坐标位置的像素点在长基线中是不连续像素点，在短基线中也是不连续像素点；第三种情况，相同坐标位置的像素点在长基线中是连续像素点，在短基线中也是连续像素点；第四种情况，相同坐标位置的像素点在长基线中是连续像素点，在短基线中是不连续像素点；

(3)根据所述四种情况，自适应采取对应的融合方法，两种融合方法具体实现过程为：

(4)对于第一种情况，即相同坐标位置的像素点在长基线中是不连续像素点，在短基线中是连续像素点，将满足第一种情况的像素点所对应的短基线移除参考地形相位后的绝对相位数值乘以长基线与短基线的长度比例值，得到的绝对相位替换满足第一种情况的像素点所对应的长基线移除参考地形相位后原本的绝对相位；

(5)对于第二、三、四种情况，利用步骤三中所述的长短基线各自移除参考地形相位后残余相位的相干系数和长基线与短基线的长度比例值，得到长短基线各自移除参考地形相位后残余相位的归一化概率密度函数，最终得到长短基线的联合概率密度函数，采用最大似然法对所述联合概率密度函数进行峰值搜索，搜索到的峰值为最大似然估计值，将所述最大似然估计值替换满足第二、三、四种情况的像素点所对应的长基线移除参考地形相位后原本的绝对相位；

(6)采用以上两种融合方法即实现长基线移除参考地形相位后相位精度更高的绝对融合相位估计。

本发明与现有技术相比的优点在于：

(1)适用性强。不管是平坦地形还是复杂地形，本发明提出的移除参考地形相位方法只需要干涉图同一方位向上近距端和远距端的两个控制点，利用控制点坐标和SAR主辅图像对应的天线位置坐标求出两个控制点的斜距差即可拟合出整幅图像的线性参考地形相位，将干涉相位移除参考地形相位即可实现移除参考地形相位。该方法只需要两个控制点即实现干涉相位的移除参考地形相位，降低了控制点数据获取的难度，适用性较强。

(2)运算效率快。与传统的多基线最大似然相位解缠算法相比，本发明首次提出了将相位导数方差、跳变点判别以及图像膨胀原理结合起来的不连续区域判别方法，在融合处理前长短基线各自移除参考地形相位后的绝对相位进行不连续区域的检测标记，从而为下一步长基线融合相位估计提供预处理数据，为长基线相位的每一个像素点选取不同的融合方法提供了依据，降低了运算的复杂度，有效的提高了算法的运算效率。

(3)相位估计精度高。传统的多基线最大似然相位解缠算法是对多基线的干涉相位进行最大似然估计，从而得到解缠相位，其相位精度受到噪声的影响较大并且运算速度较慢，而本发明中首先通过移除参考地形相位减小了干涉条纹密度，有效地抑制了噪声引起的解缠错误，显著减少了残差点的个数以及基线比例误差对融合相位估计精度的影响，提升了相位估计精度；另一方面，针对单基线InSAR系统下容易解缠失败的复杂陡峭地形，本发明充分利用长短基线能够获取同一场景互补信息的优势，根据长短基线域移除参考地形相位后绝对相位中每一个像素点的连续特性和基于相位导数方差的不连续区域判别方法，自适应选取有效的长短基线观测信息，对长基线移除参考相位后的绝对相位中每一个像素点采取合适的数据融合处理方法，解决了长基线测高精度高但是相位解缠误差较大的问题，与原始的长基线绝对相位相比，长基线融合绝对相位的跳变点个数显著减小，估计精度和鲁棒性均得到了很大程度的提升。

附图说明

图1是本发明的方法流程图；

图2是本发明中的多基线InSAR空间几何关系；

图3是本发明实施示例中生成的长短基线SAR主图像；

图4是本发明实施示例中生成的长基线主辅图像干涉相位；

图5是本发明实施示例中生成的短基线主辅图像干涉相位；

图6是本发明实施示例中生成的长基线参考地形相位；

图7是本发明实施示例中生成的短基线参考地形相位；

图8是本发明实施示例中生成的长基线移除参考地形后的残余缠绕相位；

图9是本发明实施示例中生成的短基线移除参考地形后的残余缠绕相位；

图10是本发明实施示例中生成的长基线滤波之后的移除参考地形相位缠绕相位；

图11是本发明实施示例中生成的短基线滤波之后的移除参考地形相位缠绕相位；

图12是本发明实施示例中生成的长基线移除参考地形相位绝对相位；

图13是本发明实施示例中生成的短基线移除参考地形相位绝对相位；

图14是利用本发明方法得到的长基线域移除参考地形相位融合绝对相位。

具体实施方式

下面将结合附图和实施案例对本发明做进一步的详细说明。

如图1所示，本发明是一种基于相位导数方差的InSAR长短基线融合相位估计方法，具体实现包括以下几个步骤：

步骤一：利用长基线和短基线的SAR主辅图像分别进行图像配准，得到同一场景下长短基线各自的SAR主辅图像干涉相位；

对于InSAR系统，平台在飞行过程中容易受到大气湍流等复杂环境的影响，从而导致天线的运动轨迹会产生横滚、俯仰、偏航等姿态变化，因此该场景下的同一目标点在SAR主辅图像中会产生方位向或距离向的偏移伸缩等现象，如果直接将主辅图像进行共轭相乘得到干涉相位，由于相干性较差很有可能无法获得具有一定规律的干涉条纹，甚至没有干涉条纹出现，因此对SAR主辅图像进行图像配准是InSAR数据处理中的关键步骤。

通常可以从配准精度的角度将图像配准分为粗配准(像素级配准)和精配准(亚像素级配准)。一般来说，多基线InSAR处理中为得到清晰的干涉条纹，所需的图像配准精度需要达到0.1-0.01个像素级。

本发明中首先对SAR主辅图像进行全局相干粗配准处理，粗配准时选择较为稳健的实相关函数法，粗配准采取的实相关函数的数学原理如下式：

其中：ρ_r表示实相关函数，u₁和u₂分别表示SAR主图像和SAR辅图像，(m+x,m+y)表示滑动窗口中心点坐标，(m,n)为像素中心点的坐标，M,N是图像的尺寸大小；|·|表示取幅度值。通过计算SAR主辅图像之间的全局相关函数，再根据相关函数的峰值位置，确定辅图像的像素级偏移量，从而得到粗配准后的SAR辅图像。

粗配准精度无法满足本发明中的处理要求，因此利用相干系数法进行精配准使得干涉条纹更加清晰，精配准的原理为在粗配准后将辅图像进行分块处理，然后对每个小块插值计算实复相干系数，获取相关函数峰值位置，从而得到每个数据块的精配准偏移量，再选用精确的插值核函数分块插值，拼接所有数据块，最终获取精配准后的SAR辅图像。

精配准采用的复相关函数的数学原理如下式：

其中：ρ_c表示复相关函数，u'₁和u'₂分别表示粗配准之后的主图像和辅图像，(m+x,m+y)表示滑动窗口中心点坐标，(m,n)表示小块图像插值后的像素中心点坐标，M,N是图像的尺寸。

根据图2所示的多基线InSAR空间几何关系，按照表1所示的系统参数仿真得到长短基线各自的SAR主辅图像。其中长短基线的SAR主图像为同一幅图像，如图3所示。从图中可以看到，主辅图像对应的场景地形特点是左侧区域为平坦地形，右侧区域的地形起伏变化较大，因此数据处理难度较大。

表1多基线InSAR系统仿真参数

SAR主辅图像经图像配准之后，将不同基线下精配准后的主辅图像进行共轭相乘处理，再求取其结果的相位值，即可得到长短基线的原始干涉相位。图4和图5分别为长短基线的SAR主辅图像配准之后共轭相乘取辐角的干涉相位。对比两幅图可以看出，图4长基线干涉相位图中的干涉条纹较密集，因此相位解缠的难度较大，而图5短基线干涉相位图中的干涉条纹比图4稀疏，但是细节信息没有长基线丰富，这与理论分析一致。

步骤二：选择步骤一中所述场景下同一方位向上近距端和远距端的两个控制点，分别计算所述两个控制点的坐标和长短基线各自的SAR辅图像对应的天线位置坐标之间的距离，即两个控制点的斜距，再计算出所述两个控制点的斜距差，最终拟合出长短基线各自的参考地形相位；

InSAR系统特有的空间几何关系决定了主辅图像共轭相乘得到的干涉相位是目标高程变化和水平斜距变化这两部分相位的叠加。平地效应是指在没有高程变化的平地上也会产生明暗相间的干涉条纹，通常参考地形相位大于反映高程变化的相位，这样会导致反映地形变化的干涉条纹被淹没在较密集的干涉条纹中，从而增加相位解缠的难度和高程反演的精度，特别是在长基线下平地效应的影响更为明显。因此，在对干涉条纹图进行相位解缠之前，通常会将参考地形相位从干涉相位中移除，获得只包含目标高程变化的地形缠绕相位后，这样有利于提升相位解缠的精度。

传统的移除参考地形相位方法有基于干涉图条纹频率的频移法和基于系统轨道参数和先验DEM的直接法。参考地形相位通常是沿距离向明暗相间的周期性变化条纹，变换到频域表现为干涉图的频谱出现偏移，因此频移法的基本原理是找出频谱中的最大值对应的频率并搬移到零频处，从而实现移除参考地形相位。频移法对条纹频率估计精度的要求较高，特别是在复杂地形情况下干涉条纹较密集，加大了移除参考地形相位的难度，运算速度较慢；直接法移除参考地形相位的基本原理为选择一个参考平面，根据天线位置和先验DEM坐标计算出在该坐标系下场景内每一点的斜距差，进而得到参考地形相位。直接法对数据的要求较高，通常需要精确的轨道数据和先验DEM信息，因此在缺乏足够的地形数据和轨道信息的情况下，很难准确估计的参考地形相位。

在本发明中提供了一种新的移除参考地形相位方法，其原理为选择步骤一中所述场景下同一方位向上近距端和远距端的两个控制点，根据所述两个控制点的坐标和长短基线各自的主辅图像对应的天线位置坐标，分别计算出两个控制点与天线之间的距离，即两个控制点的斜距，再计算出所述两个控制点的斜距差，然后从InSAR空间几何关系出发，假定在该场景距离向范围内，主辅星的斜距差变化是线性的，那么则可以利用近距端和远距端两个控制点求解线性关系表达式，最终拟合出长短基线下整幅图像沿距离向的线性参考地形相位。图6和图7分别是长短基线域利用近距端和远距端两个控制点拟合得到的线性参考地形相位。从两幅图中可以看出，线性参考地形相位是沿距离向的明暗相间的竖直条纹，其中长基线线性参考相位的条纹比短基线的条纹要密集一些，这与长短基线的干涉相位呈现出来的特性一致。

本发明提出的移除参考地形相位方法不需要先验DEM信息，也不需要进行精确的条纹频率估计，只需要两个控制点，从而大大提升了运算效率并且降低了对原始数据的要求。虽然拟合得到的参考地形与平地不一定为平行关系，但是由于长短基线域移除的是相同两个控制点拟合得到的长短基线域各自的参考地形相位，这两个参考地形相位比例值与长短基线比例相等，因此不会影响长短基线域移除参考地形相位后的绝对相位比例与长短基线比例的关系，从而可以降低长短基线融合时移除参考地形相位后绝对相位的数量级，减小误差的传递。

步骤三：将长短基线各自的SAR主辅图像干涉相位均减去长短基线各自的参考地形相位，剩下的相位即为长短基线各自移除参考地形相位后的残余相位，然后计算长短基线各自移除参考地形相位后残余相位的相干系数，再依次对长短基线各自移除参考地形相位后的残余相位进行干涉图滤波、相位解缠以及控制点校正，最终得到长短基线各自移除参考地形相位后的绝对相位；

得到长短基线域的线性参考地形相位之后，将长短基线的干涉相位减去各自的参考地形相位，剩下的残余相位是由目标高程变化引起的。图8和图9分别展示了移除参考地形相位后，再将剩下的相位缠绕至[-π,π]的残余相位，从图中可以看出，移除参考地形相位之后，条纹密度显著减小，残余相位可以大致反映出与主图像相对应地形细节。然后针对该残余相位进行干涉图滤波，在InSAR的数据处理中，相位噪声来源于多种因素，具体包括系统热噪声、平台抖动引起的基线去相干、时间去相干、大气延迟、干涉处理过程中引入的误差等。因此为了提高相位解缠的精度，干涉图滤波是至关重要的步骤，经典的相位滤波方法包括均值滤波、圆周中值滤波、Goldstein滤波、Lee滤波、斜坡自适应滤波等，在本发明中，采用的是基于深度学习的干涉图滤波方法，这样可以在有效去除相位噪声的同时保留条纹细节，进一步减小残差点，降低相位解缠的难度，有利于提高后续的相位估计精度，最终滤波之后的长短基线移除参考地形相位后的残余相位分别如图10和图11所示，与滤波之前的残余相位图相比，可以明显的看出滤波可以有效的抑制噪声。

在InSAR数据处理方法中，相位解缠是InSAR高程反演处理过程中的关键问题以及难点，它的准确性直接影响到InSAR生成数字高程模型的精确性。干涉相位是缠绕在[-π,π]之间，而相位解缠的目的则是将InSAR测量中模糊的干涉相位恢复到模糊之前的真实相位，从而保证能够用来正确的反演地面高程。目前已基于不同类型的相位解缠算法开展了研究，大致可以分为以下几类：基于路径跟踪的相位解缠算法、最小二乘法、基于网络规划的相位解缠算法。在本发明中采用的是质量图法与最小费用流结合的相位解缠方法，该方法可以最大限度的减小解缠误差的传递，将有用相位信息进行正确的解缠恢复，并且运算效率较高。

加上参考地形相位后的解缠相位与真实地形相位之间存在一个整体的偏差，需要根据已知的控制点信息求出真实相位与解缠相位之间的差值，对整幅图像利用控制点校正至绝对相位，最终得到移除参考地形相位后的绝对相位。长短基线移除参考地形相位之后的绝对相位分别如图12和图13所示，从图中可以看出虽然短基线的解缠失败区域很少，但是绝对相位的范围较小，难以精确地展现出该区域的地形变化情况，绝对相位的不连续区域较少，而长基线的绝对相位范围较大，能够较好地展现出地形起伏变化细节，但同时不连续区域也显著增加，这两幅图可以很明显的展现出长短基线相位解缠结果的优缺点，验证了理论分析的正确性，因此可以利用长短基线获取的互补信息实现长基线域融合绝对相位估计。

步骤四：针对所述长短基线各自移除参考地形相位后的绝对相位，分别求出长短基线各自移除参考地形相位后的绝对相位归一化的相位导数方差，将所述相位导数方差作为质量图引导，结合绝对相位上的跳变点和图像膨胀原理，标记出长短基线各自移除参考地形相位后的绝对相位的不连续区域；

针对一幅相位图，具有多种质量图评价方法，通常可选择相干系数图、伪相干系数图、相位导数方差图、最大相位梯度图等作为相位质量图。针对真实地形中高程急速变化导致相位不连续等现象，相位导数方差图可以较好的反映出相位变化情况，标记出不连续区域，其定义式为：

其中

为缠绕相位两个方向上的梯度，

和

分别为在k×k范围内

和

的平均值。为了更好的选取相位导数方差阈值，在本发明中将其归一化到0～1之间，它是相位图质量的逆向反映，即归一化相位导数方差数值越高，其真实地形高程变化越剧烈，通常在该区域容易引起解缠失败。

跳变点是指相位梯度差绝对值超过π的像素点，通常可以将跳变点认为是解缠失败的像素点。然而在更多的情况下，相位噪声并不是一个孤立的像素点，而是一个封闭式的孤立区域或者具有一定面积的条状区域，如果用传统的跳变点判别方法，只会将孤立区域或条状区域的边缘检测出来，而孤立区域的内部和渐变区域则没有被标记为不连续点，从而使得在利用长短基线相位信息融合时这些误差仍然存在，影响融合结果。

因此在本发明中，首次提出一种基于相位导数方差、跳变点判别以及图像膨胀原理的不连续区域判别方法，准确标记了移除参考地形相位后绝对相位中的不连续区域，从而为后续的长短基线融合相位估计提供了数据预处理。其基本思路为：首先设置一个0～1之间的相位导数方差阈值，将大于该阈值的像素点标记为不连续点，这样可以将大部分相位梯度变化剧烈的像素点mark1检测出来，然后根据跳变点定义可以标记出噪声引起的孤立区域和条状区域的边缘像素点mark2，取mark1和mark2的并集之后，利用图像膨胀原理，设置一个圆形结构元素，对以上像素点进行膨胀操作，最终可以实现将孤立区域的内部像素点和条状区域的不连续点都标记出来，得到长短基线域移除参考地形相位后的绝对相位不连续区域mark_long和mark_short，从而为下一步长短基线融合相位估计提供预处理数据和不同的相位融合方法选取依据。

步骤五：根据步骤四所述的长短基线各自移除参考地形相位后的绝对相位的不连续区域和步骤三中所述的相干系数，针对长短基线长短基线各自移除参考地形相位后的绝对相位的连续特性，自适应选取融合相位估计方法，最终得到相位精度更高的长基线移除参考相位后的绝对融合相位。

多基线最大似然相位解缠的基本原理是对利用不同长度基线提取的干涉相位进行数据融合，获取更为丰富的目标信息，从而提高解缠精度获得更为准确的真实相位。

在本发明中提出了一种基于长短基线相位连续特性的自适应融合相位估计方法，根据相同坐标位置的像素点在长基线和短基线中的连续特性分成四种情况，分别是：第一种情况，相同坐标位置的像素点在长基线中是不连续像素点，在短基线中是连续像素点；第二种情况，相同坐标位置的像素点在长基线中是不连续像素点，在短基线中也是不连续像素点；第三种情况，相同坐标位置的像素点在长基线中是连续像素点，在短基线中也是连续像素点；第四种情况，相同坐标位置的像素点在长基线中是连续像素点，在短基线中是不连续像素点；根据所述四种情况，自适应采取对应的融合方法，两种融合方法具体实现过程为：

对于第一种情况，即相同坐标位置的像素点在长基线中是不连续像素点，在短基线中是连续像素点，将满足第一种情况的像素点所对应的短基线移除参考地形相位后的绝对相位数值乘以长基线与短基线的长度比例值，得到的绝对相位替换满足第一种情况的像素点所对应的长基线移除参考地形相位后原本的绝对相位；

对于第二、三、四种情况，利用步骤三中所述的长短基线各自移除参考地形相位后残余相位的相干系数和长基线与短基线的长度比例值，得到长短基线各自移除参考地形相位后残余相位的归一化概率密度函数，最终得到长短基线的联合概率密度函数，采用最大似然法对所述联合概率密度函数进行峰值搜索，搜索到的峰值为最大似然估计值，将所述最大似然估计值替换满足第二、三、四种情况的像素点所对应的长基线移除参考地形相位后原本的绝对相位；

利用本发明提出的一种基于相位导数方差的长短基线融合相位估计方法，最终得到长基线移除参考地形相位后的融合绝对相位如图14所示，从图中可以看出，融合绝对相位与融合之前的长基线绝对相位相比，解缠失败减少的同时保持了地形起伏变化的细节信息，验证了本发明所提出方法的有效性。

下面分别从跳变点个数、相位误差均值、相位相对误差标准差三个方面对融合相位估计结果进行评估，并与原始的长基线域绝对相位进行比较，统计结果如表2所示。

表2原始长基线域绝对相位和融合估计相位的精度对比

从表2所示的统计结果可以看出，本发明所提出的移除参考地形相位方法可以快速有效地移除参考地形相位，减小基线比例误差对融合相位估计精度的影响；基于相位导数方差的不连续区域判别方法可以实现长短基线域的不连续区域的准确标记，为后续的融合方法选取提供了依据；根据长短基线绝对相位每一个像素点的连续特性选取对应的融合方法，最终融合估计绝对相位的跳变点个数显著减少，由融合之前的2763减小至150，相位误差的均值和标准差均有不同程度的减小，更加接近于真实相位，显著提升了长基线融合相位估计精度和鲁棒性。

提供以上实施例仅仅是为了描述本发明的目的，而并非要限制本发明的范围。本发明的范围由所附权利要求限定。不脱离本发明的精神和原理而做出的各种等同替换和修改，均应涵盖在本发明的范围之内。

Claims (2)

1.一种基于相位导数方差的InSAR长短基线融合相位估计方法，其特征在于：所述方法包括以下几个步骤：

步骤一：利用长基线和短基线的SAR主辅图像分别进行图像配准，得到同一场景下长短基线各自的SAR主辅图像干涉相位；

步骤二：采用移除参考地形相位方法：选择所述同一场景下的方位向上近距端和远距端的两个控制点，分别计算所述两个控制点的坐标和长短基线各自的SAR主 辅图像对应的天线位置坐标之间的距离，即两个控制点的斜距，再计算出所述两个控制点的斜距差，最终拟合出长短基线各自的参考地形相位；

步骤三：将长短基线各自的SAR主辅图像干涉相位均减去长短基线各自的参考地形相位，剩下的相位即为长短基线各自移除参考地形相位后的残余相位，然后计算长短基线各自移除参考地形相位后残余相位的相干系数，再依次对长短基线各自移除参考地形相位后的残余相位进行干涉图滤波、相位解缠以及控制点校正，最终得到长短基线各自移除参考地形相位后的绝对相位；

步骤四：采用基于相位导数方差、跳变点判别以及图像膨胀原理的不连续区域判别方法：针对所述长短基线各自移除参考地形相位后的绝对相位，分别求出长短基线各自移除参考地形相位后的绝对相位归一化的相位导数方差，将所述相位导数方差作为质量图引导，结合绝对相位上的跳变点和图像膨胀原理，标记出长短基线各自移除参考地形相位后的绝对相位的不连续区域；

步骤五：采用基于长短基线相位连续特性的自适应融合相位估计方法：根据步骤四所述的长短基线各自移除参考地形相位后的绝对相位的不连续区域和步骤三中所述的相干系数，针对长短基线长短基线各自移除参考地形相位后的绝对相位的连续特性，自适应选取融合相位估计方法，最终得到相位精度更高的长基线移除参考相位后的绝对融合相位；

所述步骤二中，移除参考地形相位方法具体实现过程为：

(1)选择所述同一场景下的方位向上近距端和远距端的两个控制点，分别计算所述两个控制点坐标和长短基线各自的SAR主 辅图像对应的天线位置坐标之间的距离，即两个控制点的斜距，再计算出所述两个控制点的斜距差；

(2)从InSAR空间几何关系出发，假定在SAR主辅图像对应场景的同一个方位向上，SAR主辅图像的斜距差沿距离向的变化是线性的，则利用所述两个控制点的斜距差和SAR主辅图像同一方位向上像素点的总个数即可求解出参考地形相位线性关系表达式y＝ax+b中的参数a和b；

(3)最终的参考地形相位通过线性关系表达式y＝ax+b计算得到，其中y代表SAR主辅图像同一方位向上每一个像素点的参考地形相位，x代表每一个像素点在同一个方位向上的顺序；

(4)按照以上方法分别求出长短基线各自的参考地形相位，将长短基线各自的SAR主辅图像干涉相位均减去长短基线各自的参考地形相位，实现了长短基线各自的移除参考地形相位；

所述步骤四中，基于相位导数方差、跳变点判别以及图像膨胀原理的不连续区域判别方法具体实现过程为：

(1)计算移除参考地形相位后绝对相位的归一化相位导数方差，设置一个0～1之间的相位导数方差阈值，将移除参考地形相位后绝对相位中归一化相位导数方差大于相位导数方差阈值的像素点标记为不连续像素点，不连续像素点代表相位梯度变化剧烈的像素点，将这些不连续像素点的集合记为mark1；

(2)根据跳变点的判别方法，将移除参考地形相位后绝对相位中与相邻像素点的相位梯度差绝对值超过π的像素点判定为跳变点，跳变点标记出噪声引起的孤立区域和条状区域的边缘不连续像素点，将这些边缘不连续像素点的集合记为mark2；

(3)取集合mark1和集合mark2的并集之后，利用图像膨胀原理，设置一个圆形结构元素，以并集里的所有像素点为圆心进行膨胀操作，最终实现将孤立区域的内部不连续像素点和条状区域附近的不连续像素点都标记出来，将标记出来的所有不连续像素点形成的不连续区域记为mark；

(4)按照所述方法分别求出长短基线各自的移除参考地形相位后绝对相位的不连续区域。

2.根据权利要求1所述的基于相位导数方差的InSAR长短基线融合相位估计方法，其特征在于：所述步骤五中，基于长短基线相位连续特性的自适应融合相位估计方法具体实现过程为：

(1)利用步骤四所述的长短基线各自的移除参考地形相位后绝对相位的不连续区域，将长短基线各自的移除参考地形相位后绝对相位中的每个像素点标记为连续像素点或不连续像素点；

(2)因为长短基线各自的移除参考地形相位后绝对相位中像素点的个数是相同的，因此可以根据相同坐标位置的像素点在长基线和短基线中的连续特性分成四种情况，分别是：第一种情况，相同坐标位置的像素点在长基线中是不连续像素点，在短基线中是连续像素点；第二种情况，相同坐标位置的像素点在长基线中是不连续像素点，在短基线中也是不连续像素点；第三种情况，相同坐标位置的像素点在长基线中是连续像素点，在短基线中也是连续像素点；第四种情况，相同坐标位置的像素点在长基线中是连续像素点，在短基线中是不连续像素点；

(3)根据所述四种情况，自适应采取对应的融合方法，两种融合方法具体实现过程为：

(4)对于第一种情况，即相同坐标位置的像素点在长基线中是不连续像素点，在短基线中是连续像素点，将满足第一种情况的像素点所对应的短基线移除参考地形相位后的绝对相位数值乘以长基线与短基线的长度比例值，得到的绝对相位替换满足第一种情况的像素点所对应的长基线移除参考地形相位后原本的绝对相位；

(5)对于第二、三、四种情况，利用步骤三中所述的长短基线各自移除参考地形相位后残余相位的相干系数和长基线与短基线的长度比例值，得到长短基线各自移除参考地形相位后残余相位的归一化概率密度函数，最终得到长短基线的联合概率密度函数，采用最大似然法对所述联合概率密度函数进行峰值搜索，搜索到的峰值为最大似然估计值，将所述最大似然估计值替换满足第二、三、四种情况的像素点所对应的长基线移除参考地形相位后原本的绝对相位；

(6)采用以上两种融合方法即实现长基线移除参考地形相位后相位精度更高的绝对融合相位估计。

Images