CN109085587A - 机载InSAR困难区域的DEM重建方法 - Google Patents

机载InSAR困难区域的DEM重建方法 Download PDF

Info

Publication number
CN109085587A
CN109085587A CN201810919415.7A CN201810919415A CN109085587A CN 109085587 A CN109085587 A CN 109085587A CN 201810919415 A CN201810919415 A CN 201810919415A CN 109085587 A CN109085587 A CN 109085587A
Authority
CN
China
Prior art keywords
insar
dem
pixel grid
data
measure
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.)
Pending
Application number
CN201810919415.7A
Other languages
English (en)
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 CN201810919415.7A priority Critical patent/CN109085587A/zh
Publication of CN109085587A publication Critical patent/CN109085587A/zh
Pending legal-status Critical Current

Links

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

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)
  • Image Processing (AREA)

Abstract

一种机载InSAR困难区域的DEM重建方法，包括：对各个方向照射的InSAR数据，分别进行初步干涉处理，生成主、辅图像、干涉相位图及相干系数图；将待测绘区域投影到高斯坐标系下，并划分像素网格；将粗DEM数据插值使其与上述像素网格一一对应；将粗DEM值为待测量高程的初值，设置该像素网格高程值的变化范围；设置每一次迭代处理的高程值；分别计算该像素网格投影到各个方向对应的主图像斜距平面中的坐标，读取对应坐标处的干涉相位和相干系数观测值，并计算迭代高程值对应的理想干涉相位；计算迭代高程值的联合似然函数；取使联合似然函数最大时对应的迭代高程值即为该像素网格的最终高程估计值，即重建得到了每个像素网格的DEM，完成机载InSAR困难区域的DEM重建。

Description

