CN108594222B - 一种双频干涉合成孔径雷达的高程重建方法和装置 - Google Patents

一种双频干涉合成孔径雷达的高程重建方法和装置 Download PDF

Info

Publication number
CN108594222B
CN108594222B CN201810237097.6A CN201810237097A CN108594222B CN 108594222 B CN108594222 B CN 108594222B CN 201810237097 A CN201810237097 A CN 201810237097A CN 108594222 B CN108594222 B CN 108594222B
Authority
CN
China
Prior art keywords
phase
wave band
interference
filtered
diagram
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.)
Active
Application number
CN201810237097.6A
Other languages
English (en)
Other versions
CN108594222A (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.)
Institute of Electronics of CAS
Original Assignee
Institute of Electronics of CAS
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 Institute of Electronics of CAS filed Critical Institute of Electronics of CAS
Priority to CN201810237097.6A priority Critical patent/CN108594222B/zh
Publication of CN108594222A publication Critical patent/CN108594222A/zh
Application granted granted Critical
Publication of CN108594222B publication Critical patent/CN108594222B/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
    • 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/904—SAR modes

Landscapes

  • Engineering & Computer Science (AREA)
  • Remote Sensing (AREA)
  • Radar, Positioning & Navigation (AREA)
  • Physics & Mathematics (AREA)
  • Electromagnetism (AREA)
  • Computer Networks & Wireless Communication (AREA)
  • General Physics & Mathematics (AREA)
  • Radar Systems Or Details Thereof (AREA)

Abstract

本发明实施例公开了一种双频干涉合成孔径雷达的高程重建方法，该方法包括：分别获取第一波段、第二波段的两幅SAR图像，对每个波段的两幅SAR图像进行处理，得到每个波段的去平干涉相位图；对每个波段的两幅SAR图像中的主图像的幅度图和去平干涉相位图分别进行滤波，得到每个波段的滤波后的主图像的幅度图和滤波后的去平干涉相位图；根据这两个波段的滤波后的主图像的幅度图，分别对第一波段、第二波段的滤波后的去平干涉相位图进行配准，得到每个波段的配准后的去平干涉相位图；对这两个波段的配准后的去平干涉相位图进行相位解缠，得到相位解缠结果；再对相位解缠结果进行高程反演，得到数字高程模型；该方法实现了对目标区域的精确的高程反演。

Description

