CN109212522B - 一种基于双基星载sar的高精度dem反演方法及装置 - Google Patents

一种基于双基星载sar的高精度dem反演方法及装置 Download PDF

Info

Publication number
CN109212522B
CN109212522B CN201810525764.0A CN201810525764A CN109212522B CN 109212522 B CN109212522 B CN 109212522B CN 201810525764 A CN201810525764 A CN 201810525764A CN 109212522 B CN109212522 B CN 109212522B
Authority
CN
China
Prior art keywords
dem
image
baseline
interference
pixel
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
CN201810525764.0A
Other languages
English (en)
Other versions
CN109212522A (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 CN201810525764.0A priority Critical patent/CN109212522B/zh
Publication of CN109212522A publication Critical patent/CN109212522A/zh
Application granted granted Critical
Publication of CN109212522B publication Critical patent/CN109212522B/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/87—Combinations of radar systems, e.g. primary radar and secondary radar
    • 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

Abstract

本发明公开了一种基于双基星载SAR的高精度DEM反演方法、装置、存储介质和信息处理装置，获取第一星载合成孔径雷达(SAR)和第二星载SAR在第一轨道对成像目标分别进行成像得到的第一主图像和第一辅图像，并获取第一星载SAR和第二星载SAR在第二轨道对所述成像目标分别进行成像得到的第二主图像和第二辅图像；生成所述第一轨道对应的第一干涉图，和第二轨道对应的第二干涉图；再进行滤波处理，并确定所述两个干涉图上各像素对应的相干系数；将所述滤波后的第一干涉图和第二干涉图，分别转化为第一单基线数字高程模型(DEM)将和第二单基线DEM；根据各像素对应的相干系数，融合所述第一单基线DEM和第二单基线DEM，得到成像目标的融合DEM。

Description

一种基于双基星载SAR的高精度DEM反演方法及装置
技术领域
本发明涉及合成孔径雷达(SAR，Synthetic Aperture Radar)技术，尤其涉及一种基于双基星载SAR的高精度DEM反演方法及装置。
背景技术
数字高程模型(DEM，Digital Elevation Model)是描述地球表面形状的三维数字模型，由一系列包含有地理平面坐标和高程的数据集组成；DEM在科学研究、经济建设和军事领域都具有重要的应用价值。在地震形变提取、地形监测等特定的应用场合中，高分辨率高精度DEM显得尤为重要；但是，DEM提取非常复杂。干涉合成孔径雷达(InSAR，Interferometric Synthetic Aperture Radar)由于其全天时、全天候的工作特性，且作为一种主动式传感器，成为精确获取数字地图的成熟方法之一。但是，在山区或者建筑物较高且分布密集的城区等地表起伏较大的复杂地形区域，在雷达侧视成像获取回波影像对过程中，容易形成叠掩和阴影等数据空洞较多的区域，严重影响了生成的DEM质量。
2007年发射的TerraSAR-X卫星和2010年发射的TanDEM-X卫星，简称TSX/TDX(TerraSAR-X/TanDEM-X)，组成了全球第一个星载双星分布式系统；TSX/TDX的首要目的是获取全球高精度DEM。相对于传统重复轨道卫星，TSX/TDX主要特点是获取“0”时间基线和高精度轨道的干涉对，能够很好的克服大气延迟和轨道误差。在生成的DEM的精度分析中，大多数针对的是陡峭地形的山区或者建筑物分布较为稀疏的居民区；对于高楼林立的密集城区则缺少定量实验和研究，生成的DEM常容易出现透视收缩、叠掩、阴影区域等成像异常区域。同时，已有的数据融合方法多基于TSX/TDX条带模式下的数据，空间分辨率为12米，而这对于城区的地形信息解读还远远不够。
因此，如何提升高落差物体密集区域DEM的成像质量，是亟待解决的问题。
发明内容
有鉴于此，本发明实施例期望提供一种基于双基星载SAR的高精度DEM反演方法及装置，能提升高落差物体密集区域DEM的分辨率和精度，提高成像质量。
为达到上述目的，本发明的技术方案是这样实现的：
本发明实施例提供了一种基于双基星载SAR的高精度DEM反演方法，所述方法包括：
获取第一星载SAR和第二星载SAR在第一轨道对成像目标分别进行成像得到的第一主图像和第一辅图像，并获取第一星载SAR和第二星载SAR在第二轨道对所述成像目标分别进行成像得到的第二主图像和第二辅图像；
采用预设图像处理规则，处理所述第一主图像和第一辅图像，生成所述第一轨道对应的第一干涉图，处理所述第二主图像和第二辅图像，生成所述第二轨道对应的第二干涉图；
根据预设滤波规则，分别对所述第一干涉图和所述第二干涉图进行滤波处理，并确定所述第一干涉图和第二干涉图上各像素对应的相干系数；
根据预设转化规则，将所述滤波后的第一干涉图和第二干涉图，分别转化为第一单基线数字高程模型DEM将和第二单基线DEM；
根据所述第一干涉图和第二干涉图上各像素对应的相干系数，采用预设融合规则，融合所述第一单基线DEM和第二单基线DEM，得到成像目标的融合DEM。
上述方案中，所述采用预设图像处理规则，处理所述第一主图像和第一辅图像，生成所述第一轨道对应的第一干涉图，处理所述第二主图像和第二辅图像，生成所述第二轨道对应的第二干涉图，包括：
对所述第一主图像和第一辅图像进行图像配准，将配准处理后的第一主图像和第一辅图像共轭相乘得到第一干涉图；
对所述第二主图像和第二辅图像进行图像配准，将配准处理后的第二主图像和第二辅图像共轭相乘得到第二干涉图。
上述方案中，所述根据预设滤波规则，分别对所述第一干涉图和第二干涉图进行滤波处理，并确定所述第一干涉图和第二干涉图上各像素对应的相干系数，包括：
采用非局部滤波法，分别对所述第一干涉图和第二干涉图中的各像素进行迭代估计；
根据所述第一干涉图各像素在迭代估计中的相似性权重，采用配准处理后的所述第一主图像和第一辅图像的Goodman模型，得到所述第一干涉图各像素对应的相干系数；
根据所述第二干涉图各像素在迭代估计中的相似性权重，采用配准处理后的所述第二主图像和第二辅图像的Goodman模型，得到所述第二干涉图各像素对应的相干系数。
上述方案中，所述根据预设转化规则，将所述滤波后的第一干涉图和第二干涉图，分别转化为第一单基线DEM将和第二单基线DEM；包括：
将预设DEM仿真干涉相位与第一干涉图的干涉相位进行差分，得到第一差分相位；
对所述第一干涉图中相干系数低于第一预设相干阈值的区域的所述第一差分相位进行掩膜，并将掩膜后的第一差分相位进行解缠；
在解缠的第一差分相位中，将掩膜区域采用预设插值规则进行插值，并将插值后的解缠第一差分相位加上预设DEM仿真干涉相位，得到第一干涉图的解缠相位；
采用干涉合成孔径雷达InSAR测高模型，根据所述第一干涉图的解缠相位，确定所述第一轨道对应的所述第一单基线DEM；
将预设DEM仿真干涉相位与第二干涉图的干涉相位进行差分，得到第二差分相位；
对所述第二干涉图中相干系数低于第一预设相干阈值的区域的所述第二差分相位进行掩膜，并将掩膜后的第二差分相位进行解缠；
在解缠的第二差分相位中，将掩膜区域采用预设插值规则进行插值，并将插值后的解缠第二差分相位加上预设DEM仿真干涉相位，得到第二干涉图的解缠相位；
采用所述InSAR测高模型，根据所述第二干涉图的解缠相位，确定所述第二轨道对应的所述第二单基线DEM。
上述方案中，所述根据所述第一干涉图和第二干涉图上各像素对应的相干系数，采用预设融合规则，融合所述第一单基线DEM和第二单基线DEM，得到成像目标的融合DEM；包括：
分别将所述第一单基线DEM和第二单基线DEM中，成像异常区域和对应相干系数低于第二预设相干阈值区域进行掩膜；
将第一单基线DEM非掩膜区域和第二单基线DEM非掩膜区域进行加权融合，并将第一单基线DEM正常成像区域补充对应位置的所述第二单基线DEM掩膜区域，将第二单基线DEM正常成像区域补充对应位置的所述第一单基线DEM掩膜区域得到成像目标的融合DEM；
将所述第一单基线DEM非掩膜区域对应相干系数占第一单基线DEM非掩膜区域和第二单基线DEM非掩膜区域各自对应相干系数之和的比重，确定为所述第一单基线DEM非掩膜区域进行加权融合的权值；
将所述第二单基线DEM非掩膜区域对应相干系数占第一单基线DEM非掩膜区域和第二单基线DEM非掩膜区域各自对应相干系数之和的比重，确定所述第二单基线DEM非掩膜区域进行加权融合的权值。
上述方案中，所述方法还包括：
根据所述第一干涉图各像素的相干系数，采用相干系数确定规则，确定第一干涉图中各区域的相干系数；
所述第一干涉图中的各区域包含一个以上的像素；
根据述第二干涉图各像素的相干系数，采用相干系数确定规则，确定第二干涉图中各区域的相干系数；
所述第二干涉图中的各区域包含一个以上的像素。
上述方案中，在对所述第一干涉图和第二干涉图进行滤波处理之前，所述方法还包括：
对所述第一干涉图和第二干涉图，分别进行去除平地效应处理。
上述方案中，所述第一星载SAR和第二星载SAR在第一轨道对目标成像时的雷达照射方向，与在第二轨道对目标成像时的雷达照射方向呈预设夹角。
本发明实施例还提供了一种基于双基星载SAR的高精度DEM反演装置，所述装置包括：获取模块、图像处理模块、滤波模块、转化模块和融合模块；其中，
所述获取模块，用于获取第一星载SAR和第二星载SAR在第一轨道对成像目标分别进行分别得到的第一主图像和第一辅图像，并获取第一星载SAR和第二星载SAR在第二轨道对所述成像目标分别进行成像得到的第二主图像和第二辅图像；
所述图像处理模块，用于采用预设图像处理规则，处理所述第一主图像和第一辅图像，生成所述第一轨道对应的第一干涉图，处理所述第二主图像和第二辅图像，生成所述第二轨道对应的第二干涉图；
所述滤波模块，用于根据预设滤波规则，分别对所述第一干涉图和第二干涉图进行滤波处理，并确定所述第一干涉图和所述第二干涉图上各像素对应的相干系数；
所述转化模块，用于根据预设转化规则，将所述滤波后的第一干涉图和第二干涉图，分别转化为第一单基线DEM将和第二单基线DEM；
所述融合模块，用于根据所述第一干涉图和第二干涉图上各像素对应的相干系数，采用预设融合规则，融合所述第一单基线DEM和第二单基线DEM，得到成像目标的融合DEM。
上述方案中，所述图像处理模块，具体用于：
对所述第一主图像和第一辅图像进行图像配准，将配准处理后的第一主图像和第一辅图像共轭相乘得到第一干涉图；
对所述第二主图像和第二辅图像进行图像配准，将配准处理后的第二主图像和第二辅图像共轭相乘得到第二干涉图。
上述方案中，所述滤波模块，具体用于：
采用非局部滤波法，分别对所述第一干涉图和第二干涉图中的各像素进行迭代估计；
根据所述第一干涉图各像素在迭代估计中的相似性权重，采用配准处理后的所述第一主图像和第一辅图像的Goodman模型，得到所述第一干涉图各像素对应的相干系数；
根据所述第二干涉图各像素在迭代估计中的相似性权重，采用配准处理后的所述第二主图像和第二辅图像的Goodman模型，得到所述第二干涉图各像素对应的相干系数。
上述方案中，所述转化模块，具体用于：
将预设DEM仿真干涉相位与第一干涉图的干涉相位进行差分，得到第一差分相位；
对所述第一干涉图中相干系数低于第一预设相干阈值的区域的所述第一差分相位进行掩膜，并将掩膜后的第一差分相位进行解缠；
在解缠的第一差分相位中，将掩膜区域采用预设插值规则进行插值，并将插值后的解缠第一差分相位加上预设DEM仿真干涉相位，得到第一干涉图的解缠相位；
采用InSAR测高模型，根据所述第一干涉图的解缠相位，确定所述第一轨道对应的所述第一单基线DEM；
将预设DEM仿真干涉相位与第二干涉图的干涉相位进行差分，得到第二差分相位；
对所述第二干涉图中相干系数低于第一预设相干阈值的区域的所述第二差分相位进行掩膜，并将掩膜后的第二差分相位进行解缠；
在解缠的第二差分相位中，将掩膜区域采用预设插值规则进行插值，并将插值后的解缠第二差分相位加上预设DEM仿真干涉相位，得到第二干涉图的解缠相位；
采用InSAR测高模型，根据所述第二干涉图的解缠相位，确定所述第二轨道对应的所述第二单基线DEM。
上述方案中，所述融合模块，具体用于：
分别将所述第一单基线DEM和第二单基线DEM中，成像异常区域和对应相干系数低于第二预设相干阈值区域进行掩膜；
将第一单基线DEM非掩膜区域和第二单基线DEM非掩膜区域进行加权融合，并将第一单基线DEM正常成像区域补充对应位置的所述第二单基线DEM掩膜区域，将第二单基线DEM正常成像区域补充对应位置的所述第一单基线DEM掩膜区域得到成像目标的融合DEM；
将所述第一单基线DEM非掩膜区域对应相干系数占第一单基线DEM和第二单基线DEM非掩膜区域各自对应相干系数之和的比重，确定为所述第一单基线DEM非掩膜区域进行加权融合的权值；
将所述第二单基线DEM非掩膜区域对应相干系数占第一单基线DEM和第二单基线DEM非掩膜区域各自对应相干系数之和的比重，确定为所述第二单基线DEM非掩膜区域进行加权融合的权值。
上述方案中，所述滤波模块，还用于：
根据所述第一干涉图各像素的相干系数，采用相干系数确定规则，确定第一干涉图中各区域的相干系数；
所述第一干涉图中的各区域包含一个以上的像素；
根据述第二干涉图各像素的相干系数，采用相干系数确定规则，确定第二干涉图中各区域的相干系数；
所述第二干涉图中的各区域包含一个以上的像素。
上述方案中，所述图像处理模块，还用于：
对所述第一干涉图和第二干涉图，分别进行去除平地效应处理。
上述方案中，所述第一星载SAR和第二星载SAR在第一轨道对目标成像时的雷达照射方向，与在第二轨道对目标成像时的雷达照射方向呈预设夹角。
本发明实施例还提供了一种存储介质，其上存储由可执行程序，所述可执行程序被处理器执行时实现上述方法中任一种所述基于双基星载SAR的高精度DEM反演方法的步骤。
本发明实施例还提供了一种基于双基星载SAR的高精度DEM反演装置，包括处理器、存储器及存储在存储器上并能够有所述处理器运行的可执行程序，所述处理器运行所述可执行程序时执行上述方法中任一种所述基于双基星载SAR的高精度DEM反演方法的步骤。
本发明实施例所提供的基于双基星载SAR的高精度DEM反演方法及装置，获取第一星载SAR和第二星载SAR在第一轨道对成像目标分别进行成像得到的第一主图像和第一辅图像，并获取第一星载SAR和第二星载SAR在第二轨道对所述成像目标分别进行成像得到的第二主图像和第二辅图像；采用预设图像处理规则，处理所述第一主图像和第一辅图像，生成所述第一轨道对应的第一干涉图，处理所述第二主图像和第二辅图像，生成所述第二轨道对应的第二干涉图；根据预设滤波规则，分别对所述第一干涉图和所述第二干涉图进行滤波处理，并确定所述第一干涉图和第二干涉图上各像素对应的相干系数；根据预设转化规则，将所述滤波后的第一干涉图和第二干涉图，分别转化为第一单基线DEM将和第二单基线DEM；根据所述第一干涉图和第二干涉图上各像素对应的相干系数，采用预设融合规则，融合所述第一单基线DEM和第二单基线DEM，得到成像目标的融合DEM。如此，根据两个轨道上的主、辅图像的相干系数，对两个轨道的单基线DEM进行融合，填补单轨成像中的空白区域，能提升高落差物体密集区域DEM的成像质量。
附图说明
图1为本发明实施例基于双基星载SAR的高精度DEM反演方法的流程示意图；
图2为本发明实施例基于双基星载SAR的高精度DEM反演方法的数据处理框架示意图；
图3为本发明实施例基于分块搜索计算的权重示意图；
图4为本发明实施例InSAR测高模型；
图5为本发明实施例升轨配准、去平地后干涉示意图；
图6为本发明实施例轨升轨非局部滤波后的干涉示意图图；
图7为本发明实施例低相干区域掩膜后的升轨DEM示意图；
图8为本发明实施例低相干区域掩膜后的降轨DEM示意图；
图9为本发明实施例融合后的DEM示意图；
图10为本发明实施例基于双基星载SAR的高精度DEM反演装置组成结构示意图。
具体实施方式
本发明实施例中，获取第一星载SAR和第二星载SAR在第一轨道对成像目标分别进行成像得到的第一主图像和第一辅图像，并获取第一星载SAR和第二星载SAR在第二轨道对所述成像目标分别进行成像得到的第二主图像和第二辅图像；采用预设图像处理规则，处理所述第一主图像和第一辅图像，生成所述第一轨道对应的第一干涉图，处理所述第二主图像和第二辅图像，生成所述第二轨道对应的第二干涉图；根据预设滤波规则，分别对所述第一干涉图和所述第二干涉图进行滤波处理，并确定所述第一干涉图和第二干涉图上各像素对应的相干系数；根据预设转化规则，将所述滤波后的第一干涉图和第二干涉图，分别转化为第一单基线DEM将和第二单基线DEM；根据所述第一干涉图和第二干涉图上各像素对应的相干系数，采用预设融合规则，融合所述第一单基线DEM和第二单基线DEM，得到成像目标的融合DEM。
本发明实施例提供的基于双基星载SAR的高精度DEM反演方法，如图1所示，所述方法包括：
步骤110：获取第一星载SAR和第二星载SAR在第一轨道对成像目标分别进行成像得到的第一主图像和第一辅图像，并获取第一星载SAR和第二星载SAR在第二轨道对所述成像目标分别进行成像得到的第二主图像和第二辅图像；
这里，所述第一星载SAR和第二星载SAR可以是双星分布式系统中的两颗卫星的星载SAR，如TerraSAR-X卫星的星载SAR和TanDEM-X卫星的星载SAR；第一卫星和第二卫星可以是双星分布式系统中的主星和辅星；主星和辅星组成的双星系统，以紧密编队方式飞行对地面进行成像，如TerraSAR-X卫星和TanDEM-X卫星之间距离为200～300；所述成像目标可以是城市等；
星载SAR的成像为单视复图像(SLC，Single Look Complex)；所述第一主图像和第二主图像为主星的SAR分别在第一轨道和第二轨道成像得到的SLC；所述第一辅图像和第二辅图像为辅星的SAR分别在第一轨道和第二轨道成像得到的SLC；
如图2所示的DEM数据处理框架示意图中，第一轨道位置相对较高，第二轨道位置相对较低，在第一轨道的成像称为升轨SLC，在第二轨道的成像称为降轨道SLC。主星和辅星分别在第一轨道对目标进行成像得到升轨SLC1和升轨SLC2，分别在第二轨道对所述目标进行成像得到降轨SLC1和降轨SLC2；获取所述SLC后，可以将所有图像传送到主星进行处理，也可以传送到辅星或地面站进行处理。
进一步的，所述第一星载SAR和第二星载SAR在第一轨道对目标成像时的雷达照射方向，与在第二轨道对目标成像时的雷达照射方向呈预设夹角；
所述第一轨道和第二轨道是指在不同时间，双星系统运行的不同轨道高度和位置；由于地球自转和卫星绕地球飞行的原因，双星系统在不同时间点飞越成像目标时轨道高度不同，卫星雷达照射成像目标的方向也不同；如在第一轨道对成像目标成像时，第一星载SAR和第二星载SAR对目标成像时的雷达照射方向可以是由东向西；在第二轨道对成像目标成像时，第一星载SAR和第二星载SAR对目标成像时的雷达照射方向可以是由西向东。如此，可以在两个方位对成像目标进行成像，避免因为轨道和成像目标的夹角引起的仅能在山峰、高楼等迎坡面成像，增加了背坡面成像的几率。
步骤120：采用预设图像处理规则，处理所述第一主图像和第一辅图像，生成所述第一轨道对应的第一干涉图，处理所述第二主图像和第二辅图像，生成所述第二轨道对应的第二干涉图；
这里，所述预设图像处理规则可以根据第一主图像和第一辅图像及第二主图像和第二辅图像的成像条件，如卫星与成像目标位置关系等，以及成像大小等进行设置；可以对主、辅图像进行图像配准后，再生成干涉图。可以采用最大干涉频谱法和相位差影像平均波动函数法等方法进行配准；可以将配准后的第一主图像和第一辅图像共轭相乘得到第一轨道对应的第一干涉图，将配准后的第二主图像和第二辅图像共轭相乘得到第二轨道对应的第二干涉图。
进一步的，对所述第一主图像和第一辅图像进行图像配准，将配准处理后的第一主图像和第一辅图像共轭相乘得到第一干涉图；对所述第二主图像和第二辅图像进行图像配准，将配准处理后的第二主图像和第二辅图像共轭相乘得到第二干涉图；这里，可以采用相关系数法对第一主图像和第一辅图像进行配准，及对第二主图像和第二辅图像进行配准，再生成第一干涉图和第二干涉图；
具体的，相关系数法采用基于窗口的配准方法，分别在主图像和辅图像上分别建立匹配窗和搜索窗进行配准。在主图像和辅图像中选取一些控制点，建立的匹配窗和搜索窗，以控制点为中心。滑动匹配窗求得匹配窗和搜索窗的窗口的复相关函数，进而得到偏移量；所述复相关函数可以用表达式(1)表示：
Figure GDA0002221819780000111
其中，ρc表示复相关函数，s1和s2表示主图像和辅图像，(u,v)表示滑动的匹配窗窗口中心点坐标，(m,n)表示像素中心点坐标，即控制点坐标，M、N是图像的尺寸大小，*表示取复共轭。|□|表示取幅度值，|□|中的“□”可以表示任意表达式。所述匹配窗的大小可以根据处理速度和精度确定，如20x20等；
复相关函数ρc的最大值处即为最佳匹配点，可以对复相关函数的最大值处进行插值，如采用0.5个像素为单元插值，获得更加精确的偏移量，根据控制点对应的配准偏移量获得其他像素处的配准偏移量，可以利用二阶多项式(2)来求取：
Figure GDA0002221819780000121
其中，(x,y)表示主图像中的像素坐标，(Δx,Δy)表示辅图像相对于主图像的偏移量，a和b表示拟合系数，可以根据公式1得到的偏移量拟合得到，其中，a0和b0表示常驻系数，下标表示维数和阶数；a10表示1维0阶系数，a21表示2维1阶系数，以此类推，在此不再赘述。
获取配准后的第一主图像和第一辅图像，以及第二主图像和第二辅图像，可以将配准第一主图像和第一辅图像进行共轭相乘，得到第一轨道对应的第一干涉图；将配准第二主图像和第二辅图像进行共轭相乘，得到第二轨道对应的第二干涉图。
更进一步的，在对所述第一干涉图和第二干涉图进行滤波处理之前，所述方法还包括：对所述第一干涉图和第二干涉图，分别进行去除平地效应处理；
这里，可以采用轨道法去除平地效应；以TSX/TDX为例，TSX/TDX具有精确的轨道参数等信息，可以根据轨道参数和场景地面入射角，计算平地相位，通过和干涉图共轭相乘去除平地效应。
步骤130：根据预设滤波规则，分别对所述第一干涉图和所述第二干涉图进行滤波处理，并确定所述第一干涉图和第二干涉图上各像素对应的相干系数；
由于相干成像系统的固有特性及星载SAR的系统噪声等，第一干涉图和第二干涉图的存在斑点噪声、系统噪声等不同噪声；因此，需要对干涉图进行滤波处理；可以采用均值滤波、Goldstein滤波等方式对所述第一干涉图和第二干涉图进行滤波；
所述第一干涉图上各像素对应的相干系数是指反应所述第一干涉图上各像素对应的配准后的第一主图像和第一辅图像各像素之间的相干性关系的相干系数；所述第二干涉图上各像素对应的相干系数是指反应所述第二干涉图上各像素对应的配准后的第二主图像和第二辅图像各像素之间的相干性关系的相干系数；所述第一干涉图上各像素对应的相干系数可以通过配准后的第一主图像和第一辅图像对应位置的相邻像素采用多视处理等方法得到；所述第二干涉图上各像素对应的相干系数可以通过配准后的第二主图像和第二辅图像对应位置的相邻像素采用多视处理等方法得到。
进一步的，采用非局部滤波法，分别对所述第一干涉图和第二干涉图中的各像素进行迭代估计；根据所述第一干涉图各像素在迭代估计中的相似性权重，采用配准处理后的所述第一主图像和第一辅图像的古德曼Goodman模型，确定所述第一干涉图各像素对应的相干系数；根据所述第二干涉图各像素在迭代估计中的相似性权重，采用配准处理后的所述第二主图像和第二辅图像的Goodman模型，确定所述第二干涉图各像素对应的相干系数；
具体的，两幅配准的SLC复图像每一点的像素模型可以用三个参数来描述：反射系数、真实的干涉相位、相干系数，假设z和z′是两幅配准图像上对应的像素点，根据Goodman模型，z与z′服从表达式(3)：
Figure GDA0002221819780000131
其中，∑表示2×2的相干系数矩阵，可以用表达式(4)表示：
Figure GDA0002221819780000132
其中，R与R′表示反射系数，β表示真实相位，D表示相干系数，E表示取期望值；
可以采用非局部滤波对三个参数R、β、D联合估计，假设两幅图像的散射特性相同，即R＝R′，这种假设在相干性较高的区域是成立的，同时令A＝|z|、A′＝|z′|；φ＝arg(zz′*)为去平地后的相位，则有表达式(5)：
Figure GDA0002221819780000133
非局部性干涉合成孔径雷达(NL-InSAR，Non Local Interferometric SyntheticAperture Radar)使用权重最大似然估计(WMLE，Weighted Maximum LikelihoodEstimation)来估计权重，在分块s中，参数估计可以用表达式(6)表示：
Figure GDA0002221819780000141
令R＝R′，将公式扩展到权重最大似然估计中，令Θs＝(Rs,βs,Ds)以及
Figure GDA0002221819780000142
估计的R、β、D可以用表达式(7)表示：
Figure GDA0002221819780000143
其中，
Figure GDA0002221819780000144
w(s,t)是基于即采用迭代估计的方法实现干涉图的滤波中分块搜索计算的相似性权重；采用迭代估计的方法实现干涉图的滤波的计算方式可以如图3所示：
假设图像I上的待估计像素点s可以用表达式(8)表示：
v＝{vs|s∈I} (8)
判断待估计像素点s，和相似像素点t是否相似，是以s和t为中心选取一定大小的窗口，判断以s，t为中心的这个窗口是否相似；将以像素点s为中心像素窗口作为目标窗31，以像素点t为中心的像素窗口作为相似窗32，在选定搜索窗33内，以预设的路线34进行搜索，确定各相似窗32与目标窗31的相似性，像素点s的估计值可以用表达式(9)表示：
Figure GDA0002221819780000145
其中，vs表示待估计的像素点，vt表示相似窗中心像素；权重w(s,t)的选择是基于s、t像素周围的像素窗口的相似性，满足0≤w(s,t)≤1，并且
Figure GDA0002221819780000146
所述目标窗和相似窗的大小可以根据处理速度和精度确定，如3X3等；
获得了各像素德点估计值，即完成了滤波，这里，对当前像素进行迭代估计就是一个滤波处理过程。
使用NL-InSAR在滤除相位噪声，降噪性能效果好，利于后续的相位解缠。
进一步的，根据述第一干涉图各像素的相干系数，采用相干系数确定规则，确定第一干涉图中各区域的相干系数；所述第一干涉图中的各区域包含一个以上的像素；根据述第二干涉图各像素的相干系数，采用相干系数确定规则，确定第二干涉图中各区域的相干系数；所述第二干涉图中的各区域包含一个以上的像素；
为了提高后续的计算速度，可以针对干涉图中不同区域的相干系数，所述区域可以是干涉图中的一个像素，也可以预设区域大小；可以根据干涉图精度需求来设置区域的大小；采用相干系数确定规则，可以根据区域与所述区域包含像素关系确定，如所述区域为一个像素时，改区域相干系数为该像素的相干系数，所述区域大于一个像素时，所述区域的相干系数可以将区域中各像素相干系数做算数平均等运算确定。
步骤140：根据预设转化规则，将所述滤波后的第一干涉图和第二干涉图，分别转化为第一单基线DEM将和第二单基线DEM；
这里，可以根据干涉图和星载SAR的姿态数据等生成DEM，根据所述第一干涉图和第二干涉图可以分别生成一个与轨道对应的单基线DEM；所述预设转化规则可以根据对第一星载SAR和第二星载SAR的轨道位置、DEM的精度需求、运算时间需求等设置；可以采用常规方式生成DEM。
进一步的，可以将预设DEM仿真干涉相位与第一干涉图的干涉相位进行差分，得到第一差分相位；对所述第一干涉图中相干系数低于第一预设相干阈值的区域的所述第一差分相位进行掩膜，并将掩膜后的第一差分相位进行解缠；在解缠的第一差分相位中，将掩膜区域采用预设插值规则进行插值，并将插值后的解缠第一差分相位加上预设DEM仿真干涉相位，得到第一干涉图的解缠相位；采用InSAR测高模型，根据所述第一干涉图的解缠相位，确定所述第一轨道对应的所述第一单基线DEM；将预设DEM仿真干涉相位与第二干涉图的干涉相位进行差分，得到第二差分相位；对所述第二干涉图中相干系数低于第一预设相干阈值的区域的所述第二差分相位进行掩膜，并将掩膜后的第二差分相位进行解缠；在解缠的第二差分相位中，将掩膜区域采用预设插值规则进行插值，并将插值后的解缠第二差分相位加上预设DEM仿真干涉相位，得到第二干涉图的解缠相位；采用所述InSAR测高模型，根据所述第二干涉图的解缠相位，确定所述第二轨道对应的所述第二单基线DEM；
这里，所述预设DEM可以是现有的成像目标的DEM，根据第一星载SAR和第二星载SAR的SAR系统参数，以及第一星载SAR和第二星载SAR的轨道参数，通过仿真的到所述预设DEM的干涉相位；采用预设DEM的仿真干涉相位与干涉图进行差分干涉处理得到缠绕的差分相位；为了保证相位解缠的可靠性，可以将干涉图中低相干区域掩膜掉，再解缠差分相位；这里可以预设一个第一预设相干阈值，将干涉图各区域中相干系数低于第一预设相干阈值的区域进行掩膜，对掩膜区域进行插值；将掩膜处理后的解缠差分相位加上预设DEM的仿真干涉相位，得到干涉图的解缠相位，这里，可以对第一干涉图和第二干涉图均进行上述处理，分别得到第一干涉图和第二干涉图的解缠相位；
得到解缠相位后，可以采用InSAR几何关系模型反演高度，InSAR数字高程重建的空间几何关系如图4所示，其中：A1和A2分别分别为第一星载SAR和第二星载SAR；H表示平台高度，A1和A2分别表示两个天线的相位中心，α为基线角，B为干涉基线的长度，P为地面上的某一个目标点，r1和r2＝r1+Δr分别为零多普勒平面内两个天线的相位中心到目标点P的距离，θ为主天线相位中心到目标的视角，地形高程用h表示，yP为P点的地距，y为地距向坐标。根据这个空间几何关系可以推导出地形的高程h可以用表达式(10)表示：
Figure GDA0002221819780000161
通过第一干涉图和第二干涉图的解缠相位分别得到第一干涉图和第二干涉图对应的高程数据后，根据成像目标的地理一起进行编码，得到第一轨道和第二轨道分别对应的第一单基线DEM和第二单基线DEM。
这里，使用外部预设DEM进行差分处理，可以降低相位解缠的难度，保证相位解缠的可靠性。
步骤150：根据所述第一干涉图和第二干涉图上各像素对应的相干系数，采用预设融合规则，融合所述第一单基线DEM和第二单基线DEM，得到成像目标的融合DEM；
所述第一单基线DEM和第二单基线DEM中有较多的数据空洞；无法反应真是的地形特征；这里，可以对第一单基线DEM和第二单基线DEM进行融合；所述预设融合规则可以根据基线DEM中不同区域进行设置，如对于成像遮挡区域等进行互相补充，对正常成像区域进行加权融合等；
进一步的，分别将所述第一单基线DEM和第二单基线DEM中，成像异常区域和对应相干系数低于第二预设相干阈值区域进行掩膜；将第一单基线DEM非掩膜区域和第二单基线DEM非掩膜区域进行加权融合，并将第一单基线DEM正常成像区域补充对应位置的所述第二单基线DEM掩膜区域，将第二单基线DEM正常成像区域补充对应位置的所述第一单基线DEM掩膜区域得到成像目标的融合DEM；将所述第一单基线DEM非掩膜区域对应相干系数占第一单基线DEM非掩膜区域和第二单基线DEM非掩膜区域各自对应相干系数之和的比重，确定为所述第一单基线DEM非掩膜区域进行加权融合的权值；将所述第二单基线DEM非掩膜区域对应相干系数占第一单基线DEM非掩膜区域和第二单基线DEM非掩膜区域各自对应相干系数之和的比重，确定为所述第二单基线DEM非掩膜区域进行加权融合的权值；
具体的，所述成像异常区域可以是单基线DEM中的透视收缩、叠掩、阴影区域；所述对应位置，是指相同地理坐标的位置。可以将第一单基线DEM和第二单基线DEM的透视收缩、叠掩、阴影进行掩膜，同时，采用NL滤波得到的无偏估计的相干系数，将相干系数低于第二预设相干阈值，即相干性较低的区域也进行掩膜；然后利用第一轨道和第二轨道相反照射方向的特点，将空值区域互相补充；将非空值区域加权融合，得到融合后的DEM。所述加权融合中的加权值可以根据第一单基线DEM和第二单基线DEM各像素或区域对应的相干系数确定，如第一单基线DEM各像素或区域对应的相干系数为DA，第二单基线DEM各像素或区域对应的相干系数为DB，那么，所述加权融合中的第一单基线DEM加权值则为第二单基线DEM加权值则为其中，为了简化运算，所述第一预设相干阈值和第二预设相干阈值可以取同一个值；
如此，利用第一轨道和第二轨道不同轨道方向上获得的影像对叠掩等“空白”区域进行互补，如利用降轨DEM上背坡面的区域填补升轨DEM上迎坡面的空白区域，重建高落差物体如建筑物密集区域，细节特征丰富的高分辨率高精度城区DEM。并且由于本发明实施例直接使用一次多视处理操作，避免了多次重复迭代处理，处理起来简单快速，同时避免多次迭代时的误差传递问题。
下面结合具体示例对本发明产生的积极效果作进一步详细的描述；
本示例采用TSX/TDX升降轨图像数据验证本发明实施例基于双基星载SAR的高精度DEM反演方法。图5为两幅升轨SLC配准、去平地之后的干涉图；图6为采用非局部均值滤波得到的干涉图，可以发现改图像噪声得到了很好的去除，并且建筑物边缘轮廓清晰。图7和图8分别为采用NL-InSAR得到的分辨率为4米的升轨DEM和降轨DEM，无效区域用空值表示。图9为升降轨融合后的DEM，可以看出，相比于升轨或者降轨DEM，融合后的DEM无效区域明显减少，统计结果表明无效点数比例由升轨、降轨的4.93％和4.52％降低到1.34％，同时，与SRTM数据作比较，融合之后均方根误差(RMSE，Root Mean Square Error)由升轨的6.37米、降轨的5.61米降低到5.06米。
本发明实施例提供的基于双基星载SAR的高精度DEM反演装置，如图10所示，所述装置包括：获取模块101、图像处理模块102、滤波模块103、转化模块104和融合模块105；其中，
所述获取模块101，用于获取第一星载SAR和第二星载SAR在第一轨道对成像目标分别进行成像得到的第一主图像和第一辅图像，并获取第一星载SAR和第二星载SAR在第二轨道对所述成像目标分别进行成像得到的第二主图像和第二辅图像；
这里，所述第一星载SAR和第二星载SAR可以是双星分布式系统中的两颗卫星的星载SAR，如TerraSAR-X卫星的星载SAR和TanDEM-X卫星的星载SAR；第一卫星和第二卫星可以是双星分布式系统中的主星和辅星；主星和辅星组成的双星系统，以紧密编队方式飞行对地面进行成像，如TerraSAR-X卫星和TanDEM-X卫星之间距离为200～300；所述成像目标可以是城市等；
星载SAR的成像为SLC；所述第一主图像和第二主图像为主星的SAR分别在第一轨道和第二轨道成像得到的SLC；所述第一辅图像和第二辅图像为辅星的SAR分别在第一轨道和第二轨道成像得到的SLC；
如图2所示的DEM数据处理框架示意图中，第一轨道位置相对较高，第二轨道位置相对较低，在第一轨道的成像称为升轨SLC，在第二轨道的成像称为降轨道SLC。主星和辅星分别在第一轨道对目标进行成像得到升轨SLC1和升轨SLC2，分别在第二轨道对所述目标进行成像得到降轨SLC1和降轨SLC2；获取所述SLC后，可以将所有图像传送到主星进行处理，也可以传送到辅星或地面站进行处理。
进一步的，所述第一星载SAR和第二星载SAR在第一轨道对目标成像时的雷达照射方向，与在第二轨道对目标成像时的雷达照射方向呈预设夹角；
所述第一轨道和第二轨道是指在不同时间，双星系统运行的不同轨道高度和位置；由于地球自转和卫星绕地球飞行的原因，双星系统在不同时间点飞越成像目标时轨道高度不同，卫星雷达照射成像目标的方向也不同；如在第一轨道对成像目标成像时，第一星载SAR和第二星载SAR对目标成像时的雷达照射方向可以是由东向西；在第二轨道对成像目标成像时，第一星载SAR和第二星载SAR对目标成像时的雷达照射方向可以是由西向东。如此，可以在两个方位对成像目标进行成像，避免因为轨道和成像目标的夹角引起的仅能在山峰、高楼等迎坡面成像，增加了背坡面成像的几率。
所述图像处理模块102，用于采用预设图像处理规则，处理所述第一主图像和第一辅图像，生成所述第一轨道对应的第一干涉图，处理所述第二主图像和第二辅图像，生成所述第二轨道对应的第二干涉图；
这里，所述预设图像处理规则可以根据第一主图像和第一辅图像及第二主图像和第二辅图像的成像条件，如卫星与成像目标位置关系等，以及成像大小等进行设置；可以对主、辅图像进行图像配准后，再生成干涉图。可以采用最大干涉频谱法和相位差影像平均波动函数法等方法进行配准；可以将配准后的第一主图像和第一辅图像共轭相乘得到第一轨道对应的第一干涉图，将配准后的第二主图像和第二辅图像共轭相乘得到第二轨道对应的第二干涉图。
进一步的，对所述第一主图像和第一辅图像进行图像配准，将配准处理后的第一主图像和第一辅图像共轭相乘得到第一干涉图；对所述第二主图像和第二辅图像进行图像配准，将配准处理后的第二主图像和第二辅图像共轭相乘得到第二干涉图；这里，可以采用相关系数法对第一主图像和第一辅图像进行配准，及对第二主图像和第二辅图像进行配准，再生成第一干涉图和第二干涉图；
具体的，相关系数法采用基于窗口的配准方法，分别在主图像和辅图像上分别建立匹配窗和搜索窗进行配准。在主图像和辅图像中选取一些控制点，建立的匹配窗和搜索窗，以控制点为中心。滑动匹配窗求得匹配窗和搜索窗的窗口的复相关函数，进而得到偏移量；所述复相关函数可以用表达式(1)表示；
其中，ρc表示复相关函数，s1和s2表示主图像和辅图像，(u,v)表示滑动的匹配窗窗口中心点坐标，(m,n)表示像素中心点坐标，即控制点坐标，M、N是图像的尺寸大小，*表示取复共轭。|□|表示取幅度值，|□|中的“□”可以表示任意表达式。所述匹配窗的大小可以根据处理速度和精度确定，如20x20等；
复相关函数ρc的最大值处即为最佳匹配点，可以对复相关函数的最大值处进行插值，如采用0.5个像素为单元插值，获得更加精确的偏移量，根据控制点对应的配准偏移量获得其他像素处的配准偏移量，可以利用二阶多项式(2)来求取；其中，(x,y)表示主图像中的像素坐标，(Δx,Δy)表示辅图像相对于主图像的偏移量，a和b表示拟合系数，可以根据公式1得到的偏移量拟合得到，其中，a0和b0表示常驻系数，下标表示维数和阶数；a10表示1维0阶系数，a21表示2维1阶系数，以此类推，在此不再赘述。
获取配准后的第一主图像和第一辅图像，以及第二主图像和第二辅图像，可以将配准第一主图像和第一辅图像进行共轭相乘，得到第一轨道对应的第一干涉图；将配准第二主图像和第二辅图像进行共轭相乘，得到第二轨道对应的第二干涉图。
更进一步的，在对所述第一干涉图和第二干涉图进行滤波处理之前，所述方法还包括：对所述第一干涉图和第二干涉图，分别进行去除平地效应处理；
这里，可以采用轨道法去除平地效应；以TSX/TDX为例，TSX/TDX具有精确的轨道参数等信息，可以根据轨道参数和场景地面入射角，计算平地相位，通过和干涉图共轭相乘去除平地效应。
所述滤波模块103，用于根据预设滤波规则，分别对所述第一干涉图和所述第二干涉图进行滤波处理，并确定所述第一干涉图和第二干涉图上各像素对应的相干系数；
由于相干成像系统的固有特性及星载SAR的系统噪声等，第一干涉图和第二干涉图的存在斑点噪声、系统噪声等不同噪声；因此，需要对干涉图进行滤波处理；可以采用均值滤波、Goldstein滤波等方式对所述第一干涉图和第二干涉图进行滤波；
所述第一干涉图上各像素对应的相干系数是指反应所述第一干涉图上各像素对应的配准后的第一主图像和第一辅图像各像素之间的相干性关系的相干系数；所述第二干涉图上各像素对应的相干系数是指反应所述第二干涉图上各像素对应的配准后的第二主图像和第二辅图像各像素之间的相干性关系的相干系数；所述第一干涉图上各像素对应的相干系数可以通过配准后的第一主图像和第一辅图像对应位置的相邻像素采用多视处理等方法得到；所述第二干涉图上各像素对应的相干系数可以通过配准后的第二主图像和第二辅图像对应位置的相邻像素采用多视处理等方法得到。
进一步的，采用非局部滤波法，分别对所述第一干涉图和第二干涉图中的各像素进行迭代估计；根据所述第一干涉图各像素在迭代估计中的相似性权重，采用配准处理后的所述第一主图像和第一辅图像的古德曼Goodman模型，确定所述第一干涉图各像素对应的相干系数；根据所述第二干涉图各像素在迭代估计中的相似性权重，采用配准处理后的所述第二主图像和第二辅图像的Goodman模型，确定所述第二干涉图各像素对应的相干系数；
具体的，两幅配准的SLC复图像每一点的像素模型可以用三个参数来描述：反射系数、真实的干涉相位、相干系数，假设z和z′是两幅配准图像上对应的像素点，根据Goodman模型，z与z′服从表达式(3)；其中，∑表示2×2的相干系数矩阵，可以用表达式(4)表示；其中，R与R′表示反射系数，β表示真实相位，D表示相干系数，E表示取期望值；
可以采用非局部滤波对三个参数R、β、D联合估计，假设两幅图像的散射特性相同，即R＝R′，这种假设在相干性较高的区域是成立的，同时令A＝|z|、A′＝|z′|；φ＝arg(zz′*)为去平地后的相位，则有表达式(5)；
NL-InSAR使用WMLE来估计权重，在分块s中，参数估计可以用表达式(6)表示；
令R＝R′，将公式扩展到权重最大似然估计中，令Θs＝(Rs,βs,Ds)以及
Figure GDA0002221819780000221
估计的R、β、D可以用表达式(7)表示；其中， w(s,t)是基于即采用迭代估计的方法实现干涉图的滤波中分块搜索计算的相似性权重；采用迭代估计的方法实现干涉图的滤波的计算方式可以如图3所示；假设图像I上的待估计像素点s可以用表达式(8)表示；判断待估计像素点s，和相似像素点t是否相似，是以s和t为中心选取一定大小的窗口，判断以s，t为中心的这个窗口是否相似；将以像素点s为中心像素窗口作为目标窗31，以像素点t为中心的像素窗口作为相似窗32，在选定搜索窗33内，以预设的路线34进行搜索，确定各相似窗32与目标窗31的相似性，像素点s的估计值可以用表达式(9)表示；其中，vs表示待估计的像素点，vt表示相似窗中心像素；权重w(s,t)的选择是基于s、t像素周围的像素窗口的相似性，满足0≤w(s,t)≤1，并且
Figure GDA0002221819780000224
所述目标窗和相似窗的大小可以根据处理速度和精度确定，如3X3等；
获得了各像素德点估计值，即完成了滤波，这里，对当前像素进行迭代估计就是一个滤波处理过程。
使用NL-InSAR在滤除相位噪声，降噪性能效果好，利于后续的相位解缠。
进一步的，根据述第一干涉图各像素的相干系数，采用相干系数确定规则，确定第一干涉图中各区域的相干系数；所述第一干涉图中的各区域包含一个以上的像素；根据述第二干涉图各像素的相干系数，采用相干系数确定规则，确定第二干涉图中各区域的相干系数；所述第二干涉图中的各区域包含一个以上的像素；
为了提高后续的计算速度，可以针对干涉图中不同区域的相干系数，所述区域可以是干涉图中的一个像素，也可以预设区域大小；可以根据干涉图精度需求来设置区域的大小；采用相干系数确定规则，可以根据区域与所述区域包含像素关系确定，如所述区域为一个像素时，改区域相干系数为该像素的相干系数，所述区域大于一个像素时，所述区域的相干系数可以将区域中各像素相干系数做算数平均等运算确定。
所述转化模块104，用于根据预设转化规则，将所述滤波后的第一干涉图和第二干涉图，分别转化为第一单基线DEM将和第二单基线DEM；
这里，可以根据干涉图和星载SAR的姿态数据等生成DEM，根据所述第一干涉图和第二干涉图可以分别生成一个与轨道对应的单基线DEM；所述预设转化规则可以根据对第一星载SAR和第二星载SAR的轨道位置、DEM的精度需求、运算时间需求等设置；可以采用常规方式生成DEM。
进一步的，可以将预设DEM仿真干涉相位与第一干涉图的干涉相位进行差分，得到第一差分相位；对所述第一干涉图中相干系数低于第一预设相干阈值的区域的所述第一差分相位进行掩膜，并将掩膜后的第一差分相位进行解缠；在解缠的第一差分相位中，将掩膜区域采用预设插值规则进行插值，并将插值后的解缠第一差分相位加上预设DEM仿真干涉相位，得到第一干涉图的解缠相位；采用InSAR测高模型，根据所述第一干涉图的解缠相位，确定所述第一轨道对应的所述第一单基线DEM；将预设DEM仿真干涉相位与第二干涉图的干涉相位进行差分，得到第二差分相位；对所述第二干涉图中相干系数低于第一预设相干阈值的区域的所述第二差分相位进行掩膜，并将掩膜后的第二差分相位进行解缠；在解缠的第二差分相位中，将掩膜区域采用预设插值规则进行插值，并将插值后的解缠第二差分相位加上预设DEM仿真干涉相位，得到第二干涉图的解缠相位；采用所述InSAR测高模型，根据所述第二干涉图的解缠相位，确定所述第二轨道对应的所述第二单基线DEM；
这里，所述预设DEM可以是现有的成像目标的DEM，根据第一星载SAR和第二星载SAR的SAR系统参数，以及第一星载SAR和第二星载SAR的轨道参数，通过仿真的到所述预设DEM的干涉相位；采用预设DEM的仿真干涉相位与干涉图进行差分干涉处理得到缠绕的差分相位；为了保证相位解缠的可靠性，可以将干涉图中低相干区域掩膜掉，再解缠差分相位；这里可以预设一个第一预设相干阈值，将干涉图各区域中相干系数低于第一预设相干阈值的区域进行掩膜，对掩膜区域进行插值；将掩膜处理后的解缠差分相位加上预设DEM的仿真干涉相位，得到干涉图的解缠相位，这里，可以对第一干涉图和第二干涉图均进行上述处理，分别得到第一干涉图和第二干涉图的解缠相位；
得到解缠相位后，可以采用InSAR几何关系模型反演高度，InSAR数字高程重建的空间几何关系如图4所示，其中：A1和A2分别分别为第一星载SAR和第二星载SAR；H表示平台高度，A1和A2分别表示两个天线的相位中心，α为基线角，B为干涉基线的长度，P为地面上的某一个目标点，r1和r2＝r1+Δr分别为零多普勒平面内两个天线的相位中心到目标点P的距离，θ为主天线相位中心到目标的视角，地形高程用h表示，yP为P点的地距，y为地距向坐标。根据这个空间几何关系可以推导出地形的高程h可以用表达式(10)表示；
通过第一干涉图和第二干涉图的解缠相位分别得到第一干涉图和第二干涉图对应的高程数据后，根据成像目标的地理一起进行编码，得到第一轨道和第二轨道分别对应的第一单基线DEM和第二单基线DEM。
这里，使用外部预设DEM进行差分处理，可以降低相位解缠的难度，保证相位解缠的可靠性。
所述融合模块105，用于：根据所述第一干涉图和第二干涉图上各像素对应的相干系数，采用预设融合规则，融合所述第一单基线DEM和第二单基线DEM，得到成像目标的融合DEM；
所述第一单基线DEM和第二单基线DEM中有较多的数据空洞；无法反应真是的地形特征；这里，可以对第一单基线DEM和第二单基线DEM进行融合；所述预设融合规则可以根据基线DEM中不同区域进行设置，如对于成像遮挡区域等进行互相补充，对正常成像区域进行加权融合等；
进一步的，分别将所述第一单基线DEM和第二单基线DEM中，成像异常区域和对应相干系数低于第二预设相干阈值区域进行掩膜；将第一单基线DEM非掩膜区域和第二单基线DEM非掩膜区域进行加权融合，并将第一单基线DEM正常成像区域补充对应位置的所述第二单基线DEM掩膜区域，将第二单基线DEM正常成像区域补充对应位置的所述第一单基线DEM掩膜区域得到成像目标的融合DEM；将所述第一单基线DEM非掩膜区域对应相干系数占第一单基线DEM非掩膜区域和第二单基线DEM非掩膜区域各自对应相干系数之和的比重，确定为所述第一单基线DEM非掩膜区域进行加权融合的权值；将所述第二单基线DEM非掩膜区域对应相干系数占第一单基线DEM非掩膜区域和第二单基线DEM非掩膜区域各自对应相干系数之和的比重，确定为所述第二单基线DEM非掩膜区域进行加权融合的权值；
具体的，所述成像异常区域可以是单基线DEM中的透视收缩、叠掩、阴影区域；所述对应位置，是指相同地理坐标的位置。可以将第一单基线DEM和第二单基线DEM的透视收缩、叠掩、阴影进行掩膜，同时，采用NL滤波得到的无偏估计的相干系数，将相干系数低于第二预设相干阈值，即相干性较低的区域也进行掩膜；然后利用第一轨道和第二轨道相反照射方向的特点，将空值区域互相补充；将非空值区域加权融合，得到融合后的DEM。所述加权融合中的加权值可以根据第一单基线DEM和第二单基线DEM各像素或区域对应的相干系数确定，如第一单基线DEM各像素或区域对应的相干系数为DA，第二单基线DEM各像素或区域对应的相干系数为DB，那么，所述加权融合中的第一单基线DEM加权值则为
Figure GDA0002221819780000261
第二单基线DEM加权值则为
Figure GDA0002221819780000262
其中，为了简化运算，所述第一预设相干阈值和第二预设相干阈值可以取同一个值；
如此，利用第一轨道和第二轨道不同轨道方向上获得的影像对叠掩等“空白”区域进行互补，如利用降轨DEM上背坡面的区域填补升轨DEM上迎坡面的空白区域，重建高落差物体如建筑物密集区域，细节特征丰富的高分辨率高精度城区DEM。并且由于本发明实施例直接使用一次多视处理操作，避免了多次重复迭代处理，处理起来简单快速，同时避免多次迭代时的误差传递问题。
在实际应用中，所述获取模块101、图像处理模块102、滤波模块103、转化模块104和融合模块105均可以由卫星系统或地面系统中的中央处理器(CPU)、微处理器(MPU)、数字信号处理器(DSP)、或现场可编程门阵列(FPGA)等实现。
本发明实施例提供的存储介质，其上存储由可执行程序，所述可执行程序被处理器执行时实现基于双基星载SAR的高精度DEM反演方法的步骤；如图1所示，所述基于双基星载SAR的高精度DEM反演方法的步骤包括：
步骤110：获取第一星载SAR和第二星载SAR在第一轨道对成像目标分别进行成像得到的第一主图像和第一辅图像，并获取第一星载SAR和第二星载SAR在第二轨道对所述成像目标分别进行成像得到的第二主图像和第二辅图像；
这里，所述第一星载SAR和第二星载SAR可以是双星分布式系统中的两颗卫星的星载SAR，如TerraSAR-X卫星的星载SAR和TanDEM-X卫星的星载SAR；第一卫星和第二卫星可以是双星分布式系统中的主星和辅星；主星和辅星组成的双星系统，以紧密编队方式飞行对地面进行成像，如TerraSAR-X卫星和TanDEM-X卫星之间距离为200～300；所述成像目标可以是城市等；
星载SAR的成像为SLC；所述第一主图像和第二主图像为主星的SAR分别在第一轨道和第二轨道成像得到的SLC；所述第一辅图像和第二辅图像为辅星的SAR分别在第一轨道和第二轨道成像得到的SLC；
如图2所示的DEM数据处理框架示意图中，第一轨道位置相对较高，第二轨道位置相对较低，在第一轨道的成像称为升轨SLC，在第二轨道的成像称为降轨道SLC。主星和辅星分别在第一轨道对目标进行成像得到升轨SLC1和升轨SLC2，分别在第二轨道对所述目标进行成像得到降轨SLC1和降轨SLC2；获取所述SLC后，可以将所有图像传送到主星进行处理，也可以传送到辅星或地面站进行处理。
进一步的，所述第一星载SAR和第二星载SAR在第一轨道对目标成像时的雷达照射方向，与在第二轨道对目标成像时的雷达照射方向呈预设夹角；
所述第一轨道和第二轨道是指在不同时间，双星系统运行的不同轨道高度和位置；由于地球自转和卫星绕地球飞行的原因，双星系统在不同时间点飞越成像目标时轨道高度不同，卫星雷达照射成像目标的方向也不同；如在第一轨道对成像目标成像时，第一星载SAR和第二星载SAR对目标成像时的雷达照射方向可以是由东向西；在第二轨道对成像目标成像时，第一星载SAR和第二星载SAR对目标成像时的雷达照射方向可以是由西向东。如此，可以在两个方位对成像目标进行成像，避免因为轨道和成像目标的夹角引起的仅能在山峰、高楼等迎坡面成像，增加了背坡面成像的几率。
步骤120：采用预设图像处理规则，处理所述第一主图像和第一辅图像，生成所述第一轨道对应的第一干涉图，处理所述第二主图像和第二辅图像，生成所述第二轨道对应的第二干涉图；
这里，所述预设图像处理规则可以根据第一主图像和第一辅图像及第二主图像和第二辅图像的成像条件，如卫星与成像目标位置关系等，以及成像大小等进行设置；可以对主、辅图像进行图像配准后，再生成干涉图。可以采用最大干涉频谱法和相位差影像平均波动函数法等方法进行配准；可以将配准后的第一主图像和第一辅图像共轭相乘得到第一轨道对应的第一干涉图，将配准后的第二主图像和第二辅图像共轭相乘得到第二轨道对应的第二干涉图。
进一步的，对所述第一主图像和第一辅图像进行图像配准，将配准处理后的第一主图像和第一辅图像共轭相乘得到第一干涉图；对所述第二主图像和第二辅图像进行图像配准，将配准处理后的第二主图像和第二辅图像共轭相乘得到第二干涉图；这里，可以采用相关系数法对第一主图像和第一辅图像进行配准，及对第二主图像和第二辅图像进行配准，再生成第一干涉图和第二干涉图；
具体的，相关系数法采用基于窗口的配准方法，分别在主图像和辅图像上分别建立匹配窗和搜索窗进行配准。在主图像和辅图像中选取一些控制点，建立的匹配窗和搜索窗，以控制点为中心。滑动匹配窗求得匹配窗和搜索窗的窗口的复相关函数，进而得到偏移量；所述复相关函数可以用表达式(1)表示；
其中，ρc表示复相关函数，s1和s2表示主图像和辅图像，(u,v)表示滑动的匹配窗窗口中心点坐标，(m,n)表示像素中心点坐标，即控制点坐标，M、N是图像的尺寸大小，*表示取复共轭。|□|表示取幅度值，|□|中的“□”可以表示任意表达式。所述匹配窗的大小可以根据处理速度和精度确定，如20x20等；
复相关函数ρc的最大值处即为最佳匹配点，可以对复相关函数的最大值处进行插值，如采用0.5个像素为单元插值，获得更加精确的偏移量，根据控制点对应的配准偏移量获得其他像素处的配准偏移量，可以利用二阶多项式(2)来求取；其中，(x,y)表示主图像中的像素坐标，(Δx,Δy)表示辅图像相对于主图像的偏移量，a和b表示拟合系数，可以根据公式1得到的偏移量拟合得到，其中，a0和b0表示常驻系数，下标表示维数和阶数；a10表示1维0阶系数，a21表示2维1阶系数，以此类推，在此不再赘述。
获取配准后的第一主图像和第一辅图像，以及第二主图像和第二辅图像，可以将配准第一主图像和第一辅图像进行共轭相乘，得到第一轨道对应的第一干涉图；将配准第二主图像和第二辅图像进行共轭相乘，得到第二轨道对应的第二干涉图。
更进一步的，在对所述第一干涉图和第二干涉图进行滤波处理之前，所述方法还包括：对所述第一干涉图和第二干涉图，分别进行去除平地效应处理；
这里，可以采用轨道法去除平地效应；以TSX/TDX为例，TSX/TDX具有精确的轨道参数等信息，可以根据轨道参数和场景地面入射角，计算平地相位，通过和干涉图共轭相乘去除平地效应。
步骤130：根据预设滤波规则，分别对所述第一干涉图和所述第二干涉图进行滤波处理，并确定所述第一干涉图和第二干涉图上各像素对应的相干系数；
由于相干成像系统的固有特性及星载SAR的系统噪声等，第一干涉图和第二干涉图的存在斑点噪声、系统噪声等不同噪声；因此，需要对干涉图进行滤波处理；可以采用均值滤波、Goldstein滤波等方式对所述第一干涉图和第二干涉图进行滤波；
所述第一干涉图上各像素对应的相干系数是指反应所述第一干涉图上各像素对应的配准后的第一主图像和第一辅图像各像素之间的相干性关系的相干系数；所述第二干涉图上各像素对应的相干系数是指反应所述第二干涉图上各像素对应的配准后的第二主图像和第二辅图像各像素之间的相干性关系的相干系数；所述第一干涉图上各像素对应的相干系数可以通过配准后的第一主图像和第一辅图像对应位置的相邻像素采用多视处理等方法得到；所述第二干涉图上各像素对应的相干系数可以通过配准后的第二主图像和第二辅图像对应位置的相邻像素采用多视处理等方法得到。
进一步的，采用非局部滤波法，分别对所述第一干涉图和第二干涉图中的各像素进行迭代估计；根据所述第一干涉图各像素在迭代估计中的相似性权重，采用配准处理后的所述第一主图像和第一辅图像的古德曼Goodman模型，确定所述第一干涉图各像素对应的相干系数；根据所述第二干涉图各像素在迭代估计中的相似性权重，采用配准处理后的所述第二主图像和第二辅图像的Goodman模型，确定所述第二干涉图各像素对应的相干系数；
具体的，两幅配准的SLC复图像每一点的像素模型可以用三个参数来描述：反射系数、真实的干涉相位、相干系数，假设z和z′是两幅配准图像上对应的像素点，根据Goodman模型，z与z′服从表达式(3)；其中，∑表示2×2的相干系数矩阵，可以用表达式(4)表示；其中，R与R′表示反射系数，β表示真实相位，D表示相干系数，E表示取期望值；
可以采用非局部滤波对三个参数R、β、D联合估计，假设两幅图像的散射特性相同，即R＝R′，这种假设在相干性较高的区域是成立的，同时令A＝|z|、A′＝|z′|；φ＝arg(zz′*)为去平地后的相位，则有表达式(5)；
NL-InSAR使用WMLE来估计权重，在分块s中，参数估计可以用表达式(6)表示；
令R＝R′，将公式扩展到权重最大似然估计中，令Θs＝(Rs,βs,Ds)以及
Figure GDA0002221819780000301
估计的R、β、D可以用表达式(7)表示；其中，
Figure GDA0002221819780000302
Figure GDA0002221819780000303
w(s,t)是基于即采用迭代估计的方法实现干涉图的滤波中分块搜索计算的相似性权重；采用迭代估计的方法实现干涉图的滤波的计算方式可以如图3所示；假设图像I上的待估计像素点s可以用表达式(8)表示；判断待估计像素点s，和相似像素点t是否相似，是以s和t为中心选取一定大小的窗口，判断以s，t为中心的这个窗口是否相似；将以像素点s为中心像素窗口作为目标窗31，以像素点t为中心的像素窗口作为相似窗32，在选定搜索窗33内，以预设的路线34进行搜索，确定各相似窗32与目标窗31的相似性，像素点s的估计值可以用表达式(9)表示；其中，vs表示待估计的像素点，vt表示相似窗中心像素；权重w(s,t)的选择是基于s、t像素周围的像素窗口的相似性，满足0≤w(s,t)≤1，并且
Figure GDA0002221819780000304
所述目标窗和相似窗的大小可以根据处理速度和精度确定，如3X3等；
获得了各像素德点估计值，即完成了滤波，这里，对当前像素进行迭代估计就是一个滤波处理过程。
使用NL-InSAR在滤除相位噪声，降噪性能效果好，利于后续的相位解缠。
进一步的，根据述第一干涉图各像素的相干系数，采用相干系数确定规则，确定第一干涉图中各区域的相干系数；所述第一干涉图中的各区域包含一个以上的像素；根据述第二干涉图各像素的相干系数，采用相干系数确定规则，确定第二干涉图中各区域的相干系数；所述第二干涉图中的各区域包含一个以上的像素；
为了提高后续的计算速度，可以针对干涉图中不同区域的相干系数，所述区域可以是干涉图中的一个像素，也可以预设区域大小；可以根据干涉图精度需求来设置区域的大小；采用相干系数确定规则，可以根据区域与所述区域包含像素关系确定，如所述区域为一个像素时，改区域相干系数为该像素的相干系数，所述区域大于一个像素时，所述区域的相干系数可以将区域中各像素相干系数做算数平均等运算确定。
步骤140：根据预设转化规则，将所述滤波后的第一干涉图和第二干涉图，分别转化为第一单基线DEM将和第二单基线DEM；
这里，可以根据干涉图和星载SAR的姿态数据等生成DEM，根据所述第一干涉图和第二干涉图可以分别生成一个与轨道对应的单基线DEM；所述预设转化规则可以根据对第一星载SAR和第二星载SAR的轨道位置、DEM的精度需求、运算时间需求等设置；可以采用常规方式生成DEM。
进一步的，可以将预设DEM仿真干涉相位与第一干涉图的干涉相位进行差分，得到第一差分相位；对所述第一干涉图中相干系数低于第一预设相干阈值的区域的所述第一差分相位进行掩膜，并将掩膜后的第一差分相位进行解缠；在解缠的第一差分相位中，将掩膜区域采用预设插值规则进行插值，并将插值后的解缠第一差分相位加上预设DEM仿真干涉相位，得到第一干涉图的解缠相位；采用InSAR测高模型，根据所述第一干涉图的解缠相位，确定所述第一轨道对应的所述第一单基线DEM；将预设DEM仿真干涉相位与第二干涉图的干涉相位进行差分，得到第二差分相位；对所述第二干涉图中相干系数低于第一预设相干阈值的区域的所述第二差分相位进行掩膜，并将掩膜后的第二差分相位进行解缠；在解缠的第二差分相位中，将掩膜区域采用预设插值规则进行插值，并将插值后的解缠第二差分相位加上预设DEM仿真干涉相位，得到第二干涉图的解缠相位；采用所述InSAR测高模型，根据所述第二干涉图的解缠相位，确定所述第二轨道对应的所述第二单基线DEM；
这里，所述预设DEM可以是现有的成像目标的DEM，根据第一星载SAR和第二星载SAR的SAR系统参数，以及第一星载SAR和第二星载SAR的轨道参数，通过仿真的到所述预设DEM的干涉相位；采用预设DEM的仿真干涉相位与干涉图进行差分干涉处理得到缠绕的差分相位；为了保证相位解缠的可靠性，可以将干涉图中低相干区域掩膜掉，再解缠差分相位；这里可以预设一个第一预设相干阈值，将干涉图各区域中相干系数低于第一预设相干阈值的区域进行掩膜，对掩膜区域进行插值；将掩膜处理后的解缠差分相位加上预设DEM的仿真干涉相位，得到干涉图的解缠相位，这里，可以对第一干涉图和第二干涉图均进行上述处理，分别得到第一干涉图和第二干涉图的解缠相位；
得到解缠相位后，可以采用InSAR几何关系模型反演高度，InSAR数字高程重建的空间几何关系如图4所示，其中：A1和A2分别分别为第一星载SAR和第二星载SAR；H表示平台高度，A1和A2分别表示两个天线的相位中心，α为基线角，B为干涉基线的长度，P为地面上的某一个目标点，r1和r2＝r1+Δr分别为零多普勒平面内两个天线的相位中心到目标点P的距离，θ为主天线相位中心到目标的视角，地形高程用h表示，yP为P点的地距，y为地距向坐标。根据这个空间几何关系可以推导出地形的高程h可以用表达式(10)表示；
通过第一干涉图和第二干涉图的解缠相位分别得到第一干涉图和第二干涉图对应的高程数据后，根据成像目标的地理一起进行编码，得到第一轨道和第二轨道分别对应的第一单基线DEM和第二单基线DEM。
这里，使用外部预设DEM进行差分处理，可以降低相位解缠的难度，保证相位解缠的可靠性。
步骤150：根据所述第一干涉图和第二干涉图上各像素对应的相干系数，采用预设融合规则，融合所述第一单基线DEM和第二单基线DEM，得到成像目标的融合DEM；
所述第一单基线DEM和第二单基线DEM中有较多的数据空洞；无法反应真是的地形特征；这里，可以对第一单基线DEM和第二单基线DEM进行融合；所述预设融合规则可以根据基线DEM中不同区域进行设置，如对于成像遮挡区域等进行互相补充，对正常成像区域进行加权融合等；
进一步的，分别将所述第一单基线DEM和第二单基线DEM中，成像异常区域和对应相干系数低于第二预设相干阈值区域进行掩膜；将第一单基线DEM非掩膜区域和第二单基线DEM非掩膜区域进行加权融合，并将第一单基线DEM正常成像区域补充对应位置的所述第二单基线DEM掩膜区域，将第二单基线DEM正常成像区域补充对应位置的所述第一单基线DEM掩膜区域得到成像目标的融合DEM；将所述第一单基线DEM非掩膜区域对应相干系数占第一单基线DEM非掩膜区域和第二单基线DEM非掩膜区域各自对应相干系数之和的比重，确定为所述第一单基线DEM非掩膜区域进行加权融合的权值；将所述第二单基线DEM非掩膜区域对应相干系数占第一单基线DEM非掩膜区域和第二单基线DEM非掩膜区域各自对应相干系数之和的比重，确定为所述第二单基线DEM非掩膜区域进行加权融合的权值；
具体的，所述成像异常区域可以是单基线DEM中的透视收缩、叠掩、阴影区域；所述对应位置，是指相同地理坐标的位置。可以将第一单基线DEM和第二单基线DEM的透视收缩、叠掩、阴影进行掩膜，同时，采用NL滤波得到的无偏估计的相干系数，将相干系数低于第二预设相干阈值，即相干性较低的区域也进行掩膜；然后利用第一轨道和第二轨道相反照射方向的特点，将空值区域互相补充；将非空值区域加权融合，得到融合后的DEM。所述加权融合中的加权值可以根据第一单基线DEM和第二单基线DEM各像素或区域对应的相干系数确定，如第一单基线DEM各像素或区域对应的相干系数为DA，第二单基线DEM各像素或区域对应的相干系数为DB，那么，所述加权融合中的第一单基线DEM加权值则为
Figure GDA0002221819780000331
第二单基线DEM加权值则为
Figure GDA0002221819780000332
其中，为了简化运算，所述第一预设相干阈值和第二预设相干阈值可以取同一个值；
如此，利用第一轨道和第二轨道不同轨道方向上获得的影像对叠掩等“空白”区域进行互补，如利用降轨DEM上背坡面的区域填补升轨DEM上迎坡面的空白区域，重建高落差物体如建筑物密集区域，细节特征丰富的高分辨率高精度城区DEM。并且由于本发明实施例直接使用一次多视处理操作，避免了多次重复迭代处理，处理起来简单快速，同时避免多次迭代时的误差传递问题。
本发明实施例提供的基于双基星载SAR的高精度DEM反演装置，包括处理器、存储器及存储在存储器上并能够有所述处理器运行的可执行程序，所述处理器运行所述可执行程序时执行基于双基星载SAR的高精度DEM反演方法的步骤；如图1所示，所述基于双基星载SAR的高精度DEM反演方法的步骤包括：
步骤110：获取第一星载SAR和第二星载SAR在第一轨道对成像目标分别进行成像得到的第一主图像和第一辅图像，并获取第一星载SAR和第二星载SAR在第二轨道对所述成像目标分别进行成像得到的第二主图像和第二辅图像；
这里，所述第一星载SAR和第二星载SAR可以是双星分布式系统中的两颗卫星的星载SAR，如TerraSAR-X卫星的星载SAR和TanDEM-X卫星的星载SAR；第一卫星和第二卫星可以是双星分布式系统中的主星和辅星；主星和辅星组成的双星系统，以紧密编队方式飞行对地面进行成像，如TerraSAR-X卫星和TanDEM-X卫星之间距离为200～300；所述成像目标可以是城市等；
星载SAR的成像为SLC；所述第一主图像和第二主图像为主星的SAR分别在第一轨道和第二轨道成像得到的SLC；所述第一辅图像和第二辅图像为辅星的SAR分别在第一轨道和第二轨道成像得到的SLC；
如图2所示的DEM数据处理框架示意图中，第一轨道位置相对较高，第二轨道位置相对较低，在第一轨道的成像称为升轨SLC，在第二轨道的成像称为降轨道SLC。主星和辅星分别在第一轨道对目标进行成像得到升轨SLC1和升轨SLC2，分别在第二轨道对所述目标进行成像得到降轨SLC1和降轨SLC2；获取所述SLC后，可以将所有图像传送到主星进行处理，也可以传送到辅星或地面站进行处理。
进一步的，所述第一星载SAR和第二星载SAR在第一轨道对目标成像时的雷达照射方向，与在第二轨道对目标成像时的雷达照射方向呈预设夹角；
所述第一轨道和第二轨道是指在不同时间，双星系统运行的不同轨道高度和位置；由于地球自转和卫星绕地球飞行的原因，双星系统在不同时间点飞越成像目标时轨道高度不同，卫星雷达照射成像目标的方向也不同；如在第一轨道对成像目标成像时，第一星载SAR和第二星载SAR对目标成像时的雷达照射方向可以是由东向西；在第二轨道对成像目标成像时，第一星载SAR和第二星载SAR对目标成像时的雷达照射方向可以是由西向东。如此，可以在两个方位对成像目标进行成像，避免因为轨道和成像目标的夹角引起的仅能在山峰、高楼等迎坡面成像，增加了背坡面成像的几率。
步骤120：采用预设图像处理规则，处理所述第一主图像和第一辅图像，生成所述第一轨道对应的第一干涉图，处理所述第二主图像和第二辅图像，生成所述第二轨道对应的第二干涉图；
这里，所述预设图像处理规则可以根据第一主图像和第一辅图像及第二主图像和第二辅图像的成像条件，如卫星与成像目标位置关系等，以及成像大小等进行设置；可以对主、辅图像进行图像配准后，再生成干涉图。可以采用最大干涉频谱法和相位差影像平均波动函数法等方法进行配准；可以将配准后的第一主图像和第一辅图像共轭相乘得到第一轨道对应的第一干涉图，将配准后的第二主图像和第二辅图像共轭相乘得到第二轨道对应的第二干涉图。
进一步的，对所述第一主图像和第一辅图像进行图像配准，将配准处理后的第一主图像和第一辅图像共轭相乘得到第一干涉图；对所述第二主图像和第二辅图像进行图像配准，将配准处理后的第二主图像和第二辅图像共轭相乘得到第二干涉图；这里，可以采用相关系数法对第一主图像和第一辅图像进行配准，及对第二主图像和第二辅图像进行配准，再生成第一干涉图和第二干涉图；
具体的，相关系数法采用基于窗口的配准方法，分别在主图像和辅图像上分别建立匹配窗和搜索窗进行配准。在主图像和辅图像中选取一些控制点，建立的匹配窗和搜索窗，以控制点为中心。滑动匹配窗求得匹配窗和搜索窗的窗口的复相关函数，进而得到偏移量；所述复相关函数可以用表达式(1)表示；
其中，ρc表示复相关函数，s1和s2表示主图像和辅图像，(u,v)表示滑动的匹配窗窗口中心点坐标，(m,n)表示像素中心点坐标，即控制点坐标，M、N是图像的尺寸大小，*表示取复共轭。|□|表示取幅度值，|□|中的“□”可以表示任意表达式。所述匹配窗的大小可以根据处理速度和精度确定，如20x20等；
复相关函数ρc的最大值处即为最佳匹配点，可以对复相关函数的最大值处进行插值，如采用0.5个像素为单元插值，获得更加精确的偏移量，根据控制点对应的配准偏移量获得其他像素处的配准偏移量，可以利用二阶多项式(2)来求取；其中，(x,y)表示主图像中的像素坐标，(Δx,Δy)表示辅图像相对于主图像的偏移量，a和b表示拟合系数，可以根据公式1得到的偏移量拟合得到，其中，a0和b0表示常驻系数，下标表示维数和阶数；a10表示1维0阶系数，a21表示2维1阶系数，以此类推，在此不再赘述。
获取配准后的第一主图像和第一辅图像，以及第二主图像和第二辅图像，可以将配准第一主图像和第一辅图像进行共轭相乘，得到第一轨道对应的第一干涉图；将配准第二主图像和第二辅图像进行共轭相乘，得到第二轨道对应的第二干涉图。
更进一步的，在对所述第一干涉图和第二干涉图进行滤波处理之前，所述方法还包括：对所述第一干涉图和第二干涉图，分别进行去除平地效应处理；
这里，可以采用轨道法去除平地效应；以TSX/TDX为例，TSX/TDX具有精确的轨道参数等信息，可以根据轨道参数和场景地面入射角，计算平地相位，通过和干涉图共轭相乘去除平地效应。
步骤130：根据预设滤波规则，分别对所述第一干涉图和所述第二干涉图进行滤波处理，并确定所述第一干涉图和第二干涉图上各像素对应的相干系数；
由于相干成像系统的固有特性及星载SAR的系统噪声等，第一干涉图和第二干涉图的存在斑点噪声、系统噪声等不同噪声；因此，需要对干涉图进行滤波处理；可以采用均值滤波、Goldstein滤波等方式对所述第一干涉图和第二干涉图进行滤波；
所述第一干涉图上各像素对应的相干系数是指反应所述第一干涉图上各像素对应的配准后的第一主图像和第一辅图像各像素之间的相干性关系的相干系数；所述第二干涉图上各像素对应的相干系数是指反应所述第二干涉图上各像素对应的配准后的第二主图像和第二辅图像各像素之间的相干性关系的相干系数；所述第一干涉图上各像素对应的相干系数可以通过配准后的第一主图像和第一辅图像对应位置的相邻像素采用多视处理等方法得到；所述第二干涉图上各像素对应的相干系数可以通过配准后的第二主图像和第二辅图像对应位置的相邻像素采用多视处理等方法得到。
进一步的，采用非局部滤波法，分别对所述第一干涉图和第二干涉图中的各像素进行迭代估计；根据所述第一干涉图各像素在迭代估计中的相似性权重，采用配准处理后的所述第一主图像和第一辅图像的古德曼Goodman模型，确定所述第一干涉图各像素对应的相干系数；根据所述第二干涉图各像素在迭代估计中的相似性权重，采用配准处理后的所述第二主图像和第二辅图像的Goodman模型，确定所述第二干涉图各像素对应的相干系数；
具体的，两幅配准的SLC复图像每一点的像素模型可以用三个参数来描述：反射系数、真实的干涉相位、相干系数，假设z和z′是两幅配准图像上对应的像素点，根据Goodman模型，z与z′服从表达式(3)；其中，∑表示2×2的相干系数矩阵，可以用表达式(4)表示；其中，R与R′表示反射系数，β表示真实相位，D表示相干系数，E表示取期望值；
可以采用非局部滤波对三个参数R、β、D联合估计，假设两幅图像的散射特性相同，即R＝R′，这种假设在相干性较高的区域是成立的，同时令A＝|z|、A′＝|z′|；φ＝arg(zz′*)为去平地后的相位，则有表达式(5)；
NL-InSAR使用WMLE来估计权重，在分块s中，参数估计可以用表达式(6)表示；
令R＝R′，将公式扩展到权重最大似然估计中，令Θs＝(Rs,βs,Ds)以及
Figure GDA0002221819780000381
估计的R、β、D可以用表达式(7)表示；其中，
Figure GDA0002221819780000382
w(s,t)是基于即采用迭代估计的方法实现干涉图的滤波中分块搜索计算的相似性权重；采用迭代估计的方法实现干涉图的滤波的计算方式可以如图3所示；假设图像I上的待估计像素点s可以用表达式(8)表示；判断待估计像素点s，和相似像素点t是否相似，是以s和t为中心选取一定大小的窗口，判断以s，t为中心的这个窗口是否相似；将以像素点s为中心像素窗口作为目标窗31，以像素点t为中心的像素窗口作为相似窗32，在选定搜索窗33内，以预设的路线34进行搜索，确定各相似窗32与目标窗31的相似性，像素点s的估计值可以用表达式(9)表示；其中，vs表示待估计的像素点，vt表示相似窗中心像素；权重w(s,t)的选择是基于s、t像素周围的像素窗口的相似性，满足0≤w(s,t)≤1，并且
Figure GDA0002221819780000384
所述目标窗和相似窗的大小可以根据处理速度和精度确定，如3X3等；
获得了各像素德点估计值，即完成了滤波，这里，对当前像素进行迭代估计就是一个滤波处理过程。
使用NL-InSAR在滤除相位噪声，降噪性能效果好，利于后续的相位解缠。
进一步的，根据述第一干涉图各像素的相干系数，采用相干系数确定规则，确定第一干涉图中各区域的相干系数；所述第一干涉图中的各区域包含一个以上的像素；根据述第二干涉图各像素的相干系数，采用相干系数确定规则，确定第二干涉图中各区域的相干系数；所述第二干涉图中的各区域包含一个以上的像素；
为了提高后续的计算速度，可以针对干涉图中不同区域的相干系数，所述区域可以是干涉图中的一个像素，也可以预设区域大小；可以根据干涉图精度需求来设置区域的大小；采用相干系数确定规则，可以根据区域与所述区域包含像素关系确定，如所述区域为一个像素时，改区域相干系数为该像素的相干系数，所述区域大于一个像素时，所述区域的相干系数可以将区域中各像素相干系数做算数平均等运算确定。
步骤140：根据预设转化规则，将所述滤波后的第一干涉图和第二干涉图，分别转化为第一单基线DEM将和第二单基线DEM；
这里，可以根据干涉图和星载SAR的姿态数据等生成DEM，根据所述第一干涉图和第二干涉图可以分别生成一个与轨道对应的单基线DEM；所述预设转化规则可以根据对第一星载SAR和第二星载SAR的轨道位置、DEM的精度需求、运算时间需求等设置；可以采用常规方式生成DEM。
进一步的，可以将预设DEM仿真干涉相位与第一干涉图的干涉相位进行差分，得到第一差分相位；对所述第一干涉图中相干系数低于第一预设相干阈值的区域的所述第一差分相位进行掩膜，并将掩膜后的第一差分相位进行解缠；在解缠的第一差分相位中，将掩膜区域采用预设插值规则进行插值，并将插值后的解缠第一差分相位加上预设DEM仿真干涉相位，得到第一干涉图的解缠相位；采用InSAR测高模型，根据所述第一干涉图的解缠相位，确定所述第一轨道对应的所述第一单基线DEM；将预设DEM仿真干涉相位与第二干涉图的干涉相位进行差分，得到第二差分相位；对所述第二干涉图中相干系数低于第一预设相干阈值的区域的所述第二差分相位进行掩膜，并将掩膜后的第二差分相位进行解缠；在解缠的第二差分相位中，将掩膜区域采用预设插值规则进行插值，并将插值后的解缠第二差分相位加上预设DEM仿真干涉相位，得到第二干涉图的解缠相位；采用所述InSAR测高模型，根据所述第二干涉图的解缠相位，确定所述第二轨道对应的所述第二单基线DEM；
这里，所述预设DEM可以是现有的成像目标的DEM，根据第一星载SAR和第二星载SAR的SAR系统参数，以及第一星载SAR和第二星载SAR的轨道参数，通过仿真的到所述预设DEM的干涉相位；采用预设DEM的仿真干涉相位与干涉图进行差分干涉处理得到缠绕的差分相位；为了保证相位解缠的可靠性，可以将干涉图中低相干区域掩膜掉，再解缠差分相位；这里可以预设一个第一预设相干阈值，将干涉图各区域中相干系数低于第一预设相干阈值的区域进行掩膜，对掩膜区域进行插值；将掩膜处理后的解缠差分相位加上预设DEM的仿真干涉相位，得到干涉图的解缠相位，这里，可以对第一干涉图和第二干涉图均进行上述处理，分别得到第一干涉图和第二干涉图的解缠相位；
得到解缠相位后，可以采用InSAR几何关系模型反演高度，InSAR数字高程重建的空间几何关系如图4所示，其中：A1和A2分别分别为第一星载SAR和第二星载SAR；H表示平台高度，A1和A2分别表示两个天线的相位中心，α为基线角，B为干涉基线的长度，P为地面上的某一个目标点，r1和r2＝r1+Δr分别为零多普勒平面内两个天线的相位中心到目标点P的距离，θ为主天线相位中心到目标的视角，地形高程用h表示，yP为P点的地距，y为地距向坐标。根据这个空间几何关系可以推导出地形的高程h可以用表达式(10)表示；
通过第一干涉图和第二干涉图的解缠相位分别得到第一干涉图和第二干涉图对应的高程数据后，根据成像目标的地理一起进行编码，得到第一轨道和第二轨道分别对应的第一单基线DEM和第二单基线DEM。
这里，使用外部预设DEM进行差分处理，可以降低相位解缠的难度，保证相位解缠的可靠性。
步骤150：根据所述第一干涉图和第二干涉图上各像素对应的相干系数，采用预设融合规则，融合所述第一单基线DEM和第二单基线DEM，得到成像目标的融合DEM；
所述第一单基线DEM和第二单基线DEM中有较多的数据空洞；无法反应真是的地形特征；这里，可以对第一单基线DEM和第二单基线DEM进行融合；所述预设融合规则可以根据基线DEM中不同区域进行设置，如对于成像遮挡区域等进行互相补充，对正常成像区域进行加权融合等；
进一步的，分别将所述第一单基线DEM和第二单基线DEM中，成像异常区域和对应相干系数低于第二预设相干阈值区域进行掩膜；将第一单基线DEM非掩膜区域和第二单基线DEM非掩膜区域进行加权融合，并将第一单基线DEM正常成像区域补充对应位置的所述第二单基线DEM掩膜区域，将第二单基线DEM正常成像区域补充对应位置的所述第一单基线DEM掩膜区域得到成像目标的融合DEM；将所述第一单基线DEM非掩膜区域对应相干系数占第一单基线DEM非掩膜区域和第二单基线DEM非掩膜区域各自对应相干系数之和的比重，确定为所述第一单基线DEM非掩膜区域进行加权融合的权值；将所述第二单基线DEM非掩膜区域对应相干系数占第一单基线DEM非掩膜区域和第二单基线DEM非掩膜区域各自对应相干系数之和的比重，确定为所述第二单基线DEM非掩膜区域进行加权融合的权值；
具体的，所述成像异常区域可以是单基线DEM中的透视收缩、叠掩、阴影区域；所述对应位置，是指相同地理坐标的位置。可以将第一单基线DEM和第二单基线DEM的透视收缩、叠掩、阴影进行掩膜，同时，采用NL滤波得到的无偏估计的相干系数，将相干系数低于第二预设相干阈值，即相干性较低的区域也进行掩膜；然后利用第一轨道和第二轨道相反照射方向的特点，将空值区域互相补充；将非空值区域加权融合，得到融合后的DEM。所述加权融合中的加权值可以根据第一单基线DEM和第二单基线DEM各像素或区域对应的相干系数确定，如第一单基线DEM各像素或区域对应的相干系数为DA，第二单基线DEM各像素或区域对应的相干系数为DB，那么，所述加权融合中的第一单基线DEM加权值则为
Figure GDA0002221819780000411
第二单基线DEM加权值则为
Figure GDA0002221819780000412
其中，为了简化运算，所述第一预设相干阈值和第二预设相干阈值可以取同一个值；
如此，利用第一轨道和第二轨道不同轨道方向上获得的影像对叠掩等“空白”区域进行互补，如利用降轨DEM上背坡面的区域填补升轨DEM上迎坡面的空白区域，重建高落差物体如建筑物密集区域，细节特征丰富的高分辨率高精度城区DEM。并且由于本发明实施例直接使用一次多视处理操作，避免了多次重复迭代处理，处理起来简单快速，同时避免多次迭代时的误差传递问题。
以上所述，仅为本发明的最佳实施例而已，并非用于限定本发明的保护范围，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。

Claims (16)

1.一种基于双基星载SAR的高精度DEM反演方法，其特征在于，所述方法包括：
获取第一星载合成孔径雷达SAR和第二星载SAR在第一轨道对成像目标分别进行成像得到的第一主图像和第一辅图像，并获取第一星载SAR和第二星载SAR在第二轨道对所述成像目标分别进行成像得到的第二主图像和第二辅图像；
采用预设图像处理规则，处理所述第一主图像和第一辅图像，生成所述第一轨道对应的第一干涉图，处理所述第二主图像和第二辅图像，生成所述第二轨道对应的第二干涉图；
根据预设滤波规则，分别对所述第一干涉图和所述第二干涉图进行滤波处理，并确定所述第一干涉图和第二干涉图上各像素对应的相干系数；
根据预设转化规则，将所述滤波后的第一干涉图和第二干涉图，分别转化为第一单基线数字高程模型DEM和第二单基线DEM；
根据所述第一干涉图和第二干涉图上各像素对应的相干系数，采用预设融合规则，融合所述第一单基线DEM和第二单基线DEM，得到成像目标的融合DEM；其中，所述根据所述第一干涉图和第二干涉图上各像素对应的相干系数，采用预设融合规则，融合所述第一单基线DEM和第二单基线DEM，得到成像目标的融合DEM；包括：分别将所述第一单基线DEM和第二单基线DEM中，成像异常区域和对应相干系数低于第一预设相干阈值区域进行掩膜；将第一单基线DEM非掩膜区域和第二单基线DEM非掩膜区域进行加权融合，并将第一单基线DEM正常成像区域补充对应位置的所述第二单基线DEM掩膜区域，将第二单基线DEM正常成像区域补充对应位置的所述第一单基线DEM掩膜区域得到成像目标的融合DEM；将所述第一单基线DEM非掩膜区域对应相干系数占第一单基线DEM非掩膜区域和第二单基线DEM非掩膜区域各自对应相干系数之和的比重，确定为所述第一单基线DEM非掩膜区域进行加权融合的权值；将所述第二单基线DEM非掩膜区域对应相干系数占第一单基线DEM非掩膜区域和第二单基线DEM非掩膜区域各自对应相干系数之和的比重，确定为所述第二单基线DEM非掩膜区域进行加权融合的权值。
2.根据权利要求1所述的方法，其特征在于，所述采用预设图像处理规则，处理所述第一主图像和第一辅图像，生成所述第一轨道对应的第一干涉图，处理所述第二主图像和第二辅图像，生成所述第二轨道对应的第二干涉图，包括：
对所述第一主图像和第一辅图像进行图像配准，将配准处理后的第一主图像和第一辅图像共轭相乘得到第一干涉图；
对所述第二主图像和第二辅图像进行图像配准，将配准处理后的第二主图像和第二辅图像共轭相乘得到第二干涉图。
3.根据权利要求2所述的方法，其特征在于，所述根据预设滤波规则，分别对所述第一干涉图和第二干涉图进行滤波处理，并确定所述第一干涉图和第二干涉图上各像素对应的相干系数，包括：
采用非局部滤波法，分别对所述第一干涉图和第二干涉图中的各像素进行迭代估计；
根据所述第一干涉图各像素在迭代估计中的相似性权重，采用配准处理后的所述第一主图像和第一辅图像的古德曼Goodman模型，得到所述第一干涉图各像素对应的相干系数；
根据所述第二干涉图各像素在迭代估计中的相似性权重，采用配准处理后的所述第二主图像和第二辅图像的Goodman模型，得到所述第二干涉图各像素对应的相干系数。
4.根据权利要求3所述的方法，其特征在于，所述根据预设转化规则，将所述滤波后的第一干涉图和第二干涉图，分别转化为第一单基线DEM和第二单基线DEM；包括：
将预设DEM仿真干涉相位与第一干涉图的干涉相位进行差分，得到第一差分相位；
对所述第一干涉图中相干系数低于第二预设相干阈值的区域的所述第一差分相位进行掩膜，并将掩膜后的第一差分相位进行解缠；
在解缠的第一差分相位中，将掩膜区域采用预设插值规则进行插值，并将插值后的解缠第一差分相位加上预设DEM仿真干涉相位，得到第一干涉图的解缠相位；
采用干涉合成孔径雷达InSAR测高模型，根据所述第一干涉图的解缠相位，确定所述第一轨道对应的所述第一单基线DEM；
将预设DEM仿真干涉相位与第二干涉图的干涉相位进行差分，得到第二差分相位；
对所述第二干涉图中相干系数低于第二预设相干阈值的区域的所述第二差分相位进行掩膜，并将掩膜后的第二差分相位进行解缠；
在解缠的第二差分相位中，将掩膜区域采用预设插值规则进行插值，并将插值后的解缠第二差分相位加上预设DEM仿真干涉相位，得到第二干涉图的解缠相位；
采用所述InSAR测高模型，根据所述第二干涉图的解缠相位，确定所述第二轨道对应的所述第二单基线DEM。
5.根据权利要求1或4所述的方法，其特征在于，所述方法还包括：
根据所述第一干涉图各像素的相干系数，采用相干系数确定规则，确定第一干涉图中各区域的相干系数；
所述第一干涉图中的各区域包含一个以上的像素；
根据所 述第二干涉图各像素的相干系数，采用相干系数确定规则，确定第二干涉图中各区域的相干系数；
所述第二干涉图中的各区域包含一个以上的像素。
6.根据权利要求1至4任一项所述的方法，其特征在于，在对所述第一干涉图和第二干涉图进行滤波处理之前，所述方法还包括：
对所述第一干涉图和第二干涉图，分别进行去除平地效应处理。
7.根据权利要求1至4任一项所述的方法，其特征在于，所述第一星载SAR和第二星载SAR在第一轨道对目标成像时的雷达照射方向，与在第二轨道对目标成像时的雷达照射方向呈预设夹角。
8.一种基于双基星载SAR的高精度DEM反演装置，其特征在于，所述装置包括：获取模块、图像处理模块、滤波模块、转化模块和融合模块；其中，
所述获取模块，用于获取第一星载SAR和第二星载SAR在第一轨道对成像目标分别进行成像得到的第一主图像和第一辅图像，并获取第一星载SAR和第二星载SAR在第二轨道对所述成像目标分别进行成像得到的第二主图像和第二辅图像；
所述图像处理模块，用于采用预设图像处理规则，处理所述第一主图像和第一辅图像，生成所述第一轨道对应的第一干涉图，处理所述第二主图像和第二辅图像，生成所述第二轨道对应的第二干涉图；
所述滤波模块，用于根据预设滤波规则，分别对所述第一干涉图和第二干涉图进行滤波处理，并确定所述第一干涉图和所述第二干涉图上各像素对应的相干系数；
所述转化模块，用于根据预设转化规则，将所述滤波后的第一干涉图和第二干涉图，分别转化为第一单基线DEM和第二单基线DEM；
所述融合模块，用于根据所述第一干涉图和第二干涉图上各像素对应的相干系数，采用预设融合规则，融合所述第一单基线DEM和第二单基线DEM，得到成像目标的融合DEM；其中，所述融合模块，具体用于：分别将所述第一单基线DEM和第二单基线DEM中，成像异常区域和对应相干系数低于第一预设相干阈值区域进行掩膜；将第一单基线DEM非掩膜区域和第二单基线DEM非掩膜区域进行加权融合，并将第一单基线DEM正常成像区域补充对应位置的所述第二单基线DEM掩膜区域，将第二单基线DEM正常成像区域补充对应位置的所述第一单基线DEM掩膜区域得到成像目标的融合DEM；将所述第一单基线DEM非掩膜区域对应相干系数占第一单基线DEM和第二单基线DEM非掩膜区域各自对应相干系数之和的比重，确定为所述第一单基线DEM非掩膜区域进行加权融合的权值；将所述第二单基线DEM非掩膜区域对应相干系数占第一单基线DEM和第二单基线DEM非掩膜区域各自对应相干系数之和的比重，确定为所述第二单基线DEM非掩膜区域进行加权融合的权值。
9.根据权利要求8所述的装置，其特征在于，所述图像处理模块，具体用于：
对所述第一主图像和第一辅图像进行图像配准，将配准处理后的第一主图像和第一辅图像共轭相乘得到第一干涉图；
对所述第二主图像和第二辅图像进行图像配准，将配准处理后的第二主图像和第二辅图像共轭相乘得到第二干涉图。
10.根据权利要求9所述的装置，其特征在于，所述滤波模块，具体用于：
采用非局部滤波法，分别对所述第一干涉图和第二干涉图中的各像素进行迭代估计；
根据所述第一干涉图各像素在迭代估计中的相似性权重，采用配准处理后的所述第一主图像和第一辅图像的Goodman模型，得到所述第一干涉图各像素对应的相干系数；
根据所述第二干涉图各像素在迭代估计中的相似性权重，采用配准处理后的所述第二主图像和第二辅图像的Goodman模型，得到所述第二干涉图各像素对应的相干系数。
11.根据权利要求10所述的装置，其特征在于，所述转化模块，具体用于：
将预设DEM仿真干涉相位与第一干涉图的干涉相位进行差分，得到第一差分相位；
对所述第一干涉图中相干系数低于第二预设相干阈值的区域的所述第一差分相位进行掩膜，并将掩膜后的第一差分相位进行解缠；
在解缠的第一差分相位中，将掩膜区域采用预设插值规则进行插值，并将插值后的解缠第一差分相位加上预设DEM仿真干涉相位，得到第一干涉图的解缠相位；
采用InSAR测高模型，根据所述第一干涉图的解缠相位，确定所述第一轨道对应的所述第一单基线DEM；
将预设DEM仿真干涉相位与第二干涉图的干涉相位进行差分，得到第二差分相位；
对所述第二干涉图中相干系数低于第二预设相干阈值的区域的所述第二差分相位进行掩膜，并将掩膜后的第二差分相位进行解缠；
在解缠的第二差分相位中，将掩膜区域采用预设插值规则进行插值，并将插值后的解缠第二差分相位加上预设DEM仿真干涉相位，得到第二干涉图的解缠相位；
采用InSAR测高模型，根据所述第二干涉图的解缠相位，确定所述第二轨道对应的所述第二单基线DEM。
12.根据权利要求8或11所述的装置，其特征在于，所述滤波模块，还用于：
根据所述第一干涉图各像素的相干系数，采用相干系数确定规则，确定第一干涉图中各区域的相干系数；
所述第一干涉图中的各区域包含一个以上的像素；
根据所 述第二干涉图各像素的相干系数，采用相干系数确定规则，确定第二干涉图中各区域的相干系数；
所述第二干涉图中的各区域包含一个以上的像素。
13.根据权利要求8至11任一项所述的装置，其特征在于，所述图像处理模块，还用于：
对所述第一干涉图和第二干涉图，分别进行去除平地效应处理。
14.根据权利要求8至11任一项所述的装置，其特征在于，所述第一星载SAR和第二星载SAR在第一轨道对目标成像时的雷达照射方向，与在第二轨道对目标成像时的雷达照射方向呈预设夹角。
15.一种存储介质，其上存储有 可执行程序，其特征在于，所述可执行程序被处理器执行时实现如权利要求1至7任一项所述基于双基星载SAR的高精度DEM反演方法的步骤。
16.一种基于双基星载SAR的高精度DEM反演装置，包括处理器、存储器及存储在存储器上并能够由 所述处理器运行的可执行程序，其特征在于，所述处理器运行所述可执行程序时执行如权利要求1至7任一项所述基于双基星载SAR的高精度DEM反演方法的步骤。
CN201810525764.0A 2018-05-28 2018-05-28 一种基于双基星载sar的高精度dem反演方法及装置 Active CN109212522B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810525764.0A CN109212522B (zh) 2018-05-28 2018-05-28 一种基于双基星载sar的高精度dem反演方法及装置

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810525764.0A CN109212522B (zh) 2018-05-28 2018-05-28 一种基于双基星载sar的高精度dem反演方法及装置

Publications (2)

Publication Number Publication Date
CN109212522A CN109212522A (zh) 2019-01-15
CN109212522B true CN109212522B (zh) 2020-01-21

Family

ID=64991141

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810525764.0A Active CN109212522B (zh) 2018-05-28 2018-05-28 一种基于双基星载sar的高精度dem反演方法及装置

Country Status (1)

Country Link
CN (1) CN109212522B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN111538006A (zh) * 2020-05-13 2020-08-14 深圳大学 基于动态基线的InSAR数字高程模型构建方法及系统
CN112419198B (zh) * 2020-11-27 2024-02-02 中国矿业大学 一种用于sar干涉图滤波的非局部均值定权方法
CN113065467A (zh) * 2021-04-01 2021-07-02 中科星图空间技术有限公司 一种基于深度学习的卫星图像低相干区域识别方法及装置
CN113466855B (zh) * 2021-07-22 2023-04-25 中国科学院空天信息创新研究院 一种信号重构方法及装置
CN113624122B (zh) * 2021-08-10 2022-09-20 中咨数据有限公司 融合GNSS数据与InSAR技术的桥梁变形监测方法
CN114236544B (zh) * 2022-02-23 2022-05-13 中国科学院空天信息创新研究院 一种基于几何匹配的升降轨星载sar三维成像方法及装置

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103424744B (zh) * 2012-05-16 2015-03-25 中国科学院电子学研究所 干涉sar叠掩区域数字高程模型重建的方法
CN105929398B (zh) * 2016-04-20 2018-11-02 中国电力工程顾问集团中南电力设计院有限公司 结合外控点的InSAR高精度高分辨率DEM获取方法
WO2018027332A1 (es) * 2016-08-08 2018-02-15 Comercial E Industrial Gesecology Limitada Método y sistema para el análisis y generación de alertas tempranas o predictivas de estabilidad de taludes en minas a rajo abierto
CN107102333B (zh) * 2017-06-27 2020-01-10 北京航空航天大学 一种星载InSAR长短基线融合解缠方法

Also Published As

Publication number Publication date
CN109212522A (zh) 2019-01-15

Similar Documents

Publication Publication Date Title
CN109212522B (zh) 一种基于双基星载sar的高精度dem反演方法及装置
Catalão et al. Merging GPS and atmospherically corrected InSAR data to map 3-D terrain displacement velocity
CN111856459B (zh) 一种改进的DEM最大似然约束多基线InSAR相位解缠方法
CN109031301A (zh) 基于PSInSAR技术的山区地形形变提取方法
CN105677942A (zh) 一种重复轨道星载自然场景sar复图像数据快速仿真方法
CN104007439B (zh) 一种干涉圆迹sar高程估计处理方法
CN113960595A (zh) 一种地表形变监测方法及系统
CN109100719A (zh) 基于星载sar影像与光学影像的地形图联合测图方法
Méric et al. A multiwindow approach for radargrammetric improvements
CN109270527B (zh) 圆迹sar子孔径图像序列联合相关dem提取方法
CN103454636A (zh) 基于多像素协方差矩阵的差分干涉相位估计方法
CN107064933A (zh) 基于循环谱估计的sar层析建筑物高度的方法
CN113238228B (zh) 基于水准约束的三维地表形变获取方法、系统及装置
Joughin Estimation of ice-sheet topography and motion using interferometric synthetic aperture radar
Feng et al. A hierarchical network densification approach for reconstruction of historical ice velocity fields in East Antarctica
CN110618409B (zh) 顾及叠掩及阴影的多通道InSAR干涉图仿真方法及系统
CN109946682B (zh) 基于ICESat/GLAS的GF3数据基线估计方法
Liu et al. A comparative study of DEM reconstruction using the single-baseline and multibaseline InSAR techniques
CN116299453A (zh) 星载sar非沿迹模式干涉高程测量方法
CN115546264A (zh) 一种星载InSAR图像精细配准与立体测量方法
CN112835043B (zh) 一种任意方向的二维形变的监测方法
Recla et al. From Relative to Absolute Heights in SAR-based Single-Image Height Prediction
CN115540908A (zh) 基于小波变换的InSAR干涉条纹匹配方法
CN111856464B (zh) 一种基于单控制点信息的车载sar的dem提取方法
Lombardi et al. Accuracy of high resolution CSK interferometric Digital Elevation Models

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