机载InSAR困难区域的DEM重建方法与流程

本公开涉及干涉合成孔径雷达信号处理领域，尤其涉及一种机载insar困难区域的dem重建方法。

背景技术

干涉合成孔径雷达(interferometricsyntheticapertureradar，insar)技术将无线电干涉测量技术与sar(合成孔径雷达)技术相结合，在获取二维sar图像的同时，能够利用sar复数据的相位信息提取地表的三维信息和变化信息，从而将sar的测量空间从二维拓展到三维，随着insar系统的迅速发展，insar技术的应用领域也得以不断扩展，其中，地形测绘一直是insar技术最直接和最主要的应用之一，与摄影测量、激光雷达等技术相比，insar具有全天时、全天候的优点，与雷达立体测量技术相比，insar具有更高的测量精度。

然而，由于sar侧视成像的特点，使得sar图像中存在几何畸变现象，对于地形起伏较大的城市、山区等区域而言，这一现象尤为严重，影响insar地形测绘的几何畸变现象包括阴影和叠掩，阴影区域由于雷达接收不到地表的有用回波信息，相应的干涉相位表现为噪声特性，其中并不包含区域本身的真实地形信息，无法用于高程反演；而叠掩区域则是不同高度区域的回波投影到同一个距离-多普勒单元内形成，其干涉相位反映的是多个散射源矢量叠加的结果，常规的insar处理方法无法分辨同一采样单元内的多个散射源的相位信息，从而也无法准确反演其高程。由此可见，在起伏的地形条件下，常规的干涉处理方法无法实现阴影、叠掩区域的高程测量，这些区域是insar处理的困难区域，也是限制insar地形测绘应用实用化的关键。

目前，对于insar困难区域的处理方法主要包括以下两个方面：首先，对于叠掩区域的处理，现有技术主要利用多基线insar技术实现多个散射源的检测分辨和高程估计，然而，这类方法要达到较高的估计精度，需要较多的基线数，因而在实际系统中难以实现。其次，对于阴影区域的处理，现有技术主要集中在对其的检测识别上，在检测的基础上，一方面，在干涉处理过程中，对这些区域进行掩膜处理，从而避免相位解缠时误差的传播；另一方面，利用外源的dem(digitalelevationmodel，数字高程模型)信息对检测出的区域进行补充，但是，可以公开获取的外源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数据包括：srtmdem数据或asterdem数据。

在本公开实施例中，所述步骤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(shuttleradartopographymission，航天飞机雷达地形测绘任务)dem数据和aster(advancedspacebornethermalemissionandreflectionradiometer，先进星载热辐射与反射辐射计)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％的变化。

再者，单词“包含”不排除存在未列在权利要求中的元件或步骤。位于元件之前的单词“一”或“一个”不排除存在多个这样的元件。

说明书与权利要求中所使用的序数例如“第一”、“第二”、“第三”等的用词，以修饰相应的元件，其本身并不意味着该元件有任何的序数，也不代表某一元件与另一元件的顺序、或是制造方法上的顺序，该些序数的使用仅用来使具有某命名的一元件得以和另一具有相同命名的元件能做出清楚区分。

此外，除非特别描述或必须依序发生的步骤，上述步骤的顺序并无限制于以上所列，且可根据所需设计而变化或重新安排。并且上述实施例可基于设计及可靠度的考虑，彼此混合搭配使用或与其他实施例混合搭配使用，即不同实施例中的技术特征可以自由组合形成更多的实施例。

本领域那些技术人员可以理解，可以对实施例中的设备中的模块进行自适应性地改变并且把它们设置在与该实施例不同的一个或多个设备中。可以把实施例中的模块或单元或组件组合成一个模块或单元或组件，以及此外可以把它们分成多个子模块或子单元或子组件。除了这样的特征和/或过程或者单元中的至少一些是相互排斥之外，可以采用任何组合对本说明书(包括伴随的权利要求、摘要和附图)中公开的所有特征以及如此公开的任何方法或者设备的所有过程或单元进行组合。除非另外明确陈述，本说明书(包括伴随的权利要求、摘要和附图)中公开的每个特征可以由提供相同、等同或相似目的的替代特征来代替。并且，在列举了若干装置的单元权利要求中，这些装置中的若干个可以是通过同一个硬件项来具体体现。

类似地，应当理解，为了精简本公开并帮助理解各个公开方面中的一个或多个，在上面对本公开的示例性实施例的描述中，本公开的各个特征有时被一起分组到单个实施例、图、或者对其的描述中。然而，并不应将该公开的方法解释成反映如下意图：即所要求保护的本公开要求比在每个权利要求中所明确记载的特征更多的特征。更确切地说，如下面的权利要求书所反映的那样，公开方面在于少于前面公开的单个实施例的所有特征。因此，遵循具体实施方式的权利要求书由此明确地并入该具体实施方式，其中每个权利要求本身都作为本公开的单独实施例。

以上所述的具体实施例，对本公开的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本公开的具体实施例而已，并不用于限制本公开，凡在本公开的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本公开的保护范围之内。