一种双频干涉合成孔径雷达的高程重建方法和装置
技术领域
本发明涉及数据处理技术，尤其涉及一种双频干涉合成孔径雷达的高程重建方法、装置及计算机可读存储介质。
背景技术
干涉合成孔径雷达(Interferometric Synthetic Aperture Radar,InSAR)是应用在地理测量和遥感领域的一种雷达技术，其使用卫星或飞机搭载的合成孔径雷达系统获取高分辨率的目标区域的SAR(Synthetic Aperture Radar，合成孔径雷达)图像，每一分辨元的SAR图像信息中不仅含有灰度信息，而且还包含干涉所需的相位信号。具体地，通过雷达向目标区域发射微波，然后接收该目标区域反射的回波，得到该目标区域成像的SAR图像对，若SAR图像对之间存在相干条件，SAR图像对共轭相乘可以得到干涉图，根据干涉图的相位值，得出两次成像中微波的路程差，从而计算出目标地区的高程信息，建立数字高程模型。
双频干涉合成孔径雷达(Dual-frequency Interferometry Synthetic ApertureRadar，DFInSAR)利用同一场景两个波段的四幅SAR图像，分别得到两个波段中每个波段的干涉图，从两个波段的干涉图中提取相位信息来获取地面的高程信息，DFInSAR技术可以更加有效地提取地形起伏比较大区域的复杂地形高度信息，解决传统InSAR存在的叠掩问题，但是现有的DFInSAR技术存在不同波段的干涉图的配准不够精确的问题，导致高程反演的精度较差。
发明内容
本发明的主要目的在于提出一种双频干涉合成孔径雷达的高程重建方法及装置，旨在解决现有双频干涉合成孔径雷达的高程重建方法中两个波段的干涉相位图的配准不够准确，降低了高程反演精度的问题。
本发明的技术方案是这样实现的：
本发明实施例提供一种双频干涉合成孔径雷达的高程重建方法，所述方法包括：
分别获取第一波段的两幅合成孔径雷达SAR图像和第二波段的两幅SAR图像，对每个波段的两幅SAR图像进行处理，得到每个波段的去平干涉相位图；
对每个波段的两幅SAR图像中的主图像的幅度图和去平干涉相位图分别进行滤波，得到每个波段的滤波后的主图像的幅度图和滤波后的去平干涉相位图；
根据第一波段的滤波后的主图像的幅度图和第二波段的滤波后的主图像的幅度图，分别对所述每个波段的滤波后的去平干涉相位图进行配准，得到每个波段的配准后的去平干涉相位图；
对第一波段的配准后的去平干涉相位图和第二波段的配准后的去平干涉相位图进行相位解缠，得到相位解缠结果；
对所述相位解缠结果进行高程反演，得到数字高程模型。
上述方案中，所述相位解缠结果包括：每个像素点的解缠后的真实相位；
相应地，所述对所述相位解缠结果进行高程反演，得到数字高程模型，包括：
对所述每个像素点的解缠后的真实相位进行加权均值滤波，得到每个像素点的滤波后的真实相位；
对所述每个像素点的滤波后的真实相位进行高程反演，得到所述数字高程模型。
上述方案中，所述相位解缠结果还包括：
每个像素点的解缠后的真实相位对应的相位周期数；
相应地，所述对所述每个像素点的滤波后的真实相位进行加权均值滤波，得到每个像素点的滤波后的真实相位，包括：
根据所有像素点的解缠后的真实相位和解缠后的真实相位对应的相位周期数、以及预设的有效像素点的确定方式，确定出所有像素点中有效的像素点；
根据所有像素点中有效的像素点，对所述每个像素点的解缠后的真实相位进行加权均值滤波，得到所述每个像素点的滤波后的真实相位。
上述方案中，所述对每个波段的两幅SAR图像中的主图像的幅度图和去平干涉相位图分别进行滤波，包括：
采用非局部滤波的干涉合成孔径雷达NL-InSAR方法，对每个波段的两幅SAR图像中的主图像的幅度图和去平干涉相位图分别进行滤波。
上述方案中，所述对每个波段的两幅SAR图像进行处理，得到每个波段的去平干涉相位图，包括：
对所述每个波段的两幅SAR图像依次进行配准，得到每个波段的配准后的两幅SAR图像；
对所述每个波段的配准后的两幅SAR图像进行干涉处理，获得每个波段的配准后的干涉图；
对所述每个波段的配准后的干涉图进行平地效应的去除，得到每个波段的去平干涉相位图。
本发明实施例还提供一种双频干涉合成孔径雷达的高程重建装置，所述装置包括存储器和处理器；其中，
所述存储器，用于存储计算机程序
所述处理器，用于在运行所述计算机程序时，执行以下步骤：
分别获取第一波段的两幅合成孔径雷达SAR图像和第二波段的两幅SAR图像；
对每个波段的两幅SAR图像进行处理，得到每个波段的去平干涉相位图；
对每个波段的两幅SAR图像中的主图像的幅度图和去平干涉相位图分别进行滤波，得到每个波段的滤波后的主图像的幅度图和滤波后的去平干涉相位图；
根据第一波段的滤波后的主图像的幅度图和第二波段的滤波后的主图像的幅度图，分别对所述每个波段的滤波后的去平干涉相位图进行配准，得到每个波段的配准后的去平干涉相位图；
对第一波段的配准后的去平干涉相位图和第二波段的配准后的去平干涉相位图进行相位解缠，得到相位解缠结果；
对所述相位解缠结果进行高程反演，得到数字高程模型。
上述方案中，所述相位解缠结果包括：每个像素点的解缠后的真实相位；
相应地，所述处理器具体用于在运行所述计算机程序时，执行以下步骤：
对所述每个像素点的解缠后的真实相位进行加权均值滤波，得到每个像素点的滤波后的真实相位；
对所述每个像素点的滤波后的真实相位进行高程反演，得到所述数字高程模型。
上述方案中，所述相位解缠结果还包括：每个像素点的解缠后的真实相位对应的相位周期数；
相应地，所述处理器具体用于在运行所述计算机程序时，执行以下步骤：
根据所有像素点的解缠后的真实相位和解缠后的真实相位对应的相位周期数、以及预设的有效像素点的确定方式，确定出所有像素点中有效的像素点；
根据所有像素点中有效的像素点，对所述每个像素点的解缠后的真实相位进行加权均值滤波，得到所述每个像素点的滤波后的真实相位。
上述方案中，所述处理器具体用于在运行所述计算机程序时，执行以下步骤：
采用非局部滤波的干涉合成孔径雷达NL-InSAR方法，对每个波段的两幅SAR图像中的主图像的幅度图和去平干涉相位图分别进行滤波。
本发明实施例还提供一种计算机可读存储介质，所述计算机可读存储介质存储有计算机程序，
当所述计算机程序被至少一个处理器执行时，导致所述至少一个处理器执行上述任意一种双频干涉合成孔径雷达的高程重建方法的步骤。
本发明实施例提供了一种双频干涉合成孔径雷达的高程重建方法，分别获取第一波段的两幅合成孔径雷达SAR图像和第二波段的两幅SAR图像，对每个波段的两幅SAR图像进行处理，得到每个波段的去平干涉相位图；对每个波段的两幅SAR图像中的主图像的幅度图和去平干涉相位图分别进行滤波，得到每个波段的滤波后的主图像的幅度图和滤波后的去平干涉相位图；根据第一波段的滤波后的主图像的幅度图和第二波段的滤波后的主图像的幅度图，分别对所述每个波段的滤波后的去平干涉相位图进行配准，得到每个波段的配准后的去平干涉相位图；对第一波段的配准后的去平干涉相位图和第二波段的配准后的去平干涉相位图进行相位解缠，得到相位解缠结果；对所述相位解缠结果进行高程反演，得到数字高程模型。如此，本发明实施例中每个波段的滤波后的主图像的幅度图和滤波后的去平干涉相位图的信噪比提高，相应地，两个波段的去平干涉相位图的配准精度也被提高，从而实现精确的高程反演。
附图说明
图1为本发明实施例提供的双频干涉合成孔径雷达的高程重建方法流程图一；
图2a为本发明实施例提供的DFInSAR系统获得的C波段的两幅SAR图像中的主图像；
图2b为本发明实施例提供的DFInSAR系统获得的C波段的两幅SAR图像中的辅图像；
图3为本发明实施例提供的C波段的配准前的干涉图；
图4为本发明实施例提供的双频干涉合成孔径雷达的高程重建方法流程图二；
图5为本发明实施例提供的C波段的配准后的干涉图；
图6为本发明实施例提供的DFInSAR系统与目标区域的成像几何关系图；
图7为本发明实施例提供的频移法得到的去平干涉相位图；
图8为本发明实施例提供的基于轨道参数直接计算法得到的去平干涉相位图；
图9a为本发明实施例提供的C波段的两幅SAR图像中的主图像的幅度图；
图9b为本发明实施例提供的采用NL-InSAR方法得到的C波段的滤波后的主图像的幅度图；
图10a为本发明实施例提供的基于轨道参数直接计算法得到的C波段的去平干涉相位图；
图10b为本发明实施例提供的采用NL-InSAR方法得到的C波段的滤波后的去平干涉相位图；
图11a为本发明实施例提供的C波段的两幅SAR图像之间的相干系数图；
图11b为本发明实施例提供的采用NL-InSAR方法得到的C波段的滤波后的两幅SAR图像之间的相干系数图；
图12a为本发明实施例提供的C波段的配准后的去平干涉相位图；
图12b为本发明实施例提供的X波段的配准后的去平干涉相位图；
图13a为本发明实施例提供的C波段和X波段的配准前的复干涉图之间的相干系数图；
图13b为本发明实施例提供的C波段和X波段的配准后的复干涉图之间的相干系数图；
图14a为本发明实施例提供的采用中国剩余定理(Chinese remainder theorem，CRT)方法得到的相位解缠结果示意图；
图14b为本发明实施例提供的采用最大似然估计(Maximum likelihoodestimate，MLE)方法得到的相位解缠结果示意图；
图14c为本发明实施例提供的采用整体变分最大后验估计(Total Variation-Maximum a posteriori estimate，TV-MAP)方法得到的相位解缠结果示意图；
图15a为本发明实施例提供的采用CRT方法得到的残差图；
图15b为本发明实施例提供的采用MLE方法得到的残差图；
图15c为本发明实施例提供的采用TV-MAP方法得到的残差图；
图16a为本发明实施例提供的所有像素点的解缠后的真实相位图；
图16b为本发明实施例提供的所有像素点的滤波后的真实相位图；
图17为本发明实施例提供的某山区场景图；
图18为本发明实施例提供的某山区场景中部分区域的光学图像的灰度图；
图19a为本发明实施例提供的单波段InSAR方法获取的相位解缠结果示意图；
图19b为本发明实施例提供的双波段DFInSAR方法获取的相位解缠结果示意图；
图20为本发明实施例提供的某区域的数字高程模型示意图；
图21为本发明实施例提供的分别采用C波段的InSAR方法、X波段的InSAR方法以及C波段和X波段的DFInSAR方法，获得的目标区域的高程信息示意图；
图22为本发明实施例提供的双频干涉合成孔径雷达的高程重建装置的组成结构示意图。
具体实施方式
下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述。
实施例一
本发明实施例提供一种双频干涉合成孔径雷达的高程重建方法，如图1所示，该方法包括：
S101：分别获取第一波段的两幅SAR图像和第二波段的两幅SAR图像，对每个波段的两幅SAR图像进行处理，得到每个波段的去平干涉相位图。
在实际实施时，卫星或飞机搭载的DFInSAR系统安装有两个天线，能够工作在两个不同的波段，分别是第一波段和第二波段，DFInSAR系统每次对同一个目标区域进行成像，可以获得第一波段的两幅SAR图像和第二波段的两幅SAR图像，SAR图像以复数形式记录了振幅和相位信息，每个波段的两幅SAR图像包括了一个主图像和一个辅图像。
其中，第一波段和第二波段可以是现有无线电波段中任一个，无线电波段包括了C波段、X波段、L波段、S波段、Ku波段、K波段和Ka波段等。
示例性地，第一波段可以是C波段(5.6GHz)，DFInSAR系统通过工作在C波段的双天线对目标区域进行成像，获得C波段的两幅SAR图像中的主图像如图2a所示、C波段的两幅SAR图像中的辅图像如图2b所示，直接对C波段的两幅SAR图像进行干涉处理，得到C波段的配准前的干涉图如图3所示，图3中数值表示相位，可以看到C波段的配准前的干涉图中观测相位为毫无意义的噪声；这是由于获得的每个波段的两幅SAR图像在距离向及方位向都存在一定偏移、错位或旋转，其相干性较差，无法形成可观察到的呈现条纹状的干涉图，也就是说，两幅图像的相干性决定了干涉图中观测相位的精度，而干涉相位的精度又决定了高程反演的精度，因此，首先应该对每个波段的两幅SAR图像进行配准，提高每个波段的两幅SAR图像之间的相干性，再对每个波段的配准后的两幅SAR图像进行干涉处理，才能得到每个波段的呈现条纹状的配准后的干涉图。
进一步地，获得每个波段的配准后的干涉图之后，由于平地效应的存在，使得每个波段的配准后的干涉图中的干涉条纹过于密集，不利于相位解缠，因此，对每个波段的配准后的干涉图进行平地效应的去除，得到每个波段的去平干涉相位图。
进一步地，对每个波段的配准后的两幅SAR图像进行相干系数的计算，得到每个波段的两幅SAR图像之间的相干系数图。
S102：对每个波段的两幅SAR图像中的主图像的幅度图和去平干涉相位图分别进行滤波，得到每个波段的滤波后的主图像的幅度图和滤波后的去平干涉相位图。
在实际实施时，由于不同波段的SAR图像散射特性不一样，通常对两个波段的干涉图的配准采用相干系数法，具体包括：分别对第一波段的两幅SAR图像中的主图像和第二波段的两幅SAR图像中的主图像进行取绝对值的操作，可以得到对应的幅度图，根据第一波段的两幅SAR图像中的主图像的幅度图和第二波段的两幅SAR图像中的主图像的幅度图，利用最小二乘法拟合出偏移系数，进而实现对这两个波段的干涉图的分别配准；其中，这两个波段的主图像的幅度图对这两个波段的干涉图的配准的精度影响较大，因此，在本发明实施例中，为了提高这两个波段的去平干涉相位图的配准精度，进而提高相位解缠的正确性，对每个波段的两幅SAR图像中的主图像的幅度图和去平干涉相位图进行滤波，从而提高每个波段的两幅SAR图像中的主图像的幅度图和去平干涉相位图的信噪比，减小噪声对这两个波段的去平干涉相位图的配准精度的影响。
进一步地，对每个波段的两幅SAR图像之间的相干系数图进行滤波，得到每个波段的滤波后的两幅SAR图像之间的相干系数图。
S103：根据第一波段的滤波后的主图像的幅度图和第二波段的滤波后的主图像的幅度图，分别对每个波段的滤波后的去平干涉相位图进行配准，得到每个波段的配准后的去平干涉相位图。
在实际实施时，在得到第一波段的配准后的去平干涉相位图和第二波段的配准后的去平干涉相位图之后，分别对第一波段的滤波后的两幅SAR图像之间的相干系数图和第二波段的滤波后的两幅SAR图像之间的相干系数图进行配准，得到第一波段的配准后的两幅SAR图像之间的相干系数图和第二波段的配准后的两幅SAR图像之间的相干系数图。
S104：对第一波段的配准后的去平干涉相位图和第二波段的配准后的去平干涉相位图进行相位解缠，得到相位解缠结果。
由于根据每个波段的两幅SAR图像获得的自身的干涉图中每个像素点的观测相位的取值范围是(-π，π)，丢失了2Kπ的相位信息，需要通过相位解缠，给干涉图中每个像素点的观测相位加上自身的正确的2Kπ，得到每个像素点的解缠后的真实相位，其中，K表示每个像素点的解缠后的真实相位对应的相位周期数，取值为1至N，N为正整数。
在实际实施时，采用第一波段的配准后的两幅SAR图像之间的相干系数和第二波段的配准后的两幅SAR图像之间的相干系数，对第一波段的配准后的去平干涉相位图和第二波段的配准后的去平干涉相位图进行相位解缠，得到相位解缠结果。
S105：对相位解缠结果进行高程反演，得到数字高程模型。
需要说明的是，步骤S101至步骤S105均可以由处理器实现。
由此可见，本发明实施例中，通过对每个波段的两幅SAR图像中的主图像的幅度图和去平干涉相位图进行滤波，根据第一波段、第二波段的滤波后的主图像的幅度图，分别对第一波段、第二波段的滤波后的去平干涉图进行配准，由于每个波段的滤波后的主图像的幅度图和滤波后的去平干涉相位图的信噪比提高，相应地，第一波段、第二波段的滤波后的去平干涉相位图的配准精度也被提高，最后对第一个波段的配准后的去平干涉相位图和第二波段的配准后的去平干涉相位图进行相位解缠，对相位解缠结果进行高程反演，能够获得更加准确的目标区域的高程数据，实现了高精度的高程反演。
实施例二
为了能够更加体现本发明的目的，在上述实施例的基础上，进行进一步的举例说明。
本发明实施例提供一种双频干涉合成孔径雷达的高程反演方法，如图4所示，所述方法包括：
S401-S403：分别与S101至S103的实现方式相同，这里不再赘述。
S404：采用第一波段的配准后的两幅SAR图像之间的相干系数和第二波段的配准后的两幅SAR图像之间的相干系数，对第一波段的配准后的去平干涉相位图和第二波段的配准后的去平干涉相位图进行相位解缠，得到相位解缠结果，其中，相位解缠结果包括：每个像素点的解缠后的真实相位。
S405：对每个像素点的解缠后的真实相位进行加权均值滤波，得到每个像素点的滤波后的真实相位。
在实际实施时，相位解缠结果还包括：每个像素点的解缠后的真实相位对应的相位周期数；
相应地，根据所有像素点的解缠后的真实相位和解缠后的真实相位对应的相位周期数、以及预设的有效像素点的确定方式，确定出所有像素点中的有效像素点；根据所有像素点中的有效像素点，对每个像素点的解缠后的真实相位进行加权均值滤波，得到每个像素点的滤波后的真实相位。
S406：对每个像素点的滤波后的真实相位进行高程反演，得到数字高程模型。
需要说明的是，步骤S401至步骤S406均可以由处理器实现。
由此可见，本发明实施例中，通过对每个波段的两幅SAR图像中的主图像的幅度图和去平干涉相位图进行滤波，根据第一波段、第二波段的滤波后的主图像的幅度图，分别对第一波段、第二波段的滤波后的去平干涉相位图进行配准，由于每个波段的滤波后的主图像的幅度图和滤波后的去平干涉相位图的信噪比提高，相应地，第一波段、第二波段的滤波后的去平干涉相位图的配准精度也被提高，再对第一波段的配准后的去平干涉相位图和第二波段的配准后的去平干涉相位图进行相位解缠，对相位解缠结果进行加权均值滤波，进一步减小噪点的干扰，提高方法的鲁棒性，最后对每个像素点的滤波后的真实相位进行高程反演，能够获得更加准确的目标区域的高程数据，从而实现了高精度、强鲁棒性的高程反演。
实施例三
为了能够更加体现本发明的目的，在上述实施例的基础上，进行进一步的举例说明。
本发明实施例提供一种双频干涉合成孔径雷达的高程重建方法，该方法包括：
步骤A1：分别获取第一波段的两幅SAR图像和第二波段的两幅SAR图像，对每个波段的两幅SAR图像依次进行配准和干涉处理，得到每个波段的配准后的两幅SAR图像和每个波段的配准后的干涉图。
在实际实施时，第一波段可以是C波段，第二波段可以是X波段；通常，两幅SAR图像的配准误差必须在1/8个像素元以下，才能保证由每个波段的配准后的两幅SAR图像生成的配准后的干涉图中的干涉条纹清晰度较佳，可选地，频域SVD(Singular ValueDecomposition，奇异值分解)在配准处理时具有较高的效率和精度，能够满足图像的配准误差在1/8个像素元以下，因此，本发明实施例中可以采用频域SVD对每个波段的两幅SAR图像进行配准，得到每个波段的配准后的两幅SAR图像；再对每个波段的配准后的两幅SAR图像进行干涉处理，具体包括：将配准后的两幅SAR图像进行相减，或者，对配准后的两幅SAR图像进行共轭相乘，得到每个波段的配准后的干涉图。
进一步地，还对每个波段的配准后的两幅SAR图像进行相干系数的计算，得到每个波段的两幅SAR图像之间的相干系数图。
示例性地，所述机载DFInSAR系统通过双天线获取的C波段的两幅SAR图像，如图2a、图2b所示，分别对C波段的两幅SAR图像中的主图像和辅图像进行配准，得到主图像和被配准到主图像的辅图像；对配准后的主图像和配准后的辅图像进行干涉处理，得到C波段的配准后的干涉图如图5所示，图5中数值表示相位，。
步骤A2：对第一波段的配准后的干涉图和第二波段的配准后的干涉图进行平地效应的去除，得到每个波段的去平干涉相位图。
在实际实施时，去除平地效应的方法主要有频移法和基于轨道参数直接计算法，其中，频移法是将干涉图变换到频域，假设平地相位在频域是唯一占优频率，将其搬移到零频率来实现平地效应的去除；
本发明实施例中的机载DFInSAR系统包括两个雷达天线A1、A2，两个雷达天线在每个波段下向目标区域发射微波，然后接收该目标区域反射的回波，获得该目标区域的每个波段的两幅复图像，示例性地，以目标区域中的任一目标点p为例，两个雷达天线和目标点p存在一定的成像几何关系，如图6所示，图6中，B表示两个雷达天线之间的基线距离，Bv表示垂直基线，Bh是水平基线，α表示基线与参考平面之间的基线角，θ表示雷达天线A1相位中心到目标的下视角，r1、r2分别表示雷达天线A1、A2到目标点的倾距，Δr表示两个斜距之间的斜距差，H表示雷达平台相对于大地参考平面的高度，H1雷达平台相对参考平面的高度，参考平面是根据谷歌地球来获取的，hflat表示参考平面相对于大地参考平面的平地高度，h表示待求的目标点p的海拔高度；基于轨道参数直接计算法是根据该成像几何关系，计算得到平地相位
Figure BDA0001604266710000131
再对干涉图减去平地相位
Figure BDA0001604266710000132
得到去平干涉相位图，其中，平地相位
Figure BDA0001604266710000133
表示为：
Figure BDA0001604266710000134
λ表示雷达波长，c表示雷达天线工作模式，当雷达天线工作在标准模式，即单发双收，c取值为1，当雷达天线工作在乒乓模式，即自发自收，c取值为2。
示例性地，分别采用频移法和基于轨道参数直接计算法，对C波段的配准后的干涉图进行平地效应的去除，其中，频移法得到的去平干涉相位图如图7所示，基于轨道参数直接计算法得到的去平干涉相位图如图8所示，图7和图8中的数值表示相位；从图7和图8中可以看出，由于实际中平地相位并不是一个单一的条纹频率，相较于基于轨道参数直接计算法得到的去平干涉相位图，频移法得到的去平干涉相位图中的条纹更细、更加密集，没有有效地去除平地效应。
步骤A3：对每个波段的两幅SAR图像中的主图像的幅度图、去平干涉相位图和两幅SAR图像之间的相干系数图，分别进行滤波，得到每个波段的滤波后的主图像的幅度图、滤波后的去平干涉相位图和滤波后的两幅SAR图像之间的相干系数图。
在实际实施时，在对X波段和C波段的去平干涉相位图进行配准之前，对X波段和C波段的主图像的幅度图进行滤波，优选地，采用NL-InSAR(Non-local filtering InSAR，非局部滤波的干涉合成孔径雷达)方法，对每个波段的去平干涉相位图进行观测相位的参数估计和滤波，得到每个波段的滤波后的去平干涉相位图；对每个波段的两幅SAR图像中的主图像的幅度图进行幅度的参数估计和滤波，得到每个波段的滤波后的幅度图；对每个波段的两幅SAR图像之间的相干系数图，进行相干系数的参数估计和滤波，得到每个波段的滤波后的两幅SAR图像之间的相干系数图。
示例性地，采用NL-InSAR方法，对图9a所示的C波段的两幅SAR图像中的主图像的幅度图进行幅度的参数估计和滤波，得到C波段的滤波后的主图像的幅度图如图9b所示，对比表明NL-InSAR方法在保持SAR图像有效信息的前提下，实现了对SAR图像中主图像的幅度图的滤波；对图10a所示的基于轨道参数直接计算法得到的C波段的去平干涉相位图，进行观测相位的参数估计和滤波，得到C波段的滤波后的去平干涉相位图如图10b所示，图10a、图10b中的数值表示相位，可以看出NL-InSAR方法在去除观测相位噪声的同时较好地保持了去平干涉相位图的纹理信息；对图11a所示的C波段的两幅SAR图像之间的相干系数图，进行相干系数的参数估计和滤波，得到C波段的滤波后的两幅SAR图像之间的相干系数图如图11b所示，图11a、图11b中的数值表示相干系数，可以看出C波段的滤波后的两幅SAR图像之间的相干系数的信息保留的更加全面、可靠性更强，可以有效辅助后续的相位解缠。
步骤A4：根据所述第一波段的滤波后的主图像的幅度图和所述第二波段的滤波后的主图像的幅度图，分别对所述第一波段的滤波后的去平干涉相位图和所述第二波段的滤波后的去平干涉相位图进行配准，得到第一波段的配准后的去平干涉相位图和第二波段的配准后的去平干涉相位图。
在实际实施时，由于相位解缠需要X波段的去平干涉相位图和C波段的去平干涉相位图配准精度达到亚像素级，因此在进行相位解缠之前需要分别对C波段、X波段的去平干涉相位图进行配准；并分别对C波段的滤波后的两幅SAR图像之间的相干系数图和X波段的滤波后的两幅SAR图像之间的相干系数图进行配准，得到C波段的配准后的两幅SAR图像之间的相干系数图和X波段的配准后的两幅SAR图像之间的相干系数图。
示例性地，采用的配准方法是对C波段的滤波后的主图像的幅度图和X波段的滤波后的主图像的幅度图拟合出多项式系数，根据该多项式系数分别对C波段、X波段滤波后的去平干涉相位图进行配准，得到C波段的配准后的去平干涉相位图和X波段的配准后的去平干涉相位图，分别如图12a、图12b所示，图12a、图12b中的数值表示相位。
对每个波段的滤波后的主图像的幅度图分别叠加自身的配准前和配准后的去平干涉相位图，得到每个波段的配准前和配准后的复干涉图，再分别计算C波段的配准前的复干涉图和X波段的配准前的复干涉图之间的相干系数、C波段的配准后的复干涉图和X波段的配准后的复干涉图之间的相干系数，得到C波段和X波段的配准前后的复干涉图之间的相干系数图，分别如图13a、图13b所示，图13a、图13b中的数值表示相干系数，可以看出，C波段和X波段的配准前的复干涉图之间的相干系数较低，均值大约为0.3，C波段和X波段的配准后的复干涉图之间的相干系数较高，均值大约为0.86，C波段和X波段的配准后的复干涉图之间的相干系数被提高，基于C波段、X波段的配准后的去平干涉相位图进行联合解缠后，得到的相位解缠结果也会更加准确。
步骤A5：采用C波段的配准后的两幅SAR图像之间的相干系数和X波段的配准后的两幅SAR图像之间的相干系数，对第一波段的配准后的去平干涉相位图和第二波段的配准后的去平干涉相位图进行相位解缠，得到相位解缠结果，其中，相位解缠结果包括：每个像素点的解缠后的真实相位。
在实际实施时，相位解缠算法主要分为统计信号处理法和非统计信号处理法，非统计信号处理主要是基于中国剩余定理CRT，统计信号处理的方法主要包括了最大似然估计MLE和整体变分最大后验估计TV-MAP两类，通过对一组无噪声无缠绕的初始相位数据进行相位缠绕，得到一组C波段的缠绕后的相位数据和一组X波段的缠绕后的相位数据，再分别采用上述三种方法对这两组缠绕后的相位数据进行相位解缠，得到CRT、MLE、TV-MAP的相位解缠结果，分别如图14a、图14b、图14c所示，图14a、图14b、图14c中的数值表示相位，最后分别对CRT、MLE、TV-MAP的相位解缠结果减去所述初始相位数据，得到对CRT、MLE、TV-MAP的残差图，分别如图15a、图15b、图15c所示，图15a、图15b、图15c中的数值表示相位，可以看出，上述三种方法中的整体变分最大后验估计的残差图中均值大约为0，即采用整体变分最大后验估计得到的相位解缠结果最接近初始相位数据，具有较好的鲁棒性。
可选地，本发明实施例中采用整体变分最大后验估计，对所述第一波段的配准后的去平干涉相位图和所述第二波段的配准后的去平干涉相位图进行相位解缠，得到相位解缠结果，其中，所述相位解缠结果包括：每个像素点的解缠后的真实相位，每个像素点的解缠后的真实相位对应的相位周期数。
具体地，整体变分最大后验估计使用的整体变分(Total Variation，TV)模型，其能量函数Epriori定义为：
Figure BDA0001604266710000161
公式(2)中β是一个固定的常量，用来调节真实相位的光滑性及观测相位与真实相位之间相似性之间的平衡的常数，Np为像素点p的邻域的集合，像素点q为集合Np中的任一个像素点，
Figure BDA0001604266710000162
和
Figure BDA0001604266710000163
分别表示像素点p、q的先验相位值，N×M表示图像大小，wp,q表征邻域连通性，当Np为4邻域模型，wp,q值等于1；当Np为8邻域模型时，wp,q值等于
Figure BDA0001604266710000164
上述TV模型是一种非局部模型，其特点是运算速度快，无需对超参数β进行多次估计，并能够有效应对断层和平滑的区域；
利用最大后验估计可以得到整体变分最大后验估计在TV模型下的解缠后的真实相位
Figure BDA0001604266710000165
解缠后的真实相位
Figure BDA0001604266710000166
表示为：
Figure BDA0001604266710000167
公式(3)中，Φ＝Ψ+2Kπ，进一步地，解缠后的真实相位
Figure BDA0001604266710000168
表示为：
Figure BDA0001604266710000169
公式(3)和(4)中，K表示每个像素点的解缠后的真实相位对应的相位周期数，I表示波段数，本发明实施例中I等于2，当i＝1时，表示C波段或X波段，对应地，当i＝2时，表示X波段或C波段；Ψ表示C波段或X波段的配准后的去平干涉相位图中的观测相位，Ψpi(p)表示i对应的波段的配准后的去平干涉相位图中像素点p的观测相位，f(Ψpi(p)/Φpi)表示采用i对应的波段的配准后的两幅SAR图像之间的相干系数的概率密度函数；
可以看出，上述优化模型将求解每个像素点的解缠后的真实相位
Figure BDA0001604266710000171
的问题，转化为采用图割法求解使能量最小化的相位周期数K的问题，具体地，对每个像素点改变相位周期数K值，计算得到一个使公式(4)右侧的值为最小值的相位周期数K，即为该像素点的解缠后的真实相位对应的相位周期数K，进而得到该像素点的解缠后的真实相位
Figure BDA0001604266710000172
进一步地，得到的每个像素点的真实相位对应的相位周期数K可能是该像素点的C波段的相位周期数KC，可以根据该像素点的C波段的相位周期数KC，获得该像素点的X波段的相位周期数KX；同样地，得到的每个像素点的真实相位对应的相位周期数K可能是该像素点的X波段的相位周期数KC，可以根据该像素点的X波段的相位周期数KX，获得该像素点的C波段的相位周期数KC。
步骤A6：对所述每个像素点的解缠后的真实相位进行加权均值滤波，得到每个像素点的滤波后的真实相位。
在实际实施时，对所述每个像素点的解缠后的真实相位进行加权均值滤波之前，根据每个像素点的解缠后的真实相位、C波段的相位周期数KC和X波段的相位周期数KX，对每个像素点进行有效点的判断，确定出所有像素点中的有效像素点，具体包括以下几步：
预设模糊矢量V表示为：
V＝[KC,KX] (5)
计算每个像素点的总相位周期数n表示为：
Figure BDA0001604266710000173
公式(5)和(6)中，Vi为V的第i个分量，I表示波段数，本发明实施例中I＝2，V1就是V的第1个分量KC，V2就是V的第2个分量KX，所述每个像素点的C波段的相位周期数KC和X波段的相位周期数KX相加，得到每个像素点的总相位周期数n；
根据所述每个像素点的总相位周期数，对所有像素点进行聚类分组，得到包含不同像素点的分组，当每个分组内的像素点点数小于预设第一阈值th1时，确定该分组内的所有像素点为噪点，当每个分组内的像素点点数不小于预设第一阈值th1时，确定该分组内的所有像素点为第一有效点；
每个第一有效点的解缠后的真实相位减去自身邻域内的所有第一有效点的解缠后的真实相位的均值，得到每个第一有效点的相位差值，每个第一有效点的邻域内的所有第一有效点的解缠后的真实相位的均值乘以预设第二阈值th2，得到每个第一有效点的均值倍数，当每个第一有效点的相位差值大于自身的均值倍数时，确定该第一有效点为噪点，当每个第一有效点的相位差值不大于自身的均值倍数时，确定该第一有效点为第二有效点；
每个第二有效点的总相位周期数，减去自身邻域内第二有效点的解缠后的真实相位的均值，得到每个第二有效点的相位差值，当每个第二有效点的相位差值大于预设第三阈值th3时，确定该第二有效点为噪点，当每个第二有效点的相位差值不大于预设第三阈值th3时，确定该第二有效点为第三有效点，将所有第三有效点作为有效像素点。
进一步地，完成有效点判断之后，根据所有有效像素点，对所述每个像素点的解缠后的真实相位进行加权均值滤波，得到所述每个像素点的滤波后的真实相位Φ’，任一个像素点p的滤波后的真实相位Φp’表示为：
Figure BDA0001604266710000181
公式(7)中，
Figure BDA0001604266710000182
为像素点q的解缠后的真实相位，q为像素点p的邻域Np中的任一个像素点，当q为噪点时，mq＝0；当q为有效像素点时，mq＝1。
示例性地，所有像素点的解缠后的真实相位图16a所示，采用上述方法对每个像素点的解缠后的真实相位进行加权均值滤波，得到所有像素点的滤波后的真实相位图如图16b所示，图16a和图16b中的数值表示相位、数值单位为rad，可以看出，对于未滤波之前所有像素点的解缠后的真实相位图中的孤立的黑色噪点，该方法具有较好去噪效果。
步骤A7：对所述每个像素点的滤波后的真实相位进行高程反演，得到所述数字高程模型。
在实际实施时，将滤波后的真实相位Φ’加上平地相位
Figure BDA0001604266710000191
得到带平地的真实相位
Figure BDA0001604266710000192
再根据图6所示的几何关系，对带平地的真实相位
Figure BDA0001604266710000193
进行高程反演，得到目标区域的高程h，目标区域的高程h表示为：
Figure BDA0001604266710000194
示例性地，获得某山区场景如图17所示，选取其中地形起伏剧烈的包含高架路段的部分区域，采用光学摄影系统获取所述部分区域的光学图像的灰度图，如图18所示，再分别采用基于MCF(Minimum cost flow，最小费用流)的单波段InSAR方法和本发明实施例中的双波段DFInSAR方法，对所述部分区域进行相位解缠，获取对应的相位解缠结果，单波段InSAR方法获取的相位解缠结果如图19a所示、本发明实施例的双波段DFInSAR方法获取的相位解缠结果如图19b所示；可以看出，参照该高架区域的光学图像，所述单波段InSAR方法获取的相位解缠结果中，所述部分区域中方框内的高架路段与方框外的高架路段的颜色深度不一样，即高程信息一样的高架路段得到了不同的相位解缠结果，采用所述单波段InSAR方法获取的高架路段的相位解缠结果不正确，而本发明实施例的双波段DFInSAR方法获取的相位解缠结果中，所述部分区域中方框内的高架路段与方框外的高架路段的颜色深度基本一样，即高程信息一样的高架路段得到了相同的相位解缠结果，实现了对该高架路段的正确相位解缠。
示例性地，在某一个已知数字高程模型的区域中，选取地形走势沿虚线指示的目标区域，如图20所示，图20中数值表示高程，单位为m，分别采用C波段的InSAR方法、X波段的InSAR方法以及本发明实施例的C波段和X波段的DFInSAR方法，对该目标区域进行高程反演，获得沿虚线指示的目标区域的高程h，结果如图21所示，纵坐标表示沿虚线指示的目标区域的高程，横坐标表示以沿虚线指示的目标区域的最左端为坐标原点、沿虚线指示的目标区域在水平方向上的路径长度；可以看出，图20中沿虚线指示的目标区域的颜色从浅色到深色、再到浅色，再结合图20中高程h与颜色的对应关系，可以获知，沿虚线指示的目标区域的实际高程信息的变化过程大致为：从650m到770m、再到650m，而图21中显示的本发明实施例的C波段和X波段的DFInSAR方法获取沿虚线指示的目标区域的高程h数据，更加符合沿虚线指示的目标区域的实际高程信息，实现了精确的高程反演。
由此可见，本发明实施例中，通过对每个波段的两幅SAR图像中的主图像的幅度图和去平干涉相位图进行滤波，根据第一波段、第二波段的滤波后的主图像的幅度图，对第一波段、第二波段的滤波后的去平干涉相位图进行配准，再对第一波段的配准后的去平干涉相位图和第二波段的配准后的去平干涉相位图进行相位解缠，对相位解缠结果进行加权均值滤波，最后对每个像素点的滤波后的真实相位进行高程反演，通过实验仿真分析，可以看到，本发明实施例DFInSAR的高程重建方法能够获得更加准确的目标区域的高程数据，实现了高精度、强鲁棒性的高程反演。
实施例四
为了能够更加体现本发明的目的，在前述方法实施例的基础上，进行进一步的举例说明。
本发明实施例提供一种双频干涉合成孔径雷达的高程重建装置，如图22所示，该高程重建装置2200包括：存储器2201和处理器2202；其中，
存储器2201用于存储计算机程序；
处理器2202用于在运行存储器2201中存储的计算机程序时，执行以下步骤：
分别获取第一波段的两幅合成孔径雷达SAR图像和第二波段的两幅SAR图像；
对每个波段的两幅SAR图像进行处理，得到每个波段的去平干涉相位图；
对每个波段的两幅SAR图像中的主图像的幅度图和去平干涉相位图分别进行滤波，得到每个波段的滤波后的主图像的幅度图和滤波后的去平干涉相位图；
根据第一波段的滤波后的主图像的幅度图和第二波段的滤波后的主图像的幅度图，分别对每个波段的滤波后的去平干涉相位图进行配准，得到每个波段的配准后的去平干涉相位图；
对第一波段的配准后的去平干涉相位图和第二波段的配准后的去平干涉相位图进行相位解缠，得到相位解缠结果；
对相位解缠结果进行高程反演，得到数字高程模型。
在上述方案中，相位解缠结果包括：每个像素点的解缠后的真实相位；相应地，处理器2202具体用于在运行存储器2201中存储的计算机程序时，执行以下步骤：
对每个像素点的解缠后的真实相位进行加权均值滤波，得到每个像素点的滤波后的真实相位；
对每个像素点的滤波后的真实相位进行高程反演，得到数字高程模型。
在上述方案中，相位解缠结果还包括：每个像素点的解缠后的真实相位对应的相位周期数；
相应地，处理器2202具体用于在运行存储器2201中存储的计算机程序时，执行以下步骤：根据所有像素点的解缠后的真实相位和解缠后的真实相位对应的相位周期数、以及预设的有效像素点的确定方式，确定出所有像素点中有效的像素点；根据所有像素点中有效的像素点，对每个像素点的解缠后的真实相位进行加权均值滤波，得到每个像素点的滤波后的真实相位。
示例性地，处理器2202具体用于在运行存储器2201中存储的计算机程序时，执行以下步骤：采用非局部滤波的干涉合成孔径雷达NL-InSAR方法，对每个波段的两幅SAR图像中的主图像的幅度图和去平干涉相位图分别进行滤波。
示例性地，处理器2202具体用于在运行存储器2201中存储的计算机程序时，执行以下步骤：对所述每个波段的两幅SAR图像依次进行配准，得到每个波段的配准后的两幅SAR图像；对所述每个波段的配准后的两幅SAR图像进行干涉处理，获得每个波段的配准后的干涉图；对所述每个波段的配准后的干涉图进行平地效应的去除，得到每个波段的去平干涉相位图。
实施例五
基于与前述实施例相同的技术构思，本发明实施例五提供了一种计算机可读存储介质，可以应用于装置中；前述实施例的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的全部或部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个计算机可读存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)或processor(处理器)执行本实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read Only Memory)、随机存取存储器(RAM，Random Access Memory)、磁碟或者光盘等各种可以存储程序代码的介质。
具体来讲，本实施例中的一种双频干涉合成孔径雷达的高程重建方法对应的计算机程序指令，可以被存储在光盘，硬盘，U盘等存储介质上，当存储介质中的与一种双频干涉合成孔径雷达的高程重建方法对应的计算机程序指令被一电子设备读取或被执行时，导致所述至少一个处理器执行本发明前述实施例的任意一种双频干涉合成孔径雷达的高程重建方法所述的步骤。
需要说明的是，在本文中，术语“包括”、“包含”或者其任何其他变体意在涵盖非排他性的包含，从而使得包括一系列要素的过程、方法、物品或者装置不仅包括那些要素，而且还包括没有明确列出的其他要素，或者是还包括为这种过程、方法、物品或者装置所固有的要素。在没有更多限制的情况下，由语句“包括一个……”限定的要素，并不排除在包括该要素的过程、方法、物品或者装置中还存在另外的相同要素。
上述本发明实施例序号仅仅为了描述，不代表实施例的优劣。
通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到上述实施例方法可借助软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件，但很多情况下前者是更佳的实施方式。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质(如ROM/RAM、磁碟、光盘)中，包括若干指令用以使得一台终端(可以是手机，计算机，服务器，空调器，或者网络设备等)执行本发明各个实施例所述的方法。
上面结合附图对本发明的实施例进行了描述，但是本发明并不局限于上述的具体实施方式，上述的具体实施方式仅仅是示意性的，而不是限制性的，本领域的普通技术人员在本发明的启示下，在不脱离本发明宗旨和权利要求所保护的范围情况下，还可做出很多形式，这些均属于本发明的保护之内。