机载InSAR困难区域的DEM重建方法
技术领域
本公开涉及干涉合成孔径雷达信号处理领域，尤其涉及一种机载InSAR困难区域的DEM重建方法。
背景技术
干涉合成孔径雷达(Interferometric Synthetic Aperture Radar，InSAR)技术将无线电干涉测量技术与SAR(合成孔径雷达)技术相结合，在获取二维SAR图像的同时，能够利用SAR复数据的相位信息提取地表的三维信息和变化信息，从而将SAR的测量空间从二维拓展到三维，随着InSAR系统的迅速发展，InSAR技术的应用领域也得以不断扩展，其中，地形测绘一直是InSAR技术最直接和最主要的应用之一，与摄影测量、激光雷达等技术相比，InSAR具有全天时、全天候的优点，与雷达立体测量技术相比，InSAR具有更高的测量精度。
然而，由于SAR侧视成像的特点，使得SAR图像中存在几何畸变现象，对于地形起伏较大的城市、山区等区域而言，这一现象尤为严重，影响InSAR地形测绘的几何畸变现象包括阴影和叠掩，阴影区域由于雷达接收不到地表的有用回波信息，相应的干涉相位表现为噪声特性，其中并不包含区域本身的真实地形信息，无法用于高程反演；而叠掩区域则是不同高度区域的回波投影到同一个距离-多普勒单元内形成，其干涉相位反映的是多个散射源矢量叠加的结果，常规的InSAR处理方法无法分辨同一采样单元内的多个散射源的相位信息，从而也无法准确反演其高程。由此可见，在起伏的地形条件下，常规的干涉处理方法无法实现阴影、叠掩区域的高程测量，这些区域是InSAR处理的困难区域，也是限制InSAR地形测绘应用实用化的关键。
目前，对于InSAR困难区域的处理方法主要包括以下两个方面：首先，对于叠掩区域的处理，现有技术主要利用多基线InSAR技术实现多个散射源的检测分辨和高程估计，然而，这类方法要达到较高的估计精度，需要较多的基线数，因而在实际系统中难以实现。其次，对于阴影区域的处理，现有技术主要集中在对其的检测识别上，在检测的基础上，一方面，在干涉处理过程中，对这些区域进行掩膜处理，从而避免相位解缠时误差的传播；另一方面，利用外源的DEM(Digital Elevation Model，数字高程模型)信息对检测出的区域进行补充，但是，可以公开获取的外源DEM数据通常精度较低，无法达到机载InSAR的测绘精度要求，从而影响整体制图的质量。
公开内容
(一)要解决的技术问题
本公开提供了一种机载InSAR困难区域DEM重建方法，其基于多视向数据拼接，以缓解现有技术中达到较高估计精度需较多基线数，及公开获取的外源DEM数据通常精度较低，无法达到机载InSAR的测绘精度要求，从而影响整体制图的质量等技术问题。
(二)技术方案
本公开提供一种机载InSAR困难区域的DEM重建方法，包括：步骤A：对各个方向照射的InSAR数据，分别进行初步的干涉处理，生成主图像和辅图像、干涉相位图以及相干系数图；步骤B：将待测绘DEM的场景区域四角点的地理坐标，利用高斯投影的方法投影到高斯坐标系下，并将高斯投影后高斯坐标系下的矩形区域按照生成DEM产品要求的分辨率划分为等间隔的像素网格；步骤C：选择待测绘DEM的场景区域对应的外部粗DEM数据，根据所选粗DEM数据的格式将粗DEM数据转换到高斯坐标系下，并将所述粗DEM数据插值到生成产品要求的分辨率，使粗DEM数据与步骤B所确定的等间隔的像素网格一一对应；步骤D：对步骤B所划分的每一个像素网格，以该像素网格对应的粗DEM值为待测量高程的初值，设置该像素网格高程值的变化范围；步骤E：对步骤D所处理后的每一个像素网格，设置每一次迭代处理的高程值；步骤F：对于步骤E所述的每次迭代的高程值，根据各个方向多次照射的InSAR数据主图像的成像几何，分别计算该像素网格投影到各个方向对应的主图像斜距平面中的距离向和方位向坐标，并读取步骤A所生成的干涉相位图中该像素网格对应坐标处的干涉相位观测值，以及步骤A所生成的相干系数图中该像素网格对应坐标处的相干系数观测值；步骤G：对于步骤E所述的每次迭代的高程值，根据各个方向多次照射的InSAR数据的成像几何，分别计算所述高程值对应的理想的干涉相位；步骤H：对于步骤E所述的每次迭代的高程值，针对各个方向照射的InSAR数据，分别计算在干涉相位观测值一定的前提下关于所述高程值的似然函数；步骤I：对于步骤E所述的每次迭代的高程值，计算各个方向多次照射的InSAR数据在联合的干涉相位观测值一定的前提下关于迭代高程值的联合似然函数；以及步骤J：迭代次数每次增加1，重复步骤E至步骤I，取使联合似然函数最大时对应的迭代高程值即为该像素网格的最终高程估计值，即重建得到了步骤E所述的每个像素网格的DEM，完成机载InSAR困难区域的DEM重建。
在本公开实施例中，所述步骤A，包括：步骤A1：对各个方向照射的InSAR数据生成的主、辅图像进行复图像配准；步骤A2：对步骤A1配准后的主、辅图像进行共轭相乘，对共轭相乘的结果取相位，生成干涉相位图；以及步骤A3：对步骤A1配准后的主、辅图像计算每一个像素的相干系数，生成相干系数图。
在本公开实施例中，步骤C所述粗DEM数据包括：SRTM DEM数据或ASTER DEM数据。
在本公开实施例中，所述步骤E在[hinit-Δhrange，hinit+Δhrange]范围内每次取高程值为hj＝hinit-Δhrange+(j-1)·δh，其中，步骤B所划分的每一个像素网格均有一个高斯坐标(Xgauss，Ygauss)，每一个像素网格的高斯坐标(Xgauss，Ygauss)处对应的粗DEM值为Zcoarse，Zcoarse为待测量高程的初值hinit，Δhrange的取值为所选择的粗DEM数据的精度，δh为高程值的变化步长，j为迭代次数， 表示向上取整。
在本公开实施例中，所述步骤F包括：步骤F1：对每个方向照射的InSAR数据，将该像素网格的高斯坐标(Xgauss，Ygauss)转换到与该方向机载InSAR飞行轨迹参数相同的地面坐标系下，得到该像素网格的地面坐标为(XP，YP，ZP)，其中ZP＝hj；步骤F2：对每个方向照射的InSAR数据，设对该像素网格成像时载机的位置坐标为(XA，YA，ZA)，即有其中(XA0，YA0，ZA0)为零时刻的载机位置，(VAX，VAY，VAZ)为载机的速度矢量，t为该像素网格的成像时刻，对机载InSAR成像时均处理到零多普勒平面，即多普勒中心频率fdc＝0，由此可列出如下的距离-多普勒方程：
距离方程：
多普勒方程：
其中，λ为波长，R1为该像素网格在主图像中的斜距，因此由多普勒方程计算：
将t代入距离方程计算出R1；步骤F3：对每个方向照射的InSAR数据，进一步由R1＝R0+ρrx，可计算出该像素网格投影到该方向照射的InSAR数据对应的主图像斜距平面中的距离向和方位向坐标(x，y)：
y＝t·PRF
其中R0为主图像的近距，ρr为主图像的距离向分辨率，PRF为脉冲重复频率，均为已知的参数；步骤F4：对每个方向照射的InSAR数据，通过插值得到(x，y)处的干涉相位值φi，j_measure和相干系数值γi，j_measure，其中i表示第i个方向照射的数据，i＝1，2，…，N，N为观测的次数，j为迭代次数。插值的方法包括：双线性插值或三次样条插值方法。
在本公开实施例中，步骤G中，对于每次迭代的高程值hj，根据各个方向多次照射的InSAR数据的成像几何，分别计算高程值为hj时理想的干涉相位φi，j_ideal，其中i表示第i个方向照射的数据，i＝1，2，…，N，N为观测的次数，j为迭代次数。φi，j_ideal的具体计算方法如下：
其中，R1和R2分别为该像素网格在主、辅图像中的斜距，R1在步骤F中计算得到，H为载机平台高度，θ为InSAR主天线的视角，B和α为InSAR基线长度和基线倾角，Q为系数，Q＝1时表示标准模式InSAR系统，Q＝2时表示乒乓模式InSAR系统，λ为波长，wrap{·}为干涉相位缠绕算子。
在本公开实施例中，所述步骤H中，对于每次迭代的高程值hj，针对各个方向照射的InSAR数据，分别计算在干涉相位观测值φi，j_measure一定的前提下关于hj的似然函数f(φi，j_measure\hj)为：
f(φi，j_measure\hj)＝pdf(φi，j_measure；γi，j_measure，φi，j_ideal)
其中，pdf(φi，j_measure；γi，j_measure，φi，j_ideal)为干涉相位期望值为φi，j_ideal时，干涉相位观测值φi，j_measure的概率密度函数，概率密度函数与该像素网格的相关系数γi，j_measure有关，计算方法如下：
其中，L为多视数，β＝|γi，j_measure|cos(φi，j_measure-φi，j_ideal)，为高斯超几何分布函数，记为：
其中，p，k分别为向量n，d的长度， 当多视数L＝1时，
在本公开实施例中，所述步骤I中，计算各个方向多次照射的InSAR数据在联合的干涉相位观测值(φ1，j_measure，φ2，j_measure，…，φN，j_measure)及关于hj的联合似然函数f(φ1，j_measure，φ2，j_measure，…，φN，j_measure\hj)，不同方向照射的InSAR数据之间相互独立，因此有：
在本公开实施例中，所述步骤J中，取使联合似然函数最大时对应的hj即为该像素网格的最终高程估计值，进而重建得到每个像素网格的DEM，完成机载InSAR困难区域的DEM重建，公式如下：
(三)有益效果
从上述技术方案可以看出，本公开机载InSAR困难区域DEM重建方法至少具有以下有益效果其中之一或其中一部分：
(1)通过将不同方向多次照射的InSAR数据进行拼接处理，实现了不同方向照射数据可测绘区域的互相补充，克服了起伏地形条件下，单次观测大面积几何畸变区域无法进行DEM重建的问题，大大提升了可测绘区域的范围，提高了机载InSAR地形测绘应用的实用化程度；
(2)对于不同方向多次照射的InSAR数据均能进行DEM重建的非几何畸变区域，通过估计联合似然函数的这一现代信号处理方法，有效地利用了多次观测之间的冗余信息，提高了DEM重建的精度；
(3)直接在统一的地理坐标系，即高斯投影坐标系下对不同方向多次照射的InSAR数据进行拼接，解决了多次照射数据无法在斜距平面直接进行拼接的问题，从而同时实现了DEM的地理编码过程和多次照射数据的拼接融合过程，避免了对多次照射的数据分别进行地理编码处理后再进行拼接处理的复杂性。
附图说明
图1为本公开实施例机载InSAR困难区域DEM重建方法流程示意图。
图2为本公开实施例机载InSAR困难区域DEM重建方法构架示意图。
图3为本公开实施例机载InSAR成像几何示意图。
具体实施方式
本公开提供了一种机载InSAR困难区域的DEM重建方法，所述机载InSAR困难区域的DEM重建方法，通过将不同方向多次照射的InSAR数据进行拼接处理，实现了不同方向照射数据可测绘区域的互相补充，克服了起伏地形条件下，单次观测大面积几何畸变区域无法进行DEM重建的问题，大大提升了可测绘区域的范围，提高了机载InSAR地形测绘应用的实用化程度。
为使本公开的目的、技术方案和优点更加清楚明白，以下结合具体实施例，并参照附图，对本公开进一步详细说明。
在本公开实施例中，提供一种机载InSAR困难区域的DEM重建方法，图1为所述机载InSAR困难区域DEM重建方法流程示意图，图2为本公开实施例机载InSAR困难区域DEM重建方法构架示意图。结合图1和图2所示，所述机载InSAR困难区域的DEM重建方法，包括如下步骤：
步骤A：对各个方向照射的InSAR数据，分别进行初步的干涉处理。
所述步骤A，包括：
步骤A1：对各个方向照射的InSAR数据生成的主、辅图像进行复图像配准；
步骤A2：对步骤A1配准后的主、辅图像进行共轭相乘，对共轭相乘的结果取相位，生成干涉相位图；
步骤A3：对步骤A1配准后的主、辅图像计算每一个像素的相干系数，生成相干系数图；
对于相同的待测绘场景，从不同的方向进行照射的InSAR数据中，所形成的几何畸变的区域不相同，因此可以通过不同方向数据的拼接实现互相补充，从而大大降低最终无法测绘的区域范围。
步骤B：将待测绘DEM的场景区域四角点的地理坐标，利用高斯投影的方法投影到高斯坐标系下，并将由此确定的矩形区域按照生成DEM产品要求的分辨率(如1m×1m)划分为等间隔的像素网格；
其中，每一个像素均有一个高斯坐标(Xgauss，Ygauss)；
步骤C：选择待测绘DEM的场景区域对应的合适的外部粗DEM数据，根据所选粗DEM数据的格式将其转换到高斯坐标系下，并将其插值到生成产品要求的分辨率，使其与步骤B所确定的等间隔的像素网格一一对应；
目前可以公开并广泛应用的粗精度DEM数据主要有SRTM(Shuttle RadarTopography Mission，航天飞机雷达地形测绘任务)DEM数据和ASTER(AdvancedSpaceborne Thermal Emission and Reflection Radiometer，先进星载热辐射与反射辐射计)DEM数据。
在本公开实施例中，每一个像素网格(Xgauss，Ygauss)处对应的高程值为Zcoarse。
步骤D：对步骤B所述的每一个像素网格，以该像素网格对应的粗DEM值Zcoarse为待测量高程的初值hinit，设置该像素网格高程值的变化范围为[hinit-Δhrange，hinit+Δhrang。]；
其中，粗DEM值Zcoarse为待测量高程的初值hinit，Δhrange的取值为所选择的粗DEM数据的精度。设置高程值的变化步长为δh，其取值为要求生成DEM产品的精度。
步骤E：对每一个像素网格，设置每次迭代处理的高程值；
在[hinit-Δhrange，hinit+Δhrange]范围内每次取高程值为hj＝hinit-Δhrange+(j-1)·δh，其中j为迭代次数， 表示向上取整。
步骤F：对每一个像素网格，对于步骤E所述的每次迭代的高程值hj，根据各个方向多次照射的InSAR数据主图像的成像几何，分别计算该像素网格投影到各个方向对应的主图像斜距平面中的距离向和方位向坐标为(x，y)，并读取干涉相位图中(x，y)处的干涉相位观测值φi，j_measure，以及相干系数图中(x，y)处的相干系数观测值γi，j_measure，其中i表示第i个方向照射的数据，i＝1，2，…，N，N为观测的次数，j为迭代次数。
在本公开实施例中，图3为机载InSAR成像几何示意图，如图3所示，(x，y)的具体计算方法如下：
步骤F1：将该像素网格的高斯坐标(Xgauss，Ygauss)转换到与机载InSAR飞行轨迹参数相同的地面坐标系下，如机载InSAR成像常用的北东天坐标系，该像素网格的地面坐标为(XP，YP，ZP)，其中ZP＝hi。
步骤F2：设对该像素网格成像时载机的位置坐标为(XA，YA，ZA)，由于机载InSAR成像时已经通过运动补偿将飞行轨迹拟合为匀速直线运动轨迹，即有其中(XA0，YA0，ZA0)为零时刻的载机位置，(VAX，VAY，VAZ)为载机的速度矢量，t为该像素网格的成像时刻。通常对机载InSAR成像时均处理到零多普勒平面，即多普勒中心频率fdc＝0，由此可列出如下的距离-多普勒方程：
距离方程：
多普勒方程：
其中，λ为波长，R1为该像素网格在主图像中的斜距，因此由多普勒方程可解出
将t代入距离方程计算出R1。
步骤F3：进一步由R1＝R0+ρrx，可计算出(x，y)：
y＝t·PRF
其中R0为主图像的近距，ρr为主图像的距离向分辨率，PRF为脉冲重复频率，均为已知的参数。
步骤F4：由于计算出的坐标(x，y)不一定为整数，因此需要通过插值得到该处的干涉相位值和相干系数值，插值的方法可选择双线性插值、三次样条插值等方法。
步骤G：对每一个像素网格，对于步骤E所述的每次迭代的高程值hj，根据各个方向多次照射的InSAR数据的成像几何，分别计算高程值为hj时理想的干涉相位φi，j_ideal，其中i表示第i个方向照射的数据，i＝1，2，…，N，N为观测的次数，j为迭代次数；
根据图3所示的机载InSAR成像几何示意图，φi，j_ideal的具体计算方法如下：
其中，R1和R2分别为该像素网格在主、辅图像中的斜距，R1在步骤F中计算得到，H为载机平台高度，θ为InSAR主天线的视角，B和α为InSAR基线长度和基线倾角，Q为系数，Q＝1时表示6标准模式InSAR系统，Q＝2时表示乒乓模式InSAR系统，λ为波长，wrap{·}为干涉相位缠绕算子。
步骤H：对每一个像素网格，对于步骤E所述的每次迭代的高程值hj，针对各个方向照射的InSAR数据，分别计算在干涉相位观测值φi，j_measure一定的前提下关于hj的似然函数f(φi，j_measure\hj)为：
所述似然函数f(φi，j_measure\hj)＝pdf(φi，j_measure；γi，j_measure，φi，j_ideal)，其中，pdf(φi，j_measure；γi，j_measure，φi，j_ideal)为干涉相位期望值为φi，j_ideal时，干涉相位观测值φi，j_measure的概率密度函数，概率密度函数与该像素网格的相关系数γi，j_measure有关，计算方法如下：
其中，L为多视数，β＝|γi，j_measure|cos(φi，j_measure-φi，j_ideal)，为高斯超几何分布函数，其定义如下：
其中，p，k分别为向量n，d的长度， 当多视数L＝1时，
步骤I：对每一个像素网格，对于步骤E所述的每次迭代的高程值，计算各个方向多次照射的InSAR数据在联合的干涉相位观测值(φ1，j_measure，φ2，j_measure，…，φN，j_measure)一定的前提下关于hj的联合似然函数f(φ1，j_measure，φ2，j_measure，…，φN，j_measure\hj)。
由于可以认为不同方向照射的InSAR数据之间相互独立，因此有，
步骤J：对每一个像素网格，迭代次数j每次增加1，重复步骤E至步骤I，取使联合似然函数最大时对应的hj即为该像素网格的最终高程估计值，进而重建得到步骤E所述的每个像素网格的DEM，完成机载InSAR困难区域的DEM重建。
至此，即重建得到了各个像素网格的DEM，也就是得到整个待测绘场景区域在高斯投影坐标系下的DEM。
至此，已经结合附图对本公开实施例进行了详细描述。需要说明的是，在附图或说明书正文中，未绘示或描述的实现方式，均为所属技术领域中普通技术人员所知的形式，并未进行详细说明。此外，上述对各元件和方法的定义并不仅限于实施例中提到的各种具体结构、形状或方式，本领域普通技术人员可对其进行简单地更改或替换。
依据以上描述，本领域技术人员应当对本公开机载InSAR困难区域DEM重建方法有了清楚的认识。
综上所述，本公开提供了一种机载InSAR困难区域的DEM重建方法，所述机载InSAR困难区域的DEM重建方法，通过将不同方向多次照射的InSAR数据进行拼接处理，实现了不同方向照射数据可测绘区域的互相补充，克服了起伏地形条件下，单次观测大面积几何畸变区域无法进行DEM重建的问题，大大提升了可测绘区域的范围，提高了机载InSAR地形测绘应用的实用化程度。
还需要说明的是，实施例中提到的方向用语，例如“上”、“下”、“前”、“后”、“左”、“右”等，仅是参考附图的方向，并非用来限制本公开的保护范围。贯穿附图，相同的元素由相同或相近的附图标记来表示。在可能导致对本公开的理解造成混淆时，将省略常规结构或构造。
并且图中各部件的形状和尺寸不反映真实大小和比例，而仅示意本公开实施例的内容。另外，在权利要求中，不应将位于括号之间的任何参考符号构造成对权利要求的限制。
除非有所知名为相反之意，本说明书及所附权利要求中的数值参数是近似值，能够根据通过本公开的内容所得的所需特性改变。具体而言，所有使用于说明书及权利要求中表示组成的含量、反应条件等等的数字，应理解为在所有情况中是受到「约」的用语所修饰。一般情况下，其表达的含义是指包含由特定数量在一些实施例中±10％的变化、在一些实施例中±5％的变化、在一些实施例中±1％的变化、在一些实施例中±0.5％的变化。
再者，单词“包含”不排除存在未列在权利要求中的元件或步骤。位于元件之前的单词“一”或“一个”不排除存在多个这样的元件。
说明书与权利要求中所使用的序数例如“第一”、“第二”、“第三”等的用词，以修饰相应的元件，其本身并不意味着该元件有任何的序数，也不代表某一元件与另一元件的顺序、或是制造方法上的顺序，该些序数的使用仅用来使具有某命名的一元件得以和另一具有相同命名的元件能做出清楚区分。
此外，除非特别描述或必须依序发生的步骤，上述步骤的顺序并无限制于以上所列，且可根据所需设计而变化或重新安排。并且上述实施例可基于设计及可靠度的考虑，彼此混合搭配使用或与其他实施例混合搭配使用，即不同实施例中的技术特征可以自由组合形成更多的实施例。
本领域那些技术人员可以理解，可以对实施例中的设备中的模块进行自适应性地改变并且把它们设置在与该实施例不同的一个或多个设备中。可以把实施例中的模块或单元或组件组合成一个模块或单元或组件，以及此外可以把它们分成多个子模块或子单元或子组件。除了这样的特征和/或过程或者单元中的至少一些是相互排斥之外，可以采用任何组合对本说明书(包括伴随的权利要求、摘要和附图)中公开的所有特征以及如此公开的任何方法或者设备的所有过程或单元进行组合。除非另外明确陈述，本说明书(包括伴随的权利要求、摘要和附图)中公开的每个特征可以由提供相同、等同或相似目的的替代特征来代替。并且，在列举了若干装置的单元权利要求中，这些装置中的若干个可以是通过同一个硬件项来具体体现。
类似地，应当理解，为了精简本公开并帮助理解各个公开方面中的一个或多个，在上面对本公开的示例性实施例的描述中，本公开的各个特征有时被一起分组到单个实施例、图、或者对其的描述中。然而，并不应将该公开的方法解释成反映如下意图：即所要求保护的本公开要求比在每个权利要求中所明确记载的特征更多的特征。更确切地说，如下面的权利要求书所反映的那样，公开方面在于少于前面公开的单个实施例的所有特征。因此，遵循具体实施方式的权利要求书由此明确地并入该具体实施方式，其中每个权利要求本身都作为本公开的单独实施例。
以上所述的具体实施例，对本公开的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本公开的具体实施例而已，并不用于限制本公开，凡在本公开的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本公开的保护范围之内。

