用于多基线干涉合成孔径雷达相位解缠的聚类校正方法、系统和介质
阅读说明：本技术 用于多基线干涉合成孔径雷达相位解缠的聚类校正方法、系统和介质 (Cluster correction method, system and medium for multi-baseline interferometric synthetic aperture radar phase unwrapping ) 是由 袁志辉 陈天骄 徐海胜 彭葳 陈立福 邢学敏 于 2021-09-17 设计创作，主要内容包括：本发明公开了用于多基线干涉合成孔径雷达相位解缠的聚类校正方法、系统和介质,方法为：根据多基线干涉合成孔径雷达的参数,计算每幅干涉图中每个像素对应的整周模糊数,从而得到每个像素的模糊矢量；将模糊矢量相同的像素聚为同一类；对干涉图中的像素进行类校正：若干涉图尺寸小于预设尺寸值,则对所有像素,均按照像素的扩展区域内的最大类别进行类校正；若干涉图尺寸超过预设尺寸值,则先根据像素的扩展区域的密度判断该像素是否需要校正,然后对需要校正的像素,按照像素的扩展区域内的最大类别进行类校正。本发明可以提高现有多基线干涉合成孔径雷达相位解缠的像素聚类结果的准确性和有效性,进而提高解缠相位的精度。(The invention discloses a clustering correction method, a system and a medium for multi-baseline interferometric synthetic aperture radar phase unwrapping, wherein the method comprises the following steps: calculating the fuzzy number of the whole circle corresponding to each pixel in each interference pattern according to the parameters of the multi-baseline interference synthetic aperture radar, thereby obtaining the fuzzy vector of each pixel; pixels with the same fuzzy vector are gathered into the same type; carrying out class correction on pixels in the interference pattern: if the sizes of the plurality of charts are smaller than the preset size value, performing class correction on all pixels according to the maximum class in the expansion area of the pixels; and if the sizes of the plurality of images exceed the preset size value, judging whether the pixel needs to be corrected according to the density of the expansion area of the pixel, and then performing class correction on the pixel needing to be corrected according to the maximum class in the expansion area of the pixel. The method can improve the accuracy and effectiveness of the pixel clustering result of phase unwrapping of the existing multi-baseline interferometric synthetic aperture radar, and further improve the accuracy of unwrapping phases.)
技术领域
本发明属于干涉合成孔径雷达技术领域，具体涉及一种用于多基线干涉合成孔径雷达相位解缠的聚类校正方法、系统和介质。
背景技术
干涉合成孔径雷达(InSAR)可根据同一场景对应的两幅SAR图像的绝对干涉相位获取地面目标的高程信息，最终可用于获取数字高程模型(DEM)[1][2]。然而，InSAR系统只能得到绝对相位的主值，即缠绕相位。为了精确估计地形高度，必须通过相位解缠(PU)来重建绝对相位[3][4]。为了解决这个问题，研究人员首先提出了单基线相位解缠方法(SBPU)[4]。对于传统的SBPU，一个方程中有两个未知数，这将导致许多解。为了保证解的唯一性，SBPU假设相邻像素之间的实际相位跳变小于π，这被称为相位连续性假设或Itoh条件[5]。为了克服SBPU的缺点，有研究者提出了以中国余数定理(Chinese remainder theorem,CRT)为理论基础的多基线相位解缠算法(MBPU)，以消除相位连续性假设的限制[6]。
在过去的几十年里，MBPU得到了广泛的研究。文献[6]首次提出了三种使用具有不同基线或频率的多干涉图MBPU方法。然后在文献[7]和[8]中提出了基于最大似然估计的方法。随后，在[9]和[10]中提出了一种基于贝叶斯框架和最大后验概率(MAP)估计的新策略。在[11]中，将最大似然准则和扩展卡尔曼滤波器(EKF)相结合，构建了多通道EKF相位解缠框架，实现了有趣的性能改进。在[12]中，非局部滤波技术被引入到基于MAP的MBPU方法中。在[13]中，提出了一种两阶段规划方法(TSPA)，它将SBPU的框架移植到MBPU。
作为最流行的MBPU方法之一，基于聚类分析的MBPU算法也得到了广泛的研究。在[14]中，提出了一种基于聚类分析的快速方法，简称CA。聚类方法根据多个干涉图的组合信息将所有像素划分为不同的类，然后逐个进行PU。在[15]中，每个像素对应的模糊向量根据闭合形式的鲁棒CRT求解，然后根据模糊向量进行聚类，并设计了一些措施来提高性能。在[16]中，提出了一种鲁棒的基于CA的多基线干涉图PU算法(缩写为CANOPUS)，其可识别的数学模型从截距维度扩展到行、列和截距维度。在[17]中，采用线性组合方法来增加高度模糊数并提高CA方法的噪声鲁棒性。文献[18]考虑了基线过多或搜索空间过大的情况，给出了类模糊矢量的闭合解公式，同时还提出了一种新的相位滤波策略，进一步提高了解缠相位的精度。
但这种基于聚类分析的MBPU方法仍然存在鲁棒性较弱的问题，容易产生错误的聚类结果:(1)直接基于CRT的PU方法对相位噪声敏感，较小的相位噪声可能造成较大的解缠误差，因此噪声鲁棒性较差；(2)虽然CA方法在一定程度上提高了噪声鲁棒性，减少了PU的时间，但是当含有像素较少的类与含有像素较多的类的截距非常接近时，由于相位噪声的存在，会将这两个类合并在一起，造成错误的聚类结果；(3)CANOPUS法仅将待聚类像素的聚类信息由截距扩展到行、列及其截距值，但仍然容易得到错误的聚类结果。
发明内容
本发明的目的在于提供一种用于多基线干涉合成孔径雷达相位解缠的聚类校正方法、系统和介质，提高基于CA的MBPU算法的噪声鲁棒性。
一种用于多基线干涉合成孔径雷达相位解缠的聚类校正方法，包括：
根据多基线干涉合成孔径雷达的参数，计算每幅干涉图中每个像素对应的整周模糊数；记第i条基线对应的第i幅干涉图中第s个像素对应的整周模糊数为ki(s)，则第s个像素的模糊矢量表示为[k1(s),k2(s),…,kM(s)]，M为基线的数量；
根据各像素对应的模糊矢量，对所有像素进行聚类：模糊矢量相同的像素属于同一类；
对干涉图中的像素进行类校正：若干涉图尺寸小于预设尺寸值，则对所有像素，均按照像素的扩展区域内的最大类别进行类校正；若干涉图尺寸超过预设尺寸值，则先根据像素的扩展区域的密度判断该像素是否需要校正，然后对需要校正的像素，按照像素的扩展区域内的最大类别进行类校正。
进一步地，像素的扩展区域是指：以当前像素为中心向周围扩展到预设尺寸的矩形窗口区域。
进一步地，像素扩展区域内的最大类别是指：在当前像素扩展区域内，包括最多像素的类别；对像素按照像素扩展区域内的最大类别进行类校正具体是指：将当前像素的类别修改为当前像素扩展区域内的最大类别。
进一步地，像素的扩展区域的密度是指：在扩展区域内，与当前像素的类别相同的像素个数；或者，像素的扩展区域的密度是指：在扩展区域内，与当前像素对应的模糊矢量截距之差小于预设截距阈值的所有像素的个数。
进一步地，根据像素的扩展区域的密度判断该像素是否需要校正的方法为：若像素的扩展区域的密度大于预设密度阈值，则当前像素的类别准确，记为核心像素，即为不需要校正的像素；否则，当前像素的类别错误，记为非核心像素，即为需要校正的像素。
进一步地，基线的数量M＝2，则两幅干涉图中各像素对应的整周模糊数之间关系表达式为：
式中，ki(s)是第i幅干涉图中第s个像素对应的整周模糊数，Bi是第i幅干涉图中的垂直基线的长度，是第i幅干涉图中第s个像素的缠绕相位，i＝1,2。
用于多基线干涉合成孔径雷达相位解缠的聚类校正系统，包括存储器及处理器，所述存储器中存储有计算机程序，其特征在于，所述计算机程序被所述处理器执行时，使得所述处理器实现如上述任一项所述的聚类校正方法。
一种可读存储介质，包括计算机程序指令，所述计算机程序指令被处理终端执行时使所述处理终端执行上述任一项所述的聚类校正方法。
有益效果
本发明对现有的多基线干涉合成孔径雷达相位解缠的像素聚类结果进行类别校正：若干涉图尺寸小于预设尺寸值，则对所有像素，均按照像素的扩展区域内的最大类别进行类校正；若干涉图尺寸超过预设尺寸值，则先根据像素的扩展区域的密度判断该像素是否需要校正，然后对需要校正的像素，按照像素的扩展区域内的最大类别进行类校正。本发明可以提高现有多基线干涉合成孔径雷达相位解缠的像素聚类结果的准确性和有效性，进而提高解缠相位的精度。
附图说明
图1是本申请第一实施例的参考EDM和使用长短两条基线对应的仿真干涉图；其中，(a)参考DEM(单位：m)，(b)短基线对应的仿真干涉图，(c)长基线对应的仿真干涉图。
图2是本申请第一实施例使用现有CA法分类、CANOPUS法分类，及使用本申请方法对现有CA法分类、CANOPUS法分类校正后的类别分布图；其中，(a)CA法对应的类别分布图(b)PPCC法校正后的CA法类别分布图(c)NPCC1法校正后的CA法类别分布图(d)NPCC2法校正后的CA法类别分布图(e)CANOPUS法对应的类别分布图(f)PPCC法校正后的CANOPUS法类别分布图(g)NPCC1法校正后的CANOPUS法类别分布图(h)NPCC2法校正后的CANOPUS法类别分布图。
图3是本申请第二实施例的参考EDM和使用长短两条基线对应的仿真干涉图；其中，(a)参考DEM(单位：m)，(b)长基线对应的无噪声仿真干涉图，(c)短基线对应的无噪声仿真干涉图，(d)长基线对应的有噪声仿真干涉图，(e)短基线对应的有噪声仿真干涉图。
图4是本申请第二实施例使用现有CA法分类、CANOPUS法分类，及使用本申请方法对现有CA法分类、CANOPUS法分类校正后的类别分布图；其中，(a)CA法对应的类别分布图，(b)PPCC法校正后的CA法类别分布图，(c)NPCC1法校正后的CA法类别分布图，(d)NPCC2法校正后的CA法类别分布图，(e)CANOPUS法对应的类别分布图，(f)PPCC法校正后的CANOPUS法类别分布图，(g)NPCC1法校正后的CANOPUS法类别分布图，(h)NPCC2法校正后的CANOPUS法类别分布图。
具体实施方式
下面对本发明的实施例作详细说明，本实施例以本发明的技术方案为依据开展，给出了详细的实施方式和具体的操作过程，对本发明的技术方案作进一步解释说明。
本发明提供一种用于多基线干涉合成孔径雷达相位解缠的聚类校正方法，包括以下步骤：
步骤1，根据多基线干涉合成孔径雷达的参数，计算每幅干涉图中每个像素对应的整周模糊数；记第i条基线对应的第i幅干涉图中第s个像素对应的整周模糊数为ki(s)，则第s个像素的模糊矢量表示为[k1(s),k2(s),…,kM(s)]，M为基线的数量。
步骤2，根据各像素对应的模糊矢量，对所有像素进行聚类：模糊矢量相同的像素属于同一类。
下面先以基于聚类分析的多基线相位解缠算法基本原理进行解释：
对于InSAR系统，由该系统获得的目标的测量相位，即缠绕相位，可以由下式(1)确定给出：
其中，是第s个像素的缠绕相位，即绝对相位的主值，ψ(s)是第s个像素的去平绝对干涉相位，k(s)是第s个像素的未知整周模糊数。从(1)中可以看出，相位解缠的目的是找到正确的k(s)的值，以便得到正确的解缠相位ψ(s)。因为在一个方程中存在两个未知数，无法直接由(1)求解两个未知数ψ(s)和k(s)。有两种方法求解(1)。第一种方法是利用传统的SBPU，即根据相位连续性假设，相邻像素之间的绝对干涉相位之差小于π，这是一个病态问题。另一种方法便是MBPU。为简便起见，在以下讨论中只考虑两个不同基线的情况。对于InSAR系统，地形高度与解缠相位的关系为
其中，h(s)是第s个像素的地形高度，λ是雷达的工作波长，r(s)表示第s个像素对应的目标点到主天线之间的斜距，θ是入射角，Bi是第i幅干涉图中的垂直基线的长度，ψi(s)是第s个像素在第i幅干涉图中的去平绝对干涉相位。
结合(1)(2)可得：
由(3)可以看出，在确定基线情况下，该直线的斜率B2/B1是常数，截距则由B1、B2、和共同决定。定义模糊矢量[k1(s),k2(s)]表示第s个像素在两幅干涉图中对应的整周模糊数，同时定义具有相同模糊矢量的像素属于同一个类。具有相同模糊矢量的像素在k1-k2平面上具有相同的斜率B2/B1，且通过相同的整数点(k1,k2)。由此可得，具有相同模糊矢量的像素对应的直线相互重叠，也即它们的截距是相同的。这就是基于聚类分析的多基线相位解缠算法的实质，即所有截距相同的直线对应的像素都被划分到同一个类中。
文献[14]和[16]中介绍的两种经典的基于CA的MBPU算法在一定程度上可以满足PU的要求，但仍存在一些不足。基于聚类分析的多基线相位解缠方法往往存在许多错误的聚类结果，需要对其进行校正。因此，本发明针对干涉图的尺寸提出两种类校正方法，通过以下步骤3，来改善现有技术中聚类结果的准确性。
步骤3，对干涉图中的像素进行类校正：若干涉图尺寸小于预设尺寸值，则对所有像素，均按照像素的扩展区域内的最大类别进行类校正；若干涉图尺寸超过预设尺寸值，则先根据像素的扩展区域的密度判断该像素是否需要校正，然后对需要校正的像素，按照像素的扩展区域内的最大类别进行类校正。
下两针对两种不同尺寸干涉图所采用的两种类校正方法进行解释说明：
A逐像素类校正(Pixel by pixel cluster correction,PPCC)法
根据相干系数等指标选择窗口的大小，将窗口中出现频率最高的类别作为窗口中心像素的类别。对所有像素逐一进行类校正，以确定每个像素对应的最终类别。具体步骤如下:
首先选择一个像素，然后以该像素为中心设计合适的矩形窗口，根据干涉图的大小设定矩形窗口的大小。统计每个矩形窗口中每个类别的出现次数，查找当前窗口中出现次数最多的类别，然后将当前中心像素的类别修改为矩形窗口中出现次数最多的类别。重复该过程，直到干涉图中所有像素的类校正完成。PPCC的伪代码如下:
将待校正像素集合定义为M；
1：WHILE M≠Φ；
2：选择M中的一个像素i，以该像素为中心选定一个合适的窗口；
3：统计窗口内各类别出现的次数；
4：将像素i的类别修改为窗口中出现次数最多的类别；
5：END WHILE
B非核心像素类校正(Non core-pixel cluster correction,NPCC)法
该方法首先利用密度来区分核心像素与非核心像素。其中，密度的定义有两种，第一种定义将选定窗口内与中心像素类别相同的像素个数定义为该中心像素的密度，第二种定义将选定窗口内与中心像素截距之差小于某个固定值的像素个数定义为该中心像素的密度。如果采用第一种定义，假设选取的矩形窗是N×N的矩阵，且有n个像素与中心像素的类别相同，则可将该中心像素的密度记为n。接着，将密度大于设定门限值的像素当作核心像素，认为是具有正确类别的像素；将密度小于设定门限值的像素当作非核心像素，认为是具有错误类别的像素。然后对非核心像素进行类校正。本发明将采用第一种密度定义的类校正方法称为NPCC1法，将采用第二种密度定义的类校正方法称为NPCC2法。下面将给出非核心像素类校正法的伪代码。
将待校正像素集合定义为M；
1：WHILE M≠Φ；
2：选择M中的一个像素i，以该像素为中心选定一个合适的窗口；
3：选择一种密度定义，计算像素i对应的密度n；
4：根据窗口内密度分布情况设定密度门限值nT；
5：若像素i的密度大于nT，则记为核心像素，反之，则记为非核心像素；
6：根据最大似然准则，将非核心像素的类别更改为窗口内出现次数最多的类别；
7：END WHILE
上述第一种是PPCC方法，将窗口中出现次数最多的类别作为中心像素的最终类别。第二种方法是NPCC方法，该方法使用密度来区分核心像素和非核心像素，然后对非核心像素进行类校正。当干涉图尺度较小，噪声较严重时，应采用PPCC方法，这样得到的PU结果具有较高的精度。在干涉图尺度大、噪声低的情况下，应采用NPCC方法，在保证PU结果准确性的同时，提高算法的计算效率。以下通过实验结果来验证本发明所述类校正方法的有效性。
第一个实验建立一个简单的仿真场景，它表现出平滑和不连续的区域。设定该仿真场景的真实高程值分别为35米和80米(具体DEM如图1(a)所示)。设定长、短基线的长度分别为500米和300米，对应干涉图的相干系数分别为0.8和0.7。该仿真场景对应的真实截距(即无噪声情况下对应的截距)分别为1和1/3，因而在无噪声情况下用聚类分析法进行聚类时应该只有2个类。仿真得到的带噪声干涉图分别如图1(b)、(c)所示。
接下来用本发明提出的类校正方法来校正已有聚类算法产生的错误聚类结果，从而验证本发明所述方法的有效性。
用CA法产生的类别分布图如图2(a)所示，可见存在很多错误的聚类结果。PPCC法应用于该图后对应的类别分布图如图2(b)所示。NPCC1法应用于该图后对应的类别分布图如图2(c)所示。NPCC2法应用于该图后对应的类别分布图如图2(d)所示。从图中可以看出，此时大部分错误类别都被校正过来了。
由于在CANOPUS法中不同的类可能具有相同的模糊矢量，因而可以存很多个类。用CANOPUS法产生的类别分布图如图2(e)所示，可见图中存在很多很小的类，有些类甚至只有几个像素，这些类显然是没有意义的。PPCC法应用于该图后对应的类别分布图如图2(f)所示。NPCC1法应用于该图后对应的类别分布图如图2(g)所示。NPCC2法应用于该图后对应的类别分布图如图2(h)所示。此时小类的数量明显减少了很多。
由该校正结果可知，本文提出的两种类校正方法都可以将错误的类别校正过来，从而提高相位解缠的准确率。
下面定量分析两种类校正方法的性能。对于CANOPUS法，由于不能直接根据类别来判断某像素所属类别的正确性，因而接下来本文将通过像素对应模糊矢量的正确性来判断解缠结果的正确性。定义相位解缠成功率(Success rate of phase unwrapping)为干涉图中被正确恢复出模糊数的像素占总像素的百分比。表1给出了实验1中CA法和CANOPUS法对应的相位解缠成功率。由该表可知，通过类校正后，相位解缠成功率有了很大的提升。
表1实验1的解缠成功率
作为第二个实施例，测试了类校正方法在更真实的高度剖面上的性能，图3(a)为美国科罗拉多州某山区(Isolation Peak)对应的DEM。相应的高度模糊数分别设定为32.1米53.5米，仿真得到的两幅无噪声干涉图如图3(b)(c)所示，对应的带噪声干涉图如图3(d)(e)所示。
图4(a)是用CA法产生的类别分布图。图4(b)是利用PPCC法对图4(a)进行校正后对应的类别分布图。图4(c)是利用NPCC1法对图4(a)进行校正后对应的类别分布图。图4(d)是利用NPCC2法对图4(a)进行校正后对应的类别分布图。图4(e)是CANOPUS法产生的类别分布图。图4(f)是利用PPCC法对图4(e)进行校正后对应的类别分布图。图4(g)是利用NPCC1法对图4(e)进行校正后对应的类别分布图。图4(h)是利用NPCC2法对图4(e)进行校正后对应的类别分布图。从图中可知，大部分错误的类别均被校正过来了。
表2给出了实验2中CA法和CANOPUS法对应的相位解缠成功率。由该表可知，通过类校正后，相位解缠成功率依然有较大的提升，可见本文所提方法对真实地形也是有效的。
表2实验2的解缠成功率
虽然两种方法都是依据最大似然准则进行类校正，但是各有优点。PPCC方法是对聚类后的所有像素进行校正，也可以对类别错误的核心像素进行校正，PU成功率高于NPCC方法。然而，当像素数相对较大或相位噪声相对较小时，由于对所有像素进行类校正，PPCC方法比NPCC方法花费更多的时间进行类校正。因此，当像素数较大或噪声影响较小时，NPCC方法比PPCC方法更有效、更节省时间。同时可以达到更好的相位解缠效果。
以上实施例为本申请的优选实施例，本领域的普通技术人员还可以在此基础上进行各种变换或改进，在不脱离本申请总的构思的前提下，这些变换或改进都应当属于本申请要求保护的范围之内。