Claims (8)

1.一种双频干涉合成孔径雷达的高程重建方法，其特征在于，所述方法包括：
分别获取第一波段的两幅合成孔径雷达SAR图像和第二波段的两幅SAR图像，对每个波段的两幅SAR图像进行处理，得到每个波段的去平干涉相位图；
对每个波段的两幅SAR图像中的主图像的幅度图和去平干涉相位图分别进行滤波，得到每个波段的滤波后的主图像的幅度图和滤波后的去平干涉相位图；
其中，所述对每个波段的两幅SAR图像中的主图像的幅度图和去平干涉相位图分别进行滤波，包括：
采用非局部滤波的干涉合成孔径雷达NL-InSAR方法，对每个波段的两幅SAR图像中的主图像的幅度图和去平干涉相位图分别进行滤波；
根据第一波段的滤波后的主图像的幅度图和第二波段的滤波后的主图像的幅度图，分别对所述每个波段的滤波后的去平干涉相位图进行配准，得到每个波段的配准后的去平干涉相位图；对第一波段的配准后的去平干涉相位图和第二波段的配准后的去平干涉相位图进行相位解缠，得到相位解缠结果；
对所述相位解缠结果进行高程反演，得到数字高程模型。
2.根据权利要求1所述的方法，其特征在于，所述相位解缠结果包括：每个像素点的解缠后的真实相位；
相应地，所述对所述相位解缠结果进行高程反演，得到数字高程模型，包括：
对所述每个像素点的解缠后的真实相位进行加权均值滤波，得到每个像素点的滤波后的真实相位；
对所述每个像素点的滤波后的真实相位进行高程反演，得到所述数字高程模型。
3.根据权利要求2所述的方法，其特征在于，所述相位解缠结果还包括：
每个像素点的解缠后的真实相位对应的相位周期数；
相应地，所述对所述每个像素点的滤波后的真实相位进行加权均值滤波，得到每个像素点的滤波后的真实相位，包括：
根据所有像素点的解缠后的真实相位和解缠后的真实相位对应的相位周期数、以及预设的有效像素点的确定方式，确定出所有像素点中有效的像素点；
根据所有像素点中有效的像素点，对所述每个像素点的解缠后的真实相位进行加权均值滤波，得到所述每个像素点的滤波后的真实相位。
4.根据权利要求1所述的方法，其特征在于，所述对每个波段的两幅SAR图像进行处理，得到每个波段的去平干涉相位图，包括：
对所述每个波段的两幅SAR图像依次进行配准，得到每个波段的配准后的两幅SAR图像；
对所述每个波段的配准后的两幅SAR图像进行干涉处理，获得每个波段的配准后的干涉图；
对所述每个波段的配准后的干涉图进行平地效应的去除，得到每个波段的去平干涉相位图。
5.一种双频干涉合成孔径雷达的高程重建装置，其特征在于，所述装置包括：存储器和处理器；其中，
所述存储器，用于存储计算机程序；
所述处理器，用于在运行所述计算机程序时，执行以下步骤：
分别获取第一波段的两幅合成孔径雷达SAR图像和第二波段的两幅SAR图像；
对每个波段的两幅SAR图像进行处理，得到每个波段的去平干涉相位图；
对每个波段的两幅SAR图像中的主图像的幅度图和去平干涉相位图分别进行滤波，得到每个波段的滤波后的主图像的幅度图和滤波后的去平干涉相位图；
其中，所述对每个波段的两幅SAR图像中的主图像的幅度图和去平干涉相位图分别进行滤波，包括：
采用非局部滤波的干涉合成孔径雷达NL-InSAR方法，对每个波段的两幅SAR图像中的主图像的幅度图和去平干涉相位图分别进行滤波；
根据第一波段的滤波后的主图像的幅度图和第二波段的滤波后的主图像的幅度图，分别对所述每个波段的滤波后的去平干涉相位图进行配准，得到每个波段的配准后的去平干涉相位图；
对第一波段的配准后的去平干涉相位图和第二波段的配准后的去平干涉相位图进行相位解缠，得到相位解缠结果；
对所述相位解缠结果进行高程反演，得到数字高程模型。
6.根据权利要求5所述的装置，其特征在于，所述相位解缠结果包括：每个像素点的解缠后的真实相位；
相应地，所述处理器具体用于在运行所述计算机程序时，执行以下步骤：
对所述每个像素点的解缠后的真实相位进行加权均值滤波，得到每个像素点的滤波后的真实相位；
对所述每个像素点的滤波后的真实相位进行高程反演，得到所述数字高程模型。
7.根据权利要求6所述的装置，其特征在于，所述相位解缠结果还包括：每个像素点的解缠后的真实相位对应的相位周期数；
相应地，所述处理器具体用于在运行所述计算机程序时，执行以下步骤：
根据所有像素点的解缠后的真实相位和解缠后的真实相位对应的相位周期数、以及预设的有效像素点的确定方式，确定出所有像素点中有效的像素点；
根据所有像素点中有效的像素点，对所述每个像素点的解缠后的真实相位进行加权均值滤波，得到所述每个像素点的滤波后的真实相位。
8.一种计算机可读存储介质，其特征在于，所述计算机可读存储介质存储有计算机程序，
当所述计算机程序被至少一个处理器执行时，导致所述至少一个处理器执行权利要求1至4任一项所述的方法的步骤。
CN201810237097.6A 2018-03-21 2018-03-21 一种双频干涉合成孔径雷达的高程重建方法和装置 Active CN108594222B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810237097.6A CN108594222B (zh) 2018-03-21 2018-03-21 一种双频干涉合成孔径雷达的高程重建方法和装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810237097.6A CN108594222B (zh) 2018-03-21 2018-03-21 一种双频干涉合成孔径雷达的高程重建方法和装置

Publications (2)

Publication Number Publication Date
CN108594222A CN108594222A (zh) 2018-09-28
CN108594222B true CN108594222B (zh) 2020-05-08

Family

ID=63627226

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810237097.6A Active CN108594222B (zh) 2018-03-21 2018-03-21 一种双频干涉合成孔径雷达的高程重建方法和装置

Country Status (1)

Country Link
CN (1) CN108594222B (zh)

Families Citing this family (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP6584369B2 (ja) * 2015-08-20 2019-10-02 三菱電機株式会社 壁の背後に位置するシーンの雑音のない画像を求めるシステムおよび方法
CN109541596B (zh) * 2018-11-28 2022-05-20 中国电子科技集团公司电子科学研究院 基于深度学习算法的InSAR图像处理方法及装置
CN111208512B (zh) * 2020-01-15 2022-06-07 电子科技大学 一种基于视频合成孔径雷达的干涉测量方法
CN113589282B (zh) * 2021-07-12 2023-08-08 中国科学院国家空间科学中心 一种基于图像域变换的星载干涉成像高度计去平地效应方法
CN113835089A (zh) * 2021-08-27 2021-12-24 中国空间技术研究院 一种基于差频InSAR的缓变高程反演方法
CN114089336B (zh) * 2021-11-19 2024-06-18 北京环境特性研究所 一种复杂起伏场景的自适应分区方法及雷达回波仿真方法
CN114609633B (zh) * 2022-03-17 2023-09-01 电子科技大学 一种圆周聚束模式干涉sar测高方法
CN114859346B (zh) * 2022-04-21 2024-05-03 桂林电子科技大学 一种基于Insarbm3d滤波算法的迭代反演滤波方法
CN115291217B (zh) * 2022-10-10 2022-12-16 中国科学院空天信息创新研究院 一种双频极化干涉四阶段树高反演方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1536247A1 (de) * 2003-11-28 2005-06-01 Deutsches Zentrum für Luft- und Raumfahrt e.V. Interferometrisches Mikrowellen-Radarverfahren
CN101866002A (zh) * 2010-06-01 2010-10-20 中国人民解放军信息工程大学 基于中国余数定理的多基线、多波段InSAR相位解缠方法
WO2014039267A1 (en) * 2012-09-04 2014-03-13 Fugro Earthdata, Inc. Method and apparatus for mapping and characterizing sea ice from airborne simultaneous dual frequency interferometric synthetic aperture radar (ifsar) measurements
CN106707281A (zh) * 2017-01-05 2017-05-24 北京航空航天大学 一种基于多频数据处理的机载D‑InSAR形变检测方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
EP1536247A1 (de) * 2003-11-28 2005-06-01 Deutsches Zentrum für Luft- und Raumfahrt e.V. Interferometrisches Mikrowellen-Radarverfahren
CN101866002A (zh) * 2010-06-01 2010-10-20 中国人民解放军信息工程大学 基于中国余数定理的多基线、多波段InSAR相位解缠方法
WO2014039267A1 (en) * 2012-09-04 2014-03-13 Fugro Earthdata, Inc. Method and apparatus for mapping and characterizing sea ice from airborne simultaneous dual frequency interferometric synthetic aperture radar (ifsar) measurements
CN106707281A (zh) * 2017-01-05 2017-05-24 北京航空航天大学 一种基于多频数据处理的机载D‑InSAR形变检测方法

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Combination of repeat and single-pass dual-frequency airborne InSAR data for accurate height estimation;M. Pinheiro等;《2014 IEEE Geoscience and Remote Sensing Symposium》;20140718;2237-2240 *
多波段InSAR差分滤波相位解缠方法;靳国旺等;《测绘学报》;20120630;第41卷(第3期);434-440 *

Also Published As

Publication number Publication date
CN108594222A (zh) 2018-09-28

Similar Documents

Publication Publication Date Title
CN108594222B (zh) 一种双频干涉合成孔径雷达的高程重建方法和装置
EP3229038B1 (en) Wavelet domain insar interferometric phase filtering method in combination with local frequency estimation
CN108564532B (zh) 大尺度地距星载sar图像镶嵌方法
CN108171647B (zh) 一种考虑地表形变的Landsat 7条带影像修复方法
CN108663678B (zh) 基于混合整数优化模型的多基线InSAR相位解缠算法
CN111856459B (zh) 一种改进的DEM最大似然约束多基线InSAR相位解缠方法
CN111768337B (zh) 一种图像处理方法、装置和电子设备
Yamaki et al. Singular unit restoration in interferograms based on complex-valued Markov random field model for phase unwrapping
CN113311433B (zh) 一种质量图和最小费用流结合的InSAR干涉相位两步解缠方法
CN109239710B (zh) 雷达高程信息的获取方法及装置、计算机可读存储介质
Fang et al. An InSAR fine registration algorithm using uniform tie points based on Voronoi diagram
CN115951350A (zh) 永久散射体点提取方法、装置、设备及介质
CN109886910B (zh) 外部数字高程模型dem修正方法及装置
CN115170406A (zh) 一种光学影像与sar强度影像高精度融合方法
CN110703252A (zh) 干涉合成孔径雷达阴影区域数字高程模型修正方法
CN108876829B (zh) 基于非线性尺度空间及径向基函数的sar高精度配准方法
CN115616575B (zh) 一种星载sar立体测量辅助的干涉相位图解缠绕方法
CN112946647A (zh) 大气误差改正InSAR干涉图堆叠地质灾害普查方法和装置
Hu et al. Planetary3D: A photogrammetric tool for 3D topographic mapping of planetary bodies
Li et al. Co-registration and residual correction of digital elevation models: A comparative study
CN115546264A (zh) 一种星载InSAR图像精细配准与立体测量方法
CN113835089A (zh) 一种基于差频InSAR的缓变高程反演方法
CN109839635B (zh) 一种通过CryoSat-2 SARIn模式L1b级波形数据提取测高脚点高程的方法
Xie Enhanced multi-baseline unscented Kalman filtering phase unwrapping algorithm
CN114609635B (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