CN110988876A - 闭合鲁棒双基线InSAR相位解缠方法、系统及可读存储介质 - Google Patents

闭合鲁棒双基线InSAR相位解缠方法、系统及可读存储介质 Download PDF

Info

Publication number
CN110988876A
CN110988876A CN201911147113.3A CN201911147113A CN110988876A CN 110988876 A CN110988876 A CN 110988876A CN 201911147113 A CN201911147113 A CN 201911147113A CN 110988876 A CN110988876 A CN 110988876A
Authority
CN
China
Prior art keywords
pixel
class
phase
intercept
fuzzy
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN201911147113.3A
Other languages
English (en)
Other versions
CN110988876B (zh
Inventor
袁志辉
尹英豪
邢学敏
陈立福
谷翼策
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Changsha University of Science and Technology
Original Assignee
Changsha University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Changsha University of Science and Technology filed Critical Changsha University of Science and Technology
Priority to CN201911147113.3A priority Critical patent/CN110988876B/zh
Publication of CN110988876A publication Critical patent/CN110988876A/zh
Application granted granted Critical
Publication of CN110988876B publication Critical patent/CN110988876B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01S—RADIO DIRECTION-FINDING; RADIO NAVIGATION; DETERMINING DISTANCE OR VELOCITY BY USE OF RADIO WAVES; LOCATING OR PRESENCE-DETECTING BY USE OF THE REFLECTION OR RERADIATION OF RADIO WAVES; ANALOGOUS ARRANGEMENTS USING OTHER WAVES
    • G01S13/00—Systems using the reflection or reradiation of radio waves, e.g. radar systems; Analogous systems using reflection or reradiation of waves whose nature or wavelength is irrelevant or unspecified
    • G01S13/88—Radar or analogous systems specially adapted for specific applications
    • G01S13/89—Radar or analogous systems specially adapted for specific applications for mapping or imaging
    • G01S13/90—Radar or analogous systems specially adapted for specific applications for mapping or imaging using synthetic aperture techniques, e.g. synthetic aperture radar [SAR] techniques
    • G01S13/9021—SAR image post-processing techniques
    • G01S13/9023—SAR image post-processing techniques combined with interferometric techniques
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F18/00—Pattern recognition
    • G06F18/20—Analysing
    • G06F18/23—Clustering techniques

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Physics & Mathematics (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • Electromagnetism (AREA)
  • Theoretical Computer Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • General Engineering & Computer Science (AREA)
  • Evolutionary Computation (AREA)
  • Evolutionary Biology (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Artificial Intelligence (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Image Processing (AREA)

Abstract

本发明公开了一种闭合鲁棒双基线InSAR相位解缠方法、系统及可读存储介质，该方法首先根据双基线InSAR系统参数得到各干涉图对应的高度模糊数，并分解为公因子和整数相乘的形式；接着把所有像素划分为不同的类；然后利用闭合求解公式计算每个类对应的类模糊矢量，进而得到各像素对应的像素模糊矢量；再利用所提出的类相位滤波方法将带有噪声的缠绕干涉相位矢量点投影到其对应的类线段上进行滤波；再根据像素模糊矢量和滤波后的缠绕干涉相位得到解缠后的绝对干涉相位；最后通过绝对干涉相位和其他双基线InSAR系统参数重建地形高程。既能解决相位不连续问题又能提高高程重建精度，还可以提高双基线相位解缠算法的鲁棒性和效率。

Description

闭合鲁棒双基线InSAR相位解缠方法、系统及可读存储介质
技术领域
本发明属于干涉合成孔径雷达领域，特别涉及一种闭合鲁棒双基线InSAR相位解缠方法、 系统和可读存储介质。
背景技术
干涉合成孔径雷达(Interferometric Synthetic Aperture Radar,InSAR)测量技术是以合成孔径 雷达复数据提取的相位信息为信息源获取地表的三维信息和变化信息的一项技术。它可以在 全天时、全天候条件下获取宽测绘范围的地面三维高分辨率图像，在地形制图、地表形变场 探测、速度场探测、森林制图和土地利用分类等方面有着非常广泛的应用。而多基线InSAR 技术可以获得比传统单基线InSAR技术更为精确的地面信息，并能有效克服或减少由目标高 度的急剧变化、较大噪声干扰及叠掩等因素带来的不利影响，提高InSAR系统对复杂地形的 测量精度，是近年来的研究热点，也是未来的发展方向。
根据对缠绕干涉相位处理方式的不同，可将现有的多基线InSAR相位解缠方法分为两大 类：统计估计法及非统计估计法。统计估计法把像素对应的高程或相邻像素间的高程差视为 统计分布框架中的一个待估参数，该框架根据干涉相位的概率密度函数对高程与缠绕干涉相 位之间的关系进行建模，并通过最大似然(ML)准则或最大后验概率(MAP)准则对其进 行估计，从而最终实现相位解缠。而非统计估计法则直接利用具有不同基线长度的多个InSAR 干涉图的组合信息来估计绝对干涉相位，而不需要干涉相位的概率密度函数。
发明内容
本发明的目的在于提供一种闭合鲁棒双基线InSAR相位解缠方法、系统及可读存储介质， 能改善基于聚类分析的双基线InSAR相位解缠算法的效率，增强该算法的鲁棒性，提高高程 重建精度。
一种闭合鲁棒的双基线InSAR相位解缠方法，包括：
根据双基线InSAR系统参数求出各干涉图对应的高度模糊数，并将它们分解为公因子和 整数相乘的形式，所述整数被称为各干涉图对应的分解整数；
根据各干涉图对应的分解整数及由双基线InSAR系统测得的带噪声缠绕干涉相位矢量求 出各像素对应的像素截距；
根据各像素对应的像素截距对所有像素进行聚类，从而获得像素的总类数、各像素所属 的类和每个类对应的类截距；
利用闭合求解公式计算所述每个类对应的类模糊矢量，进而得到各像素对应的像素模糊 矢量；
采用类相位滤波法，将各像素对应的带噪声缠绕干涉相位矢量点沿投影直线投影到该像 素对应的类线段上进行滤波，获得滤波后的缠绕干涉相位矢量点；所述投影直线是过带噪声 缠绕干涉相位矢量点的一条直线，其方向向量由各干涉图对应的相干系数决定；所述类线段 是指同类像素在以两缠绕干涉相位为坐标轴所构成的平面中对应的一条线段；
根据各像素对应的像素模糊矢量和滤波后的缠绕干涉相位得到解缠后的绝对干涉相位， 并依据解缠后的绝对干涉相位和双基线InSAR系统的其它参数重建地形高程。
整个方法求解简单，且解决相位不连续问题又能提高高程重建精度，还可以提高双基线 相位解缠算法的鲁棒性和效率。
进一步地，当各干涉图对应的高度模糊数均为整数时，公因子为所有高度模糊数的最大 公约数；当各干涉图对应的高度模糊数不全为整数时，首先将所有高度模糊数保留相同的n 位小数点位数，接着将所有高度模糊数与10n相乘，然后求取最大公约数，最后将所得的最 大公约数与10n相除，从而得到公因子。各干涉图对应的分解整数为各干涉图对应的高度模 糊数除以公因子后所得的整数。
进一步地，所述采用类相位滤波法，将每个像素对应的带噪声缠绕干涉相位矢量点沿投 影直线投影到该像素对应的类线段上进行滤波，获得滤波后的缠绕干涉相位矢量点的具体过 程如下：
把
Figure BDA0002282516060000021
平面上通过像素s对应的带噪声缠绕干涉相位矢量点
Figure BDA0002282516060000022
且方向向量为 [-1/|γ1|，1/||γ2|]的投影直线与该像素对应的类线段的交点作为滤波后的缠绕干涉相位矢量点
Figure BDA0002282516060000023
类线段所在的直线经过一个与类截距interceptc(l)有关的特定点 [0，-2π·interceptc(l)]，且其方向向量为[1/Γ1，1/Γ2]；γ1和γ2分别为两干涉图对应的相干系数， Γ1和Γ2分别表示两干涉图对应的分解整数；
从而可获得以下两条点向式直线方程：
Figure BDA0002282516060000024
经过变换后可得：
Figure BDA0002282516060000031
其中，
Figure BDA0002282516060000032
为自变量，
Figure BDA0002282516060000033
为因变量，分别对应
Figure BDA0002282516060000034
平面上的横坐标轴和纵坐标轴，
Figure BDA0002282516060000035
为像素s对应的带噪声缠绕干涉相位矢量点，
Figure BDA0002282516060000036
和
Figure BDA0002282516060000037
分别为像素s在两干 涉图上的带噪声缠绕干涉相位；interceptc(l)为像素s对应的类截距，像素s假设属于第l类；
求解上述方程组便可获得滤波后的缠绕干涉相位矢量点
Figure BDA0002282516060000038
其表达式如下：
Figure BDA0002282516060000039
进一步地，当|γ1|/||γ2|＝Γ2/Γ1时，采用垂直投影；当|γ1|/|γ2|＝0时，采用水平投影；当|γ1|/||γ2|＝∞ |γ1|/|γ2|＝∞时，采用竖直投影。
进一步地，所述利用闭合求解公式计算每个类对应的类模糊矢量的过程如下：
首先，构建每个类对应的中心点余数矢量计算方程；
Figure BDA00022825160600000310
其中，
Figure BDA00022825160600000311
和
Figure BDA00022825160600000312
分别表示第l类像素对应的类中心点缠绕干涉相位矢量
Figure BDA00022825160600000313
中的元素；
Figure BDA00022825160600000314
L表示像素的总类数；
然后，根据Γ1、Γ2、qc1(l)和qc2(l)构建同余方程组：
Figure BDA00022825160600000315
其中，kc1(l)和kc2(l)是待求的类模糊矢量[kc1(l),kc2(l)]中的元素，N是待求的整数解，Γ1和Γ2分别表示两干涉图对应的分解整数；
由中国余数定理可知，该同余方程组的解为：
Figure BDA0002282516060000041
其中，Γ＝Γ1·Γ2、ρi＝Γ/Γi,i＝1,2，
Figure BDA0002282516060000042
是ρi关于模Γi的乘法逆元；
最后，根据公式
Figure BDA0002282516060000043
得到每个类对应的类模糊矢量[kc1(l),kc2(l)]。
进一步地，像素s对应的像素截距采用以下公式计算获得：
Figure BDA0002282516060000044
其中，
Figure BDA0002282516060000045
和
Figure BDA0002282516060000046
分别为像素s在基线B1和B2对应干涉图上的带噪声缠绕干涉相位，Γ1和Γ2分别表示两干涉图对应的分解整数；
进一步地，根据像素截距对所有像素进行聚类，从而获得像素的总类数、各像素所属的 类和每个类对应的类截距的过程如下：
首先根据截距统计直方图包络的峰值个数确定像素的总类数；接着将根据类截距计算公 式算出的所有类截距作为候选类截距；然后将距离各峰值最近的候选类截距作为各峰值所在 类对应的类截距；最后以相邻两类截距的算术平均值作为相邻类之间的分界线，从而可按截 距的大小确定各像素所属的类。
进一步地，将各像素的行、列位置和截距信息作为三维空间中的位置坐标，然后根据各 像素在三维空间中的距离远近，按照设定的距离阈值，对像素进行聚类。
一种闭合鲁棒的双基线InSAR相位解缠系统，包括：
高度模糊数求解及分解模块：根据双基线InSAR系统参数求出各干涉图对应的高度模糊 数，并将两高度模糊数分解为一个公因子和整数相乘的形式；
像素截距计算模块：利用双基线InSAR系统测得的各像素对应的带噪声缠绕干涉相位矢 量，以及各干涉图对应的分解整数，求出各像素对应的像素截距；
像素聚类模块：根据各像素对应的像素截距或各像素对应的行、列位置和像素截距对所 有像素进行聚类，从而获得像素的总类数、各像素所属的类和每个类对应的类截距；
类模糊矢量求解模块：采用一种闭合鲁棒的双基线InSAR相位解缠方法，利用闭合求解 公式计算每个类对应的类模糊矢量，进而得到各像素对应的像素模糊矢量；
类相位滤波模块，采用一种闭合鲁棒的双基线InSAR相位解缠方法，将各像素对应的带 噪声缠绕干涉相位矢量点投影到该像素对应的类线段上进行滤波，获得滤波后的缠绕干涉相 位矢量点；
解缠与重建模块，根据各像素对应的像素模糊矢量和滤波后的缠绕干涉相位得到解缠后 的绝对干涉相位，并依据解缠后的绝对干涉相位和双基线InSAR系统的其它参数重建地形高 程。
一种可读存储介质，包括计算机程序指令，其特征在于，所述计算机程序指令被处理终 端执行时使所述处理终端执行上述一种闭合鲁棒的双基线InSAR相位解缠方法。
在本方案中的投影法具体实现过程如下：通过将带有噪声的缠绕干涉相位投影到其对应 的类线段上来实现一定程度的滤波。最直观的投影方式便是垂直投影，它将带有噪声的缠绕 干涉相位矢量点
Figure BDA0002282516060000051
垂直投影到其对应的类线段上而获得滤波后的缠绕干涉相位矢 量点
Figure BDA0002282516060000052
当然，也可以采用其他投影方式。事实上，我们可以把具有任意斜率的过
Figure BDA0002282516060000053
的直线与其对应类线段的交点作为滤波后的缠绕相位矢量
Figure BDA0002282516060000054
不同的 斜率会导致不同的相位滤波结果。例如，如果
Figure BDA0002282516060000055
远比
Figure BDA0002282516060000056
更准确可靠，则可以选择如图 2(b)的水平投影方式，相反，如果
Figure BDA0002282516060000057
远比
Figure BDA0002282516060000058
更准确，则可以选择如图2(c)所示的竖直投 影，以获得最佳的相位滤波结果。可以根据两幅干涉图之间的相干性来确定该投影线的斜率， 为-|γ1|/|γ2|，其中|·|表示取绝对值。
有益效果
本发明提供了一种闭合鲁棒双基线InSAR相位解缠方法、系统及可读存储介质，该方法 首先根据双基线InSAR系统参数得到各干涉图对应的高度模糊数，并分解为公因子和整数相 乘的形式；接着把所有像素划分为不同的类；然后利用闭合求解公式计算每个类对应的类模 糊矢量，进而得到各像素对应的像素模糊矢量；再利用所提出的类相位滤波方法将带有噪声 的缠绕干涉相位矢量点投影到其对应的类线段上进行滤波；再根据像素模糊矢量和滤波后的 缠绕干涉相位得到解缠后的绝对干涉相位；最后通过绝对干涉相位和其他双基线InSAR系统 参数重建地形高程。本技术方案给出了类模糊矢量的闭合求解公式，可以提高聚类算法的计 算效率，为大规模InSAR数据的处理奠定技术基础。同时，本技术方案还提出了类相位滤波 的方法，可以滤除部分相位噪声的影响，从而提高高程重建精度。另外，本技术方案还可以 消除传统单基线InSAR相位解缠技术中存在的相位连续性假设的限制，提高多基线相位解缠 技术的鲁棒性。
附图说明
图1是B2/B1＝Γ1/Γ2的一般情况下缠绕干涉相位
Figure BDA0002282516060000061
和
Figure BDA0002282516060000062
模糊数k1和k2之间的关系示意 图，其中，(a)为缠绕相位
Figure BDA0002282516060000063
和
Figure BDA0002282516060000064
之间的关系示意图；(b)模糊数k1和k2之间的关系示意图。；
图2是类相位滤波示意图，其中，(a)当|γ1|/|γ2|＝Γ2/Γ1时，采用垂直投影；(b)当|γ1|/||γ2|＝0 时，采用水平投影；(c)当|γ1|/||γ2|＝∞时，采用竖直投影；
图3是当B2/B1＝Γ1/Γ2＝5/3的特定情况下缠绕干涉相位
Figure BDA0002282516060000065
和
Figure BDA0002282516060000066
模糊数k1和k2之间的关 系示意图，其中，(a)缠绕相位
Figure BDA0002282516060000067
和
Figure BDA0002282516060000068
之间的关系示意图；(b)模糊数k1和k2之间的关系示 意图；
图4是当B2/B1＝Γ1/Γ2＝5/3时求出的类截距interceptc(l)和类模糊矢量[kc1(l),kc2(l)]之间 的对应关系；
图5是根据一个简单场景对应的仿真干涉图进行双基线相位解缠获得的结果与对比示意 图，其中，(a)是参考DEM(单位：m)，(b)是短基线对应的仿真干涉图，(c)是长基线对应的仿真干涉图，(d)是通过(b)和(c)的线性组合获得的所有像素对应的截距，(e)是 截距统计直方图，(f)是截距统计直方图包络，(g)是短基线对应的滤波后干涉图，(h)是 长基线对应的滤波后干涉图，(i)是类滤波之前的高程重建结果，(j)是类滤波之后的高程重 建结果，(k)是类滤波之前的高程误差，(l)是类滤波之后的高程误差；
图6是根据美国科罗拉多州某山区的真实DEM对应的仿真干涉图进行双基线相位解缠 获得的结果，其中，(a)是参考DEM(单位：m)，(b)是基线1对应的仿真干涉图，(c)是 基线2对应的仿真干涉图，(d)是基线3对应的仿真干涉，(e)是根据(b)和(c)的线性 组合获得的所有像素对应的截距，(f)是根据(b)和(c)进行双基线相位解缠获得的高程 重建结果，(g)是(f)对应的高程误差，(h)是根据(b)和(s)的线性组合获得的所有像 素对应的截距，(i)是根据(b)和(d)进行双基线相位解缠获得的高程重建结果，(j)是(i) 对应的高程误差；
图7是根据中国陕西省铜川市某山区对应的TanDEM-X系统实测双基线InSAR数据集进 行双基线相位解缠获得的结果与对比示意图，其中，(a)是SRTM参考DEM(单位：m)，(b)是长基线对应的干涉图，(c)是通过短基线对应的干涉图和(b)的线性组合获得的所有像素的截距，(d)是短基线对应的滤波后干涉图，(e)是长基线对应的滤波后干涉图，(f)是类 滤波之前的高程重建结果，(g)是(f)对应的高程误差，(h)是类滤波之后的高程重建结果 对应的高程误差。
具体实施方式
下面将结合附图和实例对本发明做进一步的说明。
一种闭合鲁棒双基线InSAR相位解缠方法，包括以下步骤：
(1)根据双基线InSAR系统参数得到两干涉图对应的高度模糊数：
Figure BDA0002282516060000071
其中，B1和B2分别为两干涉图对应的有效基线长度，λ为雷达工作波长，r是主天 线对应的景中心斜距，θ为主天线的视角。并将H1和H2进行分解，得到公因子M和两 干涉图对应的分解整数Γ1和Γ2。
(2)根据InSAR系统测得的每个像素s对应的带噪声缠绕干涉相位矢量
Figure BDA0002282516060000072
求 出每个像素s对应的像素截距：
Figure BDA0002282516060000073
其中，
Figure BDA0002282516060000074
和
Figure BDA0002282516060000075
分别为基线B1和B2对应干涉图上的带噪声缠绕干涉相位。
(3)利用传统的聚类算法根据各像素对应的行、列和截距等信息对所有像素进行聚类，从 而获得像素的总类数L、各像素所属的类l和每个类对应的类截距interceptc(l)(见图1 (b))，其中，l＝1，2，...，L。首先根据截距统计直方图包络的峰值个数确定像素的总 类数；接着将根据类截距计算公式算出的所有类截距作为候选类截距；然后将距离各峰值最近的候选类截距作为各峰值所在类对应的类截距；最后以相邻两类截距的算术平均值作为相邻类之间的分界线，从而可按截距的大小确定各像素所属的类。
当然，在聚类时，也可以采用将各像素的行、列位置和截距信息作为三维空间中的位置坐标，然后根据各像素在三维空间中的距离远近，按照设定的距离阈值，对像素进 行聚类。
(4)根据公式
Figure BDA0002282516060000081
计算每个类对应的类中心点缠绕干涉相位矢量
Figure BDA0002282516060000082
(5)根据公式
Figure BDA0002282516060000083
计算每个类对应的中心点余数矢量[qc1(l)，qc2(l)]。
(6)根据Γ1、Γ2、qc1(l)和qc2(l)构建同余方程组：
Figure BDA0002282516060000084
其中，kc1(l)和kc2(l)是待求的类模糊矢量[kc1(l),kc2(l)]中的元素，N是待求的整数解。
由中国余数定理可知，该同余方程组的解为：
Figure BDA0002282516060000085
其中，Γ＝Γ1·Γ2，ρi＝Γ/Γi,i＝1,2，
Figure BDA0002282516060000086
是ρi关于模Γi的乘法逆元。
(7)根据公式
Figure BDA0002282516060000087
得到每个类对应的类模糊矢量[kc1(l),kc2(l)]。
(8)根据步骤(3)获得的聚类结果计算每个像素对应的像素模糊矢量[k1(s),k2(s)]，即如果 像素s属于第l类，那么就令[k1(s),k2(s)]＝[kc1(l),kc2(l)]。
(9)利用类相位滤波对各像素对应的缠绕干涉相位矢量进行滤波，具体方法如下：
把
Figure BDA0002282516060000088
平面上通过像素s对应的带噪声缠绕干涉相位矢量点
Figure BDA0002282516060000089
且方向向量 为[-1/||γ1|，1/||γ2|]的投影直线与该像素对应的类线段的交点作为滤波后的缠绕干涉相位矢 量点
Figure BDA00022825160600000810
类线段所在的直线经过一个与类截距interceptc(l)有关的特定点[0，-2π·interceptc(l)]，且其方向向量为[1/Γ1，1/Γ2](见图1(a))；γ1和γ2分别为两干涉 图对应的相干系数，Γ1和Γ2分别表示两干涉图对应的分解整数。
从而可获得以下两条点向式直线方程：
Figure BDA0002282516060000091
经过变换后可得：
Figure BDA0002282516060000092
其中，
Figure BDA0002282516060000093
为自变量，
Figure BDA0002282516060000094
为因变量，分别对应
Figure BDA0002282516060000095
平面上的横坐标轴和纵坐标轴，
Figure BDA0002282516060000096
为像素s对应的带噪声缠绕干涉相位矢量点，
Figure BDA0002282516060000097
和
Figure BDA0002282516060000098
分别为像素s在 两干涉图上的带噪声缠绕干涉相位；interceptc(l)为像素s对应的类截距，像素s假设属 于第l类；
求解上述方程组便可获得滤波后的缠绕干涉相位矢量点
Figure BDA0002282516060000099
其表达式如下：
Figure BDA00022825160600000910
如果利用上式求出的滤波后的缠绕干涉相位超过了[0，2π)的取值范围，则应将其缠 绕回[0，2π)范围内。
(10)根据公式
Figure BDA00022825160600000911
求出解缠后的绝对干涉相位，其中ψ1(s)和ψ1(s)分别为两干涉图对应的绝对干涉相 位。
(11)根据公式求出每个像素s对应的高程：
Figure BDA0002282516060000101
利用三个针对仿真和实测InSAR数据的实验来测试本发明提出的闭合鲁棒双基线InSAR 相位解缠方法。
第一个实验是在一个简单的仿真场景下进行的，该场景只有两种高程值(50米和150米)。 图5(a)是仿真的数字高程模型(DEM)，图5(b)和(c)分别是短基线和长基线(300米和500米)对应的仿真干涉图，图5(b)的平均相干系数是0.8，图5(c)的平均相干系数 是0.7，其相应的高度模糊数分别为H1＝73.0米和H2＝43.8米，从而根据步骤(1)有M＝gcd(H1*10, H2*10)/10＝14.6，Γ1＝5，Γ2＝3；根据步骤(2)可以得到如图5(d)所示的截距分布图，对 该分布图进行统计，从而可以得到图5(e)和(f)所示的截距统计直方图及其包络；根据步 骤(3)，由图5(f)可知一共有两个峰值，因此总类数为L＝2，根据类截距计算公式可知候选类截距为图4所示集合S中的元素，由于两个峰值的位置分别为-0.3和1.05，因而每个类对应的类截距分别为interceptc(l＝1)＝-1/3和interceptc(l＝2)＝1。从而可将截距取值在小于 (-1/3+1)/2＝1/3范围内的像素划分为第一类，将截距取值在大于等于1/3范围内的像素划分为 第二类；根据步骤(4)有
Figure BDA0002282516060000102
根据步 骤(5)，我们可以得到[q1(1)，q2(1)]＝[1，2]，[q1(2)，q2(2)]＝[3,0]；根据步骤(6)和(7)，我们 可以求出[k1(1),k2(1)]＝[2，3]，[k1(2),k2(2)]＝[0，1]；根据步骤(8)得到每个像素所属的类；根 据步骤(9)得到如图5(g)和(h)所示的完成类相位滤波后的短基线和长基线对应的滤波 干涉图，与图5(b)和(c)相比，图5(g)和(h)的相位噪声明显减弱，从而验证了类相 位滤波的有效性。根据步骤(10)和(11)得到如图5(j)所示的类相位滤波之后的DEM， 图5(i)为类相位滤波之前的DEM，实际上，图5(i)就是利用传统的基于聚类分析的双基 线相位解缠方法获得的高程重建结果，图5(j)是本发明提出的基于聚类分析的双基线相位 解缠及滤波方法获得的高程重建结果，两种方法的结果对比见表1。
在第二个实验中，我们根据美国科罗拉多州某山区的真实DEM对应的仿真干涉图进行 双基线相位解缠，获得了满意的实验结果。图6(a)为参考DEM，对应的最小和最大高程分别为0米和136.7米。图6(b),(c)和(d)分别是通过仿真得到的三个不同基线长度(60 米,200米和320米)对应的带噪声干涉图，相应的高度模糊数分别设定为93.0米，27.9米和17.4米，图6(b)和(c)形成第一组基线组合(B2/B1＝10/3)，图6(b)和(d)形成第二组 基线组合(B2/B1＝16/3)。三幅干涉图对应的相干系数都设定为0.9，以便比较不同基线组合在 相同相位噪声条件下分别进行双基线相位解缠得到的高程重建结果。图6(e)是通过图6(b)和(c)的线性组合获得的截距，利用本发明所提出的基于聚类分析的双基线相位解缠及滤波 方法获得的高程重建结果如图6(f)所示，图6(g)是对应的高程误差。图6(h)是通过图 6(b)和(d)的线性组合获得的截距，图6(i)是利用本发明所提出的闭合鲁棒的双基线相 位解缠方法获得的高程重建结果，图6(j)对应的高程误差。两种基线组合对应的实验结果 见表1。
表1三个实验对应的误差指标
Figure BDA0002282516060000111
在第三个实验中，我们用实测单航过TanDEM-X双基线InSAR数据集(736像素×191像 素)来评估本发明所提方法的性能。该数据集对应于中国陕西省铜川市的一个山区，相关参 数如表2所示。
表2
TANDEM-X双基线INSAR数据集的主要参数
Figure BDA0002282516060000112
图7(a)为航天飞机雷达地形测绘任务(SRTM)对应的参考DEM，作为真实高程。图 7(b)表示长基线对应的单航过InSAR干涉图，其基线长度为361.90米。图7(c)为通过 短基线对应的单航过InSAR干涉图和(b)的线性组合获得的截距。图7(d)和(e)分别为 类相位滤波后对应的两干涉图。图7(f)为类相位滤波之前的高程重建结果，对应传统基于 聚类分析的双基线相位解缠方法的高程重建结果。图7(g)为传统基于聚类分析的双基线相 位解缠方法对应的高程误差。图7(h)为本发明所提方法对应的高程误差。实验结果见表1： 图7(g)的平均高程误差为4.16m，标准偏差为7.48m，归一化重建平方误差为5.528×10^-5； 图7(h)的平均高程误差为3.87米，标准偏差为4.93米，归一化重建平方误差为4.290×10^-5。 因此，实验结果验证了本发明所提基于聚类分析的闭合鲁棒双基线相位解缠及滤波方法对真实InSAR数据处理的有效性。
一种闭合鲁棒的双基线InSAR相位解缠系统，包括：
高度模糊数求解及分解模块：根据双基线InSAR系统参数求出各干涉图对应的高度模糊 数，并将两高度模糊数分解为一个公因子和整数相乘的形式；
像素截距计算模块：利用双基线InSAR系统测得的各像素对应的带噪声缠绕干涉相位矢 量，以及各干涉图对应的分解整数，求出各像素对应的像素截距；
像素聚类模块：根据各像素对应的像素截距或各像素对应的行、列位置和像素截距对所 有像素进行聚类，从而获得像素的总类数、各像素所属的类和每个类对应的类截距；
类模糊矢量求解模块：采用一种闭合鲁棒的双基线InSAR相位解缠方法，利用闭合求解 公式计算每个类对应的类模糊矢量，进而得到各像素对应的像素模糊矢量；
类相位滤波模块，采用一种闭合鲁棒的双基线InSAR相位解缠方法，将各像素对应的带 噪声缠绕干涉相位矢量点投影到该像素对应的类线段上进行滤波，获得滤波后的缠绕干涉相 位矢量点；
解缠与重建模块，根据各像素对应的像素模糊矢量和滤波后的缠绕干涉相位得到解缠后 的绝对干涉相位，并依据解缠后的绝对干涉相位和双基线InSAR系统的其它参数重建地形高 程。
应当理解，本发明各个实施例中的功能单元模块可以集中在一个处理单元中，也可以是 各个单元模块单独物理存在，也可以是两个或两个以上的单元模块集成在一个单元模块中， 可以采用硬件或软件的形式来实现。
本发明实施例还提供一种可读存储介质，包括计算机程序指令，其特征在于，所述计算 机程序指令被处理终端执行时使所述处理终端执行上述一种闭合鲁棒的双基线InSAR相位解 缠方法。其有益效果参见方法部分的有益效果，在此不再赘述。
本领域内的技术人员应明白，本申请的实施例可提供为方法、系统、或计算机程序产 品。因此，本申请可采用完全硬件实施例、完全软件实施例、或结合软件和硬件方面的实施 例的形式。而且，本申请可采用在一个或多个其中包含有计算机可用程序代码的计算机可用 存储介质(包括但不限于磁盘存储器、CD-ROM、光学存储器等)上实施的计算机程序产品 的形式。
本申请是参照根据本申请实施例的方法、设备(系统)、和计算机程序产品的流程图和 /或方框图来描述的。应理解可由计算机程序指令实现流程图和/或方框图中的每一流程和 /或方框、以及流程图和/或方框图中的流程和/或方框的结合。可提供这些计算机程序指 令到通用计算机、专用计算机、嵌入式处理机或其他可编程数据处理设备的处理器以产生一 个机器，使得通过计算机或其他可编程数据处理设备的处理器执行的指令产生用于实现在流 程图一个流程或多个流程和/或方框图一个方框或多个方框中指定的功能的装置。
这些计算机程序指令也可存储在能引导计算机或其他可编程数据处理设备以特定方式工 作的计算机可读存储器中，使得存储在该计算机可读存储器中的指令产生包括指令装置的制 造品，该指令装置实现在流程图一个流程或多个流程和/或方框图一个方框或多个方框中指 定的功能。
这些计算机程序指令也可装载到计算机或其他可编程数据处理设备上，使得在计算机或 其他可编程设备上执行一系列操作步骤以产生计算机实现的处理，从而在计算机或其他可编 程设备上执行的指令提供用于实现在流程图一个流程或多个流程和/或方框图一个方框或多 个方框中指定的功能的步骤。
最后应当说明的是：以上实施例仅用以说明本发明的技术方案而非对其限制，尽管参照 上述实施例对本发明进行了详细的说明，所属领域的普通技术人员应当理解：依然可以对本 发明的具体实施方式进行修改或者等同替换，而未脱离本发明精神和范围的任何修改或者等 同替换，其均应涵盖在本发明的权利要求保护范围之内。

Claims (10)

1.一种闭合鲁棒双基线InSAR相位解缠方法，其特征在于，包括：
根据双基线InSAR系统参数求出各干涉图对应的高度模糊数，并将它们分解为公因子和整数相乘的形式，所述整数被称为各干涉图对应的分解整数；
根据各干涉图对应的分解整数及由双基线InSAR系统测得的带噪声缠绕干涉相位矢量求出各像素对应的像素截距；
根据各像素对应的像素截距对所有像素进行聚类，从而获得像素的总类数、各像素所属的类和每个类对应的类截距；
利用闭合求解公式计算所述每个类对应的类模糊矢量，进而得到各像素对应的像素模糊矢量；
采用类相位滤波法，将各像素对应的带噪声缠绕干涉相位矢量点沿投影直线投影到该像素对应的类线段上进行滤波，获得滤波后的缠绕干涉相位矢量点；所述投影直线是过带噪声缠绕干涉相位矢量点的一条直线，其方向向量由各干涉图对应的相干系数决定；所述类线段是指同类像素在以两缠绕干涉相位为坐标轴所构成的平面中对应的一条线段；根据各像素对应的像素模糊矢量和滤波后的缠绕干涉相位得到解缠后的绝对干涉相位，并依据解缠后的绝对干涉相位和双基线InSAR系统的其它参数重建地形高程。
2.根据权利要求1所述的方法，其特征在于，当各干涉图对应的高度模糊数均为整数时，公因子为所有高度模糊数的最大公约数；
当各干涉图对应的高度模糊数不全为整数时，首先将所有高度模糊数保留相同的n位小数点位数，接着将所有高度模糊数与10n相乘，然后求取最大公约数，最后将所得的最大公约数与10n相除，从而得到公因子；各干涉图对应的分解整数为各干涉图对应的高度模糊数除以公因子后所得的整数。
3.根据权利要求1所述的方法，其特征在于，所述类相位滤波法，将各像素对应的带噪声缠绕干涉相位矢量点沿投影直线投影到该像素对应的类线段上进行滤波，获得滤波后的缠绕干涉相位矢量点，具体过程如下：
把
Figure FDA0002282516050000011
平面上通过像素s对应的带噪声缠绕干涉相位矢量点
Figure FDA0002282516050000012
且方向向量为[-1/|γ1|，1/|γ2|]的投影直线与该像素对应的类线段的交点作为滤波后的缠绕干涉相位矢量点
Figure FDA0002282516050000013
类线段所在的直线经过一个与类截距interceptc(l)有关的特定点[0，-2π·interceptc(l)]，且其方向向量为[1/Γ1，1/Γ2]；γ1和γ2分别为两干涉图对应的相干系数；
从而可获得以下两条点向式直线方程：
Figure FDA0002282516050000021
经过变换后可得：
Figure FDA0002282516050000022
其中，
Figure FDA0002282516050000023
为自变量，
Figure FDA0002282516050000024
为因变量，分别对应
Figure FDA0002282516050000025
平面上的横坐标轴和纵坐标轴，
Figure FDA0002282516050000026
为像素s对应的带噪声缠绕干涉相位矢量点，
Figure FDA0002282516050000027
和
Figure FDA0002282516050000028
分别为像素s在两干涉图上的带噪声缠绕干涉相位；interceptc(l)为像素s对应的类截距，像素s假设属于第l类；
求解上述方程组便可获得滤波后的缠绕干涉相位矢量点
Figure FDA0002282516050000029
其表达式如下：
Figure FDA00022825160500000210
其中，Γ1和Γ2分别表示两干涉图对应的分解整数。
4.根据权利要求3所述的方法，其特征在于，当|γ1|/|γ2|＝Γ2/Γ1时，采用垂直投影；当|γ1|/|γ2|＝0时，采用水平投影；当|γ1|/|γ2|＝∞时，采用竖直投影。
5.根据权利要求1所述的方法，其特征在于，所述利用闭合求解公式计算每个类对应的类模糊矢量的过程如下：
首先，构建每个类对应的中心点余数矢量计算方程；
Figure FDA00022825160500000211
其中，
Figure FDA00022825160500000212
和
Figure FDA00022825160500000213
分别表示第l类像素对应的类中心点缠绕干涉相位矢量
Figure FDA00022825160500000214
中的元素；
Figure FDA0002282516050000031
L表示像素的总类数；
然后，根据Γ1、Γ2、qc1(l)和qc2(l)构建同余方程组：
Figure FDA0002282516050000032
其中，kc1(l)和kc2(l)是待求的类模糊矢量[kc1(l),kc2(l)]中的元素，N是待求的整数解，Γ1和Γ2分别表示两干涉图对应的分解整数；
由中国余数定理可知，该同余方程组的解为：
Figure FDA0002282516050000033
其中，Γ＝Γ1·Γ2、ρi＝Γ/Γi,i＝1,2，
Figure FDA0002282516050000034
是ρi关于模Γi的乘法逆元；
最后，根据公式
Figure FDA0002282516050000035
得到每个类对应的类模糊矢量[kc1(l),kc2(l)]。
6.根据权利要求1所述的方法，其特征在于，所述像素s对应的像素截距采用以下公式计算获得：
Figure FDA0002282516050000036
其中，
Figure FDA0002282516050000037
和
Figure FDA0002282516050000038
分别为像素s在基线B1和B2对应干涉图上的带噪声缠绕干涉相位，Γ1和Γ2分别表示两干涉图对应的分解整数。
7.根据权利要求1-6任一项所述的方法，其特征在于，所述根据像素截距对所有像素进行聚类，从而获得像素的总类数、各像素所属的类和每个类对应的类截距的过程如下：
首先根据截距统计直方图包络的峰值个数确定像素的总类数；接着将根据类截距计算公式算出的所有类截距作为候选类截距；然后将距离各峰值最近的候选类截距作为各峰值所在类对应的类截距；最后以相邻两类截距的算术平均值作为相邻类之间的分界线，从而可按截距的大小确定各像素所属的类。
8.根据权利要求1-5任一项所述的方法，其特征在于，将各像素的行、列位置和截距信息作为三维空间中的位置坐标，然后根据各像素在三维空间中的距离远近，按照设定的距离阈值，对像素进行聚类。
9.一种闭合鲁棒双基线InSAR相位解缠系统，其特征在于，包括：
高度模糊数求解及分解模块：根据双基线InSAR系统参数求出各干涉图对应的高度模糊数，并将两高度模糊数分解为一个公因子和整数相乘的形式；
像素截距计算模块：利用双基线InSAR系统测得的各像素对应的带噪声缠绕干涉相位矢量，以及各干涉图对应的分解整数，求出各像素对应的像素截距；
像素聚类模块：根据各像素对应的像素截距或各像素对应的行、列位置和像素截距对所有像素进行聚类，从而获得像素的总类数、各像素所属的类和每个类对应的类截距；
类模糊矢量求解模块：采用权利要求1-8任一项所述的方法，利用闭合求解公式计算每个类对应的类模糊矢量，进而得到各像素对应的像素模糊矢量；
类相位滤波模块，采用权利要求1-8任一项所述的方法，将各像素对应的带噪声缠绕干涉相位矢量点投影到该像素对应的类线段上进行滤波，获得滤波后的缠绕干涉相位矢量点；
解缠与重建模块，根据各像素对应的像素模糊矢量和滤波后的缠绕干涉相位得到解缠后的绝对干涉相位，并依据解缠后的绝对干涉相位和其他双基线InSAR系统参数重建地形高程。
10.一种可读存储介质，包括计算机程序指令，其特征在于，所述计算机程序指令被处理终端执行时使所述处理终端执行权利要求1至8任一项所述的方法。
CN201911147113.3A 2019-11-21 2019-11-21 闭合鲁棒双基线InSAR相位解缠方法、系统及可读存储介质 Active CN110988876B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911147113.3A CN110988876B (zh) 2019-11-21 2019-11-21 闭合鲁棒双基线InSAR相位解缠方法、系统及可读存储介质

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911147113.3A CN110988876B (zh) 2019-11-21 2019-11-21 闭合鲁棒双基线InSAR相位解缠方法、系统及可读存储介质

Publications (2)

Publication Number Publication Date
CN110988876A true CN110988876A (zh) 2020-04-10
CN110988876B CN110988876B (zh) 2023-01-03

Family

ID=70085526

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911147113.3A Active CN110988876B (zh) 2019-11-21 2019-11-21 闭合鲁棒双基线InSAR相位解缠方法、系统及可读存储介质

Country Status (1)

Country Link
CN (1) CN110988876B (zh)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113885027A (zh) * 2021-09-17 2022-01-04 长沙理工大学 用于多基线干涉合成孔径雷达相位解缠的聚类校正方法、系统和介质
CN115015919A (zh) * 2022-02-16 2022-09-06 华北水利水电大学 双基线InSAR割平面纯整数规划相位解缠方法
CN117872366A (zh) * 2023-11-30 2024-04-12 中国科学院空天信息创新研究院 一种基于相位解缠的阵列干涉sar山区点云解模糊方法

Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101551455A (zh) * 2009-05-13 2009-10-07 西安电子科技大学 干涉合成孔径雷达三维地形成像系统及其高程测绘方法
CN103323844A (zh) * 2013-04-22 2013-09-25 中国科学院电子学研究所 一种多通道干涉合成孔径雷达高程重建方法及装置
CN103983973A (zh) * 2014-05-28 2014-08-13 西安电子科技大学 基于图像稀疏域噪声分布约束的合成孔径雷达成像方法
CN106199601A (zh) * 2016-07-01 2016-12-07 西安电子科技大学 基于粗数字高程图的InSAR绝对相位模糊估计方法
CN106249236A (zh) * 2016-07-12 2016-12-21 北京航空航天大学 一种星载InSAR长短基线图像联合配准方法
CN108008382A (zh) * 2017-10-30 2018-05-08 西安空间无线电技术研究所 一种多基地星载干涉sar系统测量陡峭地形的方法
CN108663678A (zh) * 2018-01-29 2018-10-16 西北农林科技大学 基于混合整数优化模型的多基线InSAR相位解缠算法
JPWO2018123748A1 (ja) * 2016-12-27 2019-10-31 日本電気株式会社 画像解析装置、画像解析方法及びプログラム

Patent Citations (8)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN101551455A (zh) * 2009-05-13 2009-10-07 西安电子科技大学 干涉合成孔径雷达三维地形成像系统及其高程测绘方法
CN103323844A (zh) * 2013-04-22 2013-09-25 中国科学院电子学研究所 一种多通道干涉合成孔径雷达高程重建方法及装置
CN103983973A (zh) * 2014-05-28 2014-08-13 西安电子科技大学 基于图像稀疏域噪声分布约束的合成孔径雷达成像方法
CN106199601A (zh) * 2016-07-01 2016-12-07 西安电子科技大学 基于粗数字高程图的InSAR绝对相位模糊估计方法
CN106249236A (zh) * 2016-07-12 2016-12-21 北京航空航天大学 一种星载InSAR长短基线图像联合配准方法
JPWO2018123748A1 (ja) * 2016-12-27 2019-10-31 日本電気株式会社 画像解析装置、画像解析方法及びプログラム
CN108008382A (zh) * 2017-10-30 2018-05-08 西安空间无线电技术研究所 一种多基地星载干涉sar系统测量陡峭地形的方法
CN108663678A (zh) * 2018-01-29 2018-10-16 西北农林科技大学 基于混合整数优化模型的多基线InSAR相位解缠算法

Non-Patent Citations (5)

* Cited by examiner, † Cited by third party
Title
HANWEN YU 等: "A cluster-analysis-based efficient multibaseline phase-unwrapping algorithm", 《IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING》 *
ZHIHUI YUAN 等: "Multichannel InSAR DEM Reconstruction Through Improved Closed-Form Robust Chinese Remainder Theorem", 《IEEE GEOSCIENCE AND REMOTE SENSING LETTERS》 *
刘宝泉: "干涉合成孔径雷达测量关键技术研究", 《中国优秀博硕士学位论文全文数据库(博士)信息科技辑》 *
潘舟浩 等: "三基线毫米波InSAR的相位解缠及高程反演", 《红外与毫米波学报》 *
陈立福 等: "双天线干涉SAR系统无控制点场景的高精度参考相位快速估计算法", 《国防科技大学学报》 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113885027A (zh) * 2021-09-17 2022-01-04 长沙理工大学 用于多基线干涉合成孔径雷达相位解缠的聚类校正方法、系统和介质
CN113885027B (zh) * 2021-09-17 2024-10-18 长沙理工大学 用于多基线干涉合成孔径雷达相位解缠的聚类校正方法、系统和介质
CN115015919A (zh) * 2022-02-16 2022-09-06 华北水利水电大学 双基线InSAR割平面纯整数规划相位解缠方法
CN117872366A (zh) * 2023-11-30 2024-04-12 中国科学院空天信息创新研究院 一种基于相位解缠的阵列干涉sar山区点云解模糊方法
CN117872366B (zh) * 2023-11-30 2024-06-07 中国科学院空天信息创新研究院 一种基于相位解缠的阵列干涉sar山区点云解模糊方法

Also Published As

Publication number Publication date
CN110988876B (zh) 2023-01-03

Similar Documents

Publication Publication Date Title
CN110363858B (zh) 一种三维人脸重建方法及系统
CN109658515B (zh) 点云网格化方法、装置、设备及计算机存储介质
CN110988876B (zh) 闭合鲁棒双基线InSAR相位解缠方法、系统及可读存储介质
KR101475382B1 (ko) 광학적 3차원 측량의 자기 적응 윈도우 푸리에 위상추출방법
CN103871039B (zh) 一种sar图像变化检测差异图生成方法
US11074752B2 (en) Methods, devices and computer program products for gradient based depth reconstructions with robust statistics
CN111275633A (zh) 基于图像分割的点云去噪方法、系统、装置和存储介质
CN108122280A (zh) 一种三维点云的重建方法及装置
CN111524173B (zh) 一种基于双参考平面的快速大范围相位解包裹方法
CN104318552B (zh) 基于凸包投影图匹配的模型配准方法
CN117274515A (zh) 基于ORB和NeRF映射的视觉SLAM方法及系统
CN116563493A (zh) 基于三维重建的模型训练方法、三维重建方法及装置
CN104299241A (zh) 基于 Hadoop 的遥感图像显著性目标检测方法及系统
CN115205360A (zh) 复合条纹投影钢管的三维外轮廓在线测量与缺陷检测方法及应用
CN110688440B (zh) 一种适用于子地图重叠部分较少的地图融合方法
Faion et al. Recursive Bayesian pose and shape estimation of 3D objects using transformed plane curves
CN111583323B (zh) 一种单帧结构光场三维成像方法和系统
CN117132737A (zh) 一种三维建筑模型构建方法、系统及设备
CN112802084A (zh) 基于深度学习的三维形貌测量方法、系统和存储介质
CN116958416A (zh) 一种三维建模方法、装置、系统以及存储介质
CN108680119B (zh) 一种分区域的单幅快速相位展开方法
CN107741204B (zh) 一种用于动态三维测量的条纹增强方法
CN113409457B (zh) 立体图像的三维重构与可视化方法及设备
CN115239559A (zh) 一种融合视图合成的深度图超分辨率方法及系统
CN114820921A (zh) 动态物体的三维成像方法、装置、设备及存储介质

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant