一种基于导航卫星双基地InSAR的DEM成像方法及装置

一种基于导航卫星双基地insar的dem成像方法及装置
技术领域
1.本发明涉及双基地合成孔径雷达技术领域，具体涉及一种基于导航卫星双基地insar的dem成像方法及装置。


背景技术：

2.干涉合成孔径雷达(insar，interferometric synthetic aperture radar)可用于高程反演，需要借助外源dem(digital elevation model)用于估计和补偿干涉相位中的平地相位和地形相位，从而使干涉条纹更加稀疏，有利于后续相位解缠和高精度高程反演，以进行形变监测。
3.基于导航卫星的双基地insar系统，与传统机载、星载或地基系统相比，主要具有1)回波的低信噪比；2)距离向的低分辨率；3)系统几何构型的非对称性这三大特征，导致其成像结果与传统系统存在很大差别，主要表现为强散射点的点扩展函数(psf，point spread function)的联合，分辨率差，缺乏纹理信息。将该系统应用于形变反演时，如果根据传统系统所采取的方法，即直接引入外源dem作为成像平面进行处理时，如果场景高程起伏较复杂，dem精度较高，会使成像结果中永久散射体(ps，persistent scatterer)的分辨单元的幅度和相位跳变、边缘畸变、形状不连续，从而不利于后续ps点的提取和配准，影响最终的形变反演精度。


技术实现要素：