Claims (9)

1.一种机载InSAR困难区域的DEM重建方法，包括：
步骤A：对各个方向照射的InSAR数据，分别进行初步的干涉处理，生成主图像和辅图像、干涉相位图以及相干系数图；
步骤B：将待测绘DEM的场景区域四角点的地理坐标，利用高斯投影的方法投影到高斯坐标系下，并将高斯投影后高斯坐标系下的矩形区域按照生成DEM产品要求的分辨率划分为等间隔的像素网格；
步骤C：选择待测绘DEM的场景区域对应的外部粗DEM数据，根据所选粗DEM数据的格式将粗DEM数据转换到高斯坐标系下，并将所述粗DEM数据插值到生成产品要求的分辨率，使粗DEM数据与步骤B所确定的等间隔的像素网格一一对应；
步骤D：对步骤B所划分的每一个像素网格，以该像素网格对应的粗DEM值为待测量高程的初值，设置该像素网格高程值的变化范围；
步骤E：对步骤D所处理后的每一个像素网格，设置每一次迭代处理的高程值；
步骤F：对于步骤E所述的每次迭代的高程值，根据各个方向多次照射的InSAR数据主图像的成像几何，分别计算该像素网格投影到各个方向对应的主图像斜距平面中的距离向和方位向坐标，并读取步骤A所生成的干涉相位图中该像素网格对应坐标处的干涉相位观测值，以及步骤A所生成的相干系数图中该像素网格对应坐标处的相干系数观测值；
步骤G：对于步骤E所述的每次迭代的高程值，根据各个方向多次照射的InSAR数据的成像几何，分别计算所述高程值对应的理想的干涉相位；
步骤H：对于步骤E所述的每次迭代的高程值，针对各个方向照射的InSAR数据，分别计算在干涉相位观测值一定的前提下关于所述高程值的似然函数；
步骤I：对于步骤E所述的每次迭代的高程值，计算各个方向多次照射的InSAR数据在联合的干涉相位观测值一定的前提下关于迭代高程值的联合似然函数；以及
步骤J：迭代次数每次增加1，重复步骤E至步骤I，取使联合似然函数最大时对应的迭代高程值即为该像素网格的最终高程估计值，即重建得到了步骤E所述的每个像素网格的DEM，完成机载InSAR困难区域的DEM重建。
2.根据权利要求1所述的机载InSAR困难区域的DEM重建方法，其中，所述步骤A，包括：
步骤A1：对各个方向照射的InSAR数据生成的主、辅图像进行复图像配准；
步骤A2：对步骤A1配准后的主、辅图像进行共轭相乘，对共轭相乘的结果取相位，生成干涉相位图；以及
步骤A3：对步骤A1配准后的主、辅图像计算每一个像素的相干系数，生成相干系数图。
3.如权利要求1所述的机载InSAR困难区域的DEM重建方法，其中，步骤C所述粗DEM数据包括：SRTM DEM数据或ASTER DEM数据。
4.根据权利要求1所述的机载InSAR困难区域的DEM重建方法，其中，所述步骤E在[hinit-Δhrange，hinit+Δhrange]范围内每次取高程值为hj＝hinit-Δhrange+(j-1)·δh，其中，步骤B所划分的每一个像素网格均有一个高斯坐标(Xgauss，Ygauss)，每一个像素网格的高斯坐标(Xgauss，Ygauss)处对应的粗DEM值为Zcoarse，Zcoarse为待测量高程的初值hinit，Δhrange的取值为所选择的粗DEM数据的精度，δh为高程值的变化步长，j为迭代次数， 表示向上取整。
5.根据权利要求1所述的机载InSAR困难区域的DEM重建方法，其中，所述步骤F包括：
步骤F1：对每个方向照射的InSAR数据，将该像素网格的高斯坐标(Xgauss，Ygauss)转换到与该方向机载InSAR飞行轨迹参数相同的地面坐标系下，得到该像素网格的地面坐标为(XP，YP，ZP)，其中ZP＝hj；
步骤F2：对每个方向照射的InSAR数据，设对该像素网格成像时载机的位置坐标为(XA，YA，ZA)，即有其中(XA0，YA0，ZA0)为零时刻的载机位置，(VAX，VAY，VAZ)为载机的速度矢量，t为该像素网格的成像时刻，对机载InSAR成像时均处理到零多普勒平面，即多普勒中心频率fdc＝0，由此可列出如下的距离-多普勒方程：
距离方程：
多普勒方程：
其中，λ为波长，R1为该像素网格在主图像中的斜距，因此由多普勒方程计算：
将t代入距离方程计算出R1；
步骤F3：对每个方向照射的InSAR数据，进一步由R1＝R0+ρrx，可计算出该像素网格投影到该方向照射的InSAR数据对应的主图像斜距平面中的距离向和方位向坐标(x，y)：
y＝t·PRF
其中R0为主图像的近距，ρr为主图像的距离向分辨率，PRF为脉冲重复频率，均为已知的参数；
步骤F4：对每个方向照射的InSAR数据，通过插值得到(x，y)处的干涉相位值φi，j_measure和相干系数值γi，j_measure，其中i表示第i个方向照射的数据，i＝1，2，…，N，N为观测的次数，j为迭代次数。插值的方法包括：双线性插值或三次样条插值方法。
6.根据权利要求1所述的机载InSAR困难区域的DEM重建方法，其中，步骤G中，对于每次迭代的高程值hj，根据各个方向多次照射的InSAR数据的成像几何，分别计算高程值为hj时理想的干涉相位φi，j_ideal，其中i表示第i个方向照射的数据，i＝1，2，…，N，N为观测的次数，j为迭代次数。φi，j_ideal的具体计算方法如下：
其中，R1和R2分别为该像素网格在主、辅图像中的斜距，R1在步骤F中计算得到，H为载机平台高度，θ为InSAR主天线的视角，B和α为InSAR基线长度和基线倾角，Q为系数，Q＝1时表示标准模式InSAR系统，Q＝2时表示乒乓模式InSAR系统，λ为波长，wrap{·}为干涉相位缠绕算子。
7.根据权利要求1所述的机载InSAR困难区域的DEM重建方法，其中，所述步骤H中，对于每次迭代的高程值hj，针对各个方向照射的InSAR数据，分别计算在干涉相位观测值φi，j_measure一定的前提下关于hj的似然函数f(φi，j_measure\hj)为：
f(φi，j_measure\hj)＝pdf(φi，j_measure；γi，j_measure，φi，j_ideal)
其中，pdf(φi，j_measure；γi，j_measure，φi，j_ideal)为干涉相位期望值为φi，j_ideal时，干涉相位观测值φi，j_measure的概率密度函数，概率密度函数与该像素网格的相关系数γi，j_measure有关，计算方法如下：
其中，L为多视数，为高斯超几何分布函数，记为：
其中，p，k分别为向量n，d的长度， 当多视数L＝1时，
8.根据权利要求1所述的机载InSAR困难区域的DEM重建方法，其中，所述步骤I中，计算各个方向多次照射的InSAR数据在联合的干涉相位观测值(φ1，j_measure，φ2，j_measure，…，φN，j_measure)及关于hj的联合似然函数f(φ1，j_measure，φ2，j_measure，…，φN，j_measure\hj)，不同方向照射的InSAR数据之间相互独立，因此有：
9.根据权利要求1所述的机载InSAR困难区域的DEM重建方法，其中，所述步骤J中，取使联合似然函数最大时对应的hj即为该像素网格的最终高程估计值，进而重建得到每个像素网格的DEM，完成机载InSAR困难区域的DEM重建，公式如下：
CN201810919415.7A 2018-08-13 2018-08-13 机载InSAR困难区域的DEM重建方法 Pending CN109085587A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810919415.7A CN109085587A (zh) 2018-08-13 2018-08-13 机载InSAR困难区域的DEM重建方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810919415.7A CN109085587A (zh) 2018-08-13 2018-08-13 机载InSAR困难区域的DEM重建方法

Publications (1)

Publication Number Publication Date
CN109085587A true CN109085587A (zh) 2018-12-25

Family

ID=64834637

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810919415.7A Pending CN109085587A (zh) 2018-08-13 2018-08-13 机载InSAR困难区域的DEM重建方法

Country Status (1)

Country Link
CN (1) CN109085587A (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110068817A (zh) * 2019-05-07 2019-07-30 中国科学院电子学研究所 一种基于激光测距和InSAR的地形测图方法、仪器和系统
CN110703252A (zh) * 2019-11-11 2020-01-17 中国科学院电子学研究所 干涉合成孔径雷达阴影区域数字高程模型修正方法
CN111681299A (zh) * 2020-06-16 2020-09-18 中国科学院空天信息创新研究院 基于InSAR解缠相位生成数字表面模型的方法及装置
CN111983608A (zh) * 2020-07-07 2020-11-24 西安瑞得空间信息技术有限公司 一种快速dem生成方法

Non-Patent Citations (6)

* Cited by examiner, † Cited by third party
Title
M. EINEDER等: "A maximum-likelihood estimator to simultaneously unwrap, geocode, and fuse SAR interferograms from different viewing geometries into one digital elevation model", 《 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING》 *
刘华有等: "一种有效的机载双频干涉地形高程重建方法", 《雷达学报》 *
张秋玲等: "利用多基线数据融合提高分布式卫星InSAR系统的干涉相位精度", 《电子与信息学报》 *
方东生等: "一种机载重轨干涉SAR相位偏移量估计算法", 《雷达科学与技术》 *
蔡斌等: "基于分布式星载单基线InSAR的叠掩和阴影区域检测", 《信号处理》 *
袁志辉等: "改进的基于最大似然估计的多通道InSAR高程重建方法", 《电子与信息学报》 *

Cited By (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110068817A (zh) * 2019-05-07 2019-07-30 中国科学院电子学研究所 一种基于激光测距和InSAR的地形测图方法、仪器和系统
CN110703252A (zh) * 2019-11-11 2020-01-17 中国科学院电子学研究所 干涉合成孔径雷达阴影区域数字高程模型修正方法
CN110703252B (zh) * 2019-11-11 2021-07-09 中国科学院电子学研究所 干涉合成孔径雷达阴影区域数字高程模型修正方法
CN111681299A (zh) * 2020-06-16 2020-09-18 中国科学院空天信息创新研究院 基于InSAR解缠相位生成数字表面模型的方法及装置
CN111681299B (zh) * 2020-06-16 2024-03-15 中国科学院空天信息创新研究院 基于InSAR解缠相位生成数字表面模型的方法及装置
CN111983608A (zh) * 2020-07-07 2020-11-24 西安瑞得空间信息技术有限公司 一种快速dem生成方法

Similar Documents

Publication Publication Date Title
CN109085587A (zh) 机载InSAR困难区域的DEM重建方法
Wang et al. Bistatic SAR system and signal processing technology
Chen et al. A backprojection-based imaging for circular synthetic aperture radar
JP6421395B2 (ja) Ｓａｒ図からの立体地形図形成方法
CN103543453B (zh) 一种地球同步轨道合成孔径雷达干涉的高程反演方法
CN106526593B (zh) 基于sar严密成像模型的子像素级角反射器自动定位方法
CN109800380A (zh) 星载微波遥感仪器对地探测的严密成像几何模型构建方法
CN109270527B (zh) 圆迹sar子孔径图像序列联合相关dem提取方法
KR100441590B1 (ko) 간섭측정용 합성 개구 레이다의 기하학적 특성을 이용하여지형고도를 측정하기 위한 디지털 고도모형 생성방법
CN103176171B (zh) 一种干涉sar分布目标仿真方法
Alexandrov et al. Derivation of cumulus cloud dimensions and shape from the airborne measurements by the Research Scanning Polarimeter
Liu et al. Focusing challenges of ships with oscillatory motions and long coherent processing interval
Durand et al. Qualitative assessment of four DSM generation approaches using Pléiades-HR data
CN107504919B (zh) 基于相位映射的折叠相位三维数字成像方法及装置
Qiao et al. Assessment of geo-positioning capability of high resolution satellite imagery for densely populated high buildings in metropolitan areas
Caduff et al. Registration and visualisation of deformation maps from terrestrial radar interferometry using photogrammetry and structure from Motion
CN108008367B (zh) 星载单航过InSAR系统电离层误差校正方法
CN115685200A (zh) 一种高精度大前斜视sar成像运动补偿与几何校正方法
Li et al. An advanced DSS-SAR InSAR terrain height estimation approach based on baseline decoupling
Deo et al. Evaluation of interferometric SAR DEMs generated using TanDEM-X data
Huang et al. 3D building reconstruction and visualization for single high resolution satellite image
Kuang et al. Fully three-dimensional UAV SAR imaging with multi-azimuth-angle observation
CN103197293A (zh) 一种机载干涉sar多路径误差的提取与补偿方法
CN111696207B (zh) 一种基于引导滤波的多基线dem融合方法
CARCANO Merging local DTMS: HELI-DEM project, problems and solutions

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
RJ01 Rejection of invention patent application after publication
RJ01 Rejection of invention patent application after publication

Application publication date: 20181225