4.有鉴于此，本发明提供了一种基于导航卫星双基地insar的dem成像方法及装置，能够在目标场景高程起伏较复杂、外源dem精度较高的情况下，使 ps点分辨单元内幅度和相位连续，保留其psf边缘理论特性，保证其形状连续，提高了ps点的时间相关性，提升了系统的配准精度。
5.为了解决上述技术问题，本发明是这样实现的。
6.一种基于导航卫星双基地insar的dem成像方法，所述方法包括以下步骤：
7.步骤s1：使用滑动最小二乘拟合方法拟合本地坐标系下的dem；
8.步骤s2：对拟合后的dem的各项进行二次差分，找到所述拟合后的dem 存在的所有不连续单元；
9.步骤s3：对每个拟合后的不连续单元，以所述拟合后的不连续单元为中心进行边界扩展，扩展后的结果作为与所述拟合后的不连续单元相对应的待平滑核；
10.步骤s4：对每个待平滑核，生成与所述待平滑核对应的、以所述待平滑核为基准进行中心对称变换后的高斯平滑核；
11.步骤s5：对每个待平滑核，将其与对应的高斯平滑核进行加权求和运算，进而得到每个待平滑核对应的不连续单元的拟合平滑结果；
12.步骤s6：将所有不连续单元的拟合平滑结果作为最终平滑后的dem，将所述最终平滑后的dem作为后向投影成像平面进行成像，得到最终的成像结果。
13.优选地，所述步骤s1：使用滑动最小二乘拟合方法拟合本地坐标系下的 dem，包括：
14.步骤s11：dem的n个单元，记为f
t
＝[f(x1,y1) f(x2,y2)
ꢀ…ꢀ
f(xn,yn)]，
[0015]
其中，f(xi,yi)是dem的第i单元，i＝1,2,
…
,n；xi为第i单元的横坐标，yi为第i单元的纵坐标；
[0016]
步骤s12：对于dem的n个单元中的每一个f(xi,yi)，采用二阶加权最小二乘拟合，包括：对于第i单元f(xi,yi)，设基函数为b
t
＝[1 x y x
2 xy y2]，假设拟合该单元的参考单元个数为m个,记f
t
＝[f(x1,y1) f(x2,y2)
ꢀ…ꢀ
f(xm,ym)]，每个参考单元的权重服从高斯分布，记为：
[0017][0018]
则根据加权最小二乘计算，第i单元f(xi,yi)的拟合结果为：
[0019][0020]
其中，
[0021][0022]
其中，xi和yi是基函数bi的两个变量，wi(xj,yj)为第j个参考单元的加权系数，1≤j≤m；f
t
为参与第i单元拟合的所有参考单元的真实dem构成的矩阵的转置，f为参与第i单元拟合的所有参考单元的真实dem构成的矩阵；
[0023]
步骤s13：基于dem的n个单元的拟合结果，得到dem全局拟合结果
[0024]
优选地，所述步骤s2，包括：对所述dem全局拟合结果进行二次差分处理，并以所述dem全局拟合结果的均值作为阈值；对dem的n个单元进行遍历，找到二次差分结果大于所述阈值的单元，确定为拟合不连续单元；
[0025]
所述二次差分处理为：求所述dem全局拟合结果的二阶偏导数，对x方向求偏导或对y方向求偏导。
[0026]
优选地，所述步骤s4，所述生成与所述待平滑核对应的、以所述待平滑核为基准进行中心对称变换后的高斯平滑核的方法为：产生一个与待平滑核大小相等的矩阵，该矩阵内的数值服从二维高斯分布，然后以该矩阵中心元素为中心，对该矩阵所有元素进行中心对称变换，产生一个中心值最小，周围值逐渐增大，且服从高斯分布的矩阵，即为所需要的高斯平滑核。
[0027]
本发明所提供的一种基于导航卫星双基地insar的dem成像装置，包括：
[0028]
拟合模块：配置为使用滑动最小二乘拟合方法拟合本地坐标系下的dem；
[0029]
不连贯单元确定模块：配置为对拟合后的dem的各项进行二次差分，找到所述拟合后的dem存在的所有不连续单元；
[0030]
扩展模块：配置为对每个拟合后的不连续单元，以所述拟合后的不连续单元为中心进行边界扩展，扩展后的结果作为与所述拟合后的不连续单元相对应的待平滑核；
[0031]
高斯平滑核确定模块：配置为对每个待平滑核，生成与所述待平滑核对应的、以所述待平滑核为基准进行中心对称变换后的高斯平滑核；
[0032]
平滑结果生成模块：配置为对每个待平滑核，将其与对应的高斯平滑核进行加权求和运算，进而得到每个待平滑核对应的不连续单元的拟合平滑结果；
[0033]
成像结果确定模块：配置为将所有不连续单元的拟合平滑结果作为最终平滑后的dem，将所述最终平滑后的dem作为后向投影成像平面进行成像，得到最终的成像结果。
[0034]
优选地，所述拟合模块，包括：
[0035]
第一拟合子模块：配置为dem的n个单元，记为f
t
＝[f(x1,y1)f(x2,y2)
…
f(xn,yn)]，
[0036]
其中，f(xi,yi)是dem的第i单元，i＝1,2,
…
,n；xi为第i单元的横坐标，yi为第i单元的纵坐标；
[0037]
第二拟合子模块：配置为对于dem的n个单元中的每一个f(xi,yi)，采用二阶加权最小二乘拟合，包括：对于第i单元f(xi,yi)，设基函数为b
t
＝[1xyx2xyy2]，假设拟合该单元的参考单元个数为m个,记f
t
＝[f(x1,y1)f(x2,y2)
…
f(xm,ym)]，每个参考单元的权重服从高斯分布，记为：
[0038][0039]
则根据加权最小二乘计算，第i单元f(xi,yi)的拟合结果为：
[0040][0041]
其中，
[0042][0043]
其中，xi和yi是基函数bi的两个变量，wi(xj,yj)为第j个参考单元的加权系数，1≤j≤m；f
t
为参与第i单元拟合的所有参考单元的真实dem构成的矩阵的转置，f为参与第i单元拟合的所有参考单元的真实dem构成的矩阵；
[0044]
第三拟合子模块：配置为基于dem的n个单元的拟合结果，得到dem全局拟合结果
[0045]
优选地，所述不连贯单元确定模块，包括：对所述dem全局拟合结果进行二次差分
处理，并以所述dem全局拟合结果的均值作为阈值；对dem的n 个单元进行遍历，找到二次差分结果大于所述阈值的单元，确定为拟合不连续单元；
[0046]
所述二次差分处理为：求所述dem全局拟合结果的二阶偏导数，对x方向求偏导或对y方向求偏导。
[0047]
优选地，所述高斯平滑核确定模块，所述生成与所述待平滑核对应的、以所述待平滑核为基准进行中心对称变换后的高斯平滑核的方式为：产生一个与待平滑核大小相等的矩阵，该矩阵内的数值服从二维高斯分布，然后以该矩阵中心元素为中心，对该矩阵所有元素进行中心对称变换，产生一个中心值最小，周围值逐渐增大，且服从高斯分布的矩阵，即为所需要的高斯平滑核。
[0048]
有益效果：该发明解决了导航卫星双基地insar系统在进行复杂场景成像时，直接引入外源高精度dem信息所带来的ps点分辨单元内部幅度和相位跳变、边缘畸变、形状不连续问题，提高了ps点的时间相关性，有利于提高该系统应用于形变反演的精度。
附图说明
[0049]
图1为本发明所举实施例的基于导航卫星双基地insar的dem成像方法流程示意图；
[0050]
图2为本发明所举实施例的导航卫星双基地sar系统构型示意图；
[0051]
图3为本发明所举实施例的高斯对称变换二维示意图。
[0052]
图4为本发明所举实施例的拟合边缘不连续平滑处理流程图。
[0053]
图5为本发明所举实施例的导航卫星双基地insar系统的复杂场景高精度 dem成像方法。
[0054]
图6为本发明所举实施例的仿真原始dem直接成像结果。
[0055]
图7为本发明所举实施例的实测场景下的原始地形。
[0056]
图8为本发明所举实施例的相位误差幅度拟合结果。
[0057]
图9为本发明所举实施例的实测场景下的dem拟合结果。
[0058]
图10为本发明所举实施例的实测场景下的dem拟合和平滑结果。
[0059]
图11为本发明所举实施例的实测场景下原始dem成像结果。
[0060]
图12为本发明所举实施例的实测场景下拟合dem成像结果。
[0061]
图13为本发明所举实施例的实测场景下拟合和平滑处理后的dem成像结果。
具体实施方式
[0062]
下面结合附图和实施例，对本发明进行详细描述。
[0063]
本发明提出了一种基于导航卫星双基地insar的dem成像方法及装置，该方法对本地坐标系下的dem进行滑动最小二乘二阶拟合后，做二次差分，以拟合后的dem的均值作为阈值，各单元二次差分结果大于该阈值即判定为拟合后的dem的拟合不连续单元。以每个拟合不连续单元为中心，结合dem覆盖范围和不连续单元所占比重，确定扩展单元数量，从而获得所有拟合不连续单元对应的矩形待平滑核。以待平滑核大小为准，生成与之对应的进行中心对称变换后的高斯平滑核。对所有待平滑核和它所对应的高斯平滑核进行加权求和运算，所得结果作为待平滑核对应不连续单元的拟合平滑结果。将该结果作为后向投影成
像平面进行成像，最终获得ps点分辨单元内部幅度和相位无跳变、边缘不畸变、形状连续的成像结果。
[0064]
本发明提供的基于导航卫星双基地insar的dem成像方法，如图1-2所示，包括如下步骤：
[0065]
步骤s1：使用滑动最小二乘拟合方法拟合本地坐标系下的dem；
[0066]
步骤s2：对拟合后的dem的各项进行二次差分，找到所述拟合后的dem 存在的所有不连续单元；
[0067]
步骤s3：对每个拟合后的不连续单元，以所述拟合后的不连续单元为中心进行边界扩展，扩展后的结果作为与所述拟合后的不连续单元相对应的待平滑核；
[0068]
步骤s4：对每个待平滑核，生成与所述待平滑核对应的、以所述待平滑核为基准进行中心对称变换后的高斯平滑核；
[0069]
步骤s5：对每个待平滑核，将其与对应的高斯平滑核进行加权求和运算，进而得到每个待平滑核对应的不连续单元的拟合平滑结果；
[0070]
步骤s6：将所有不连续单元的拟合平滑结果作为最终平滑后的dem，将所述最终平滑后的dem作为后向投影成像平面进行成像，得到最终的成像结果。
[0071]
所述步骤s1：使用滑动最小二乘拟合方法拟合本地坐标系下的dem，包括：
[0072]
步骤s11：dem的n个单元，记为f
t
＝[f(x1,y1) f(x2,y2)
ꢀ…ꢀ
f(xn,yn)]，
[0073]
其中，f(xi,yi)是dem的第i单元，i＝1,2,
…
,n；xi为第i单元的横坐标，yi为第i单元的纵坐标；
[0074]
步骤s12：对于dem的n个单元中的每一个f(xi,yi)，采用二阶加权最小二乘拟合，包括：对于第i单元f(xi,yi)，设基函数为b
t
＝[1 x y x
2 xy y2]，假设拟合该单元的参考单元个数为m个,记f
t
＝[f(x1,y1) f(x2,y2)
ꢀ…ꢀ
f(xm,ym)]，每个参考单元的权重服从高斯分布，记为：
[0075][0076]
则根据加权最小二乘计算，第i单元f(xi,yi)的拟合结果为：
[0077][0078]
其中，
[0079][0080]
其中，xi和yi是基函数bi的两个变量，wi(xj,yj)为第j个参考单元的加权系数，1≤j≤m；f
t
为参与第i单元拟合的所有参考单元的真实dem构成的矩阵的转置，f为参与第i单元拟合的所有参考单元的真实dem构成的矩阵。
[0081]
步骤s13：基于dem的n个单元的拟合结果，得到dem全局拟合结果
[0082]
本发明得到的dem全局拟合结果不仅能保留局部dem的分布特性，还能减小相邻dem数值强跳变情况，可以在保证精度的前提下，实现对复杂场景下高精度dem的拟合效果。
[0083]
然而，由于该方法在进行局部最小二乘拟合时，不能考虑到参考点外的数值分布情况，当全场景dem变化范围较大时，会存在局部拟合边缘不连续问题，从而造成成像结果的不连续，因此，还需引入平滑处理。
[0084]
所述步骤s2，包括：对所述dem全局拟合结果进行二次差分处理，并以所述dem全局拟合结果的均值作为阈值；对dem的n个单元进行遍历，找到二次差分结果大于所述阈值的单元，确定为拟合不连续单元。所述二次差分处理为：求所述dem全局拟合结果的二阶偏导数，可以对x方向求偏导，也可以对y方向求偏导，目的是找到相邻dem变化较大的位置，即为拟合结果中不连续的位置。
[0085]
所述步骤s3：以原始dem分辨率为参考依据，以所述拟合后的不连续单元为中心，作为待平滑单元向四周扩展能确保所得待平滑核内存在连续单元的距离，作为与所述拟合后的不连续单元相对应的待平滑核。
[0086]
所述步骤s4，如图3所示，直接引入高斯平滑核在不连续点进行加权求和平滑处理时，由于常规高斯平滑核的中心权值最大，周围权值呈现高斯分布，如图3左侧所示，在实际处理时，平滑效果较差，因此，将该高斯平滑核进行中心对称变换，即中心权值最小，边缘权值最大，中间部分呈现对称高斯分布，如图 3右侧所示，即为所需要的高斯平滑核。经过处理后的高斯平滑核能够充分考虑拟合边缘不连续像素点周围数据，平滑效果更好。
[0087]
所述生成与所述待平滑核对应的、以所述待平滑核为基准进行中心对称变换后的高斯平滑核的方法为：产生一个与待平滑核大小相等的矩阵，该矩阵内的数值服从二维高斯分布，然后以该矩阵中心元素为中心，对该矩阵所有元素进行中心对称变换，产生一个中心值最小，周围值逐渐增大，且服从高斯分布的矩阵，即为所需要的高斯平滑核。
[0088]
本实施例中，拟合边缘不连续平滑处理整体流程如图4所示。
[0089]
所述步骤s6，基于得到的一个或多个拟合平滑结果，将所述一个或多个拟合平滑结果作为后向投影成像平面进行成像，最终获得ps点分辨单元内部幅度和相位无跳变、边缘不畸变、形状连续的成像结果。
[0090]
本实施例中，复杂场景下高精度dem成像算法的主要流程如图5所示。
[0091]
在本实施例中，根据表1中的参数，以一个450米
×
450米的边坡地形作为成像场景，取其中9个点作为目标点，如图6所示，对其进行仿真，得到的结果如图6所示。
[0092]
表1仿真系统参数
[0093][0094]
图6中，9个目标点的回波ps点分辨单元边缘畸变严重，且存在不连续问题。采用本发明所述方法进行成像，结果如图7所示。图7中，各目标点的分辨单元不再存在边缘畸变问题，且内部也都是连续的。
[0095]
使用北斗meo卫星在11月22日、11月29日、12月6日、12月13日四天的轨道数据进行成像，得到后3天与第1天的3组配准对，将上述两种情况下的成像结果进行ps点提取，分别计算其ps点的相关系数并取均值，所得结果如表2所示。
[0096]
表2仿真结果ps点对相关系数对比
[0097][0098]
表2结果表明，对高精度dem进行拟合处理后，图像对ps点相关性会有所提高，从而提高后续的形变反演精度。
[0099]
下面对实测数据处理结果进行说明。选用马兰庄矿区某侧地形作为成像场景，如图8所示。使用滑动最小二乘拟合结果如图9所示。其中箭头所指区域存在拟合边缘不连续问题，因此使用经过中心对称变换后的高斯平滑核对dem 拟合结果进行平滑处理，处理结果如图10所示。图9和图10对比后可以看到，拟合不连续问题得已解决。
[0100]
图11、图12、图13分别时原始dem、拟合dem、拟合加平滑处理的dem 的成像结果。其中图11的小框共计10个待评估点，用于之后的相关系数计算。将图12和图13的左侧和右侧框分别进行对比可以看到，未平滑处理将带来成像不连续的问题；中间框对比可以看到，处理后的dem成像结果中ps点分辨单元的边缘畸变问题有所缓解。使用北斗meo卫星在11月22日、11月29日、 12月6日、12月13日四天的轨道数据进行成像，得到后3天与第1天的3组配准对，提取原始dem和拟合平滑dem成像结果中所选取的10个待评估ps 点，分别计算这3个配准对的相关系数并取均值，所得结果如表3所示。
[0101][0102][0103]
根据表3结果可以看到，经处理后的dem成像结果中相关系数比处理前有所提升，证明了本发明的有效性。
[0104]
本发明还提供了一种基于导航卫星双基地insar的dem成像装置，该装置包括：
[0105]
拟合模块：配置为使用滑动最小二乘拟合方法拟合本地坐标系下的dem；
[0106]
不连贯单元确定模块：配置为对拟合后的dem的各项进行二次差分，找到所述拟合后的dem存在的所有不连续单元；
[0107]
扩展模块：配置为对每个拟合后的不连续单元，以所述拟合后的不连续单元为中心进行边界扩展，扩展后的结果作为与所述拟合后的不连续单元相对应的待平滑核；
[0108]
高斯平滑核确定模块：配置为对每个待平滑核，生成与所述待平滑核对应的、以所述待平滑核为基准进行中心对称变换后的高斯平滑核；
[0109]
平滑结果生成模块：配置为对每个待平滑核，将其与对应的高斯平滑核进行加权求和运算，进而得到每个待平滑核对应的不连续单元的拟合平滑结果；
[0110]
成像结果确定模块：配置为将所有不连续单元的拟合平滑结果作为最终平滑后的dem，将所述最终平滑后的dem作为后向投影成像平面进行成像，得到最终的成像结果。
[0111]
以上的具体实施例仅描述了本发明的设计原理，该描述中的部件形状，名称可以不同，不受限制。所以，本发明领域的技术人员可以对前述实施例记载的技术方案进行修改或等同替换；而这些修改和替换未脱离本发明创造宗旨和技术方案，均应属于本发明的保护范围。