基于星载双站InSAR和激光雷达数据估计DEM的方法、装置、设备及介质

本发明属于合成孔径雷达干涉测量(insar)的大地测量领域，具体涉及一种基于星载双站insar和激光雷达数据估计dem的方法、装置、设备及介质。

背景技术：

1、数字高程模型(digital elevation model,dem)是自然资源调查、灾害监测以及气候和生态变化分析不可或缺的数据源，同时也是经济建设和军事国防必不可少的基础图件。合成孔径雷达干涉测量(interferometric synthetic aperture radar，insar)技术，具备全天时、全天候和大范围等优势，特别是双站insar技术，不受时间去相干和大气延迟的影响，是一种获取和更新大范围dem的强有力工具。

2、目前，常用的insar技术估计dem的方法主要包括以下两种：1)基于差分干涉技术(d-insar)的dem估计，原理为使用insar干涉相位与现有的dem产品进行差分，然后对差分相位进行解缠和相高转换后加回到现有dem上生成新的dem，该方法简便易实现，但严重依赖于现有的外部dem产品，且在大范围估计时易将现有dem的趋势性误差引入到新的dem中；2)基于原始insar技术的dem估计，该类方法估计的dem精度较高，但是需要依赖大量的地面控制数据(如角反射器位置、gnss cors站位置等)对轨道进行精炼以及进行绝对相位校正，然后再进行相高转换估计dem。此外，以上两类方法都通常需要使用外部dem辅助进行地理编码。

3、综上所述，现有的insar dem估计方法严重依赖于现有地形产品和地面控制数据，如何仅依赖星载平台估计dem在军事国防中具有重大意义。因此，有必要设计一种不依赖现有dem产品和地面数据，仅基于星载平台的大范围、高精度dem的估计方法。

技术实现思路

1、本发明提供一种基于星载双站insar和激光雷达数据估计dem的方法、装置、设备及介质，可以在不依赖于现有的外部dem产品或地面控制数据的条件下，实现大范围、高精度的dem估计。

2、为实现上述技术目的，本发明采用如下技术方案：

3、一种基于星载双站insar和激光雷达数据估计dem的方法，包括以下步骤：

4、步骤1，获取待研究区域的双站insar影像数据，并进行干涉处理，获取insar数据的相干性和由以下几个部分组成的干涉相位：

5、φint＝φflat+φtop+φorb+φnoise# (1)

6、式中φflat表示平地相位，φtop是地形相位，φorb是轨道误差相位，φnoise是噪声相位；

7、步骤2，对步骤1得到的干涉相位，首先根据sar卫星的轨道参数进行去平地相位处理，然后进行滤波和解缠处理得到解缠后的相对相位；

8、步骤3，获取待研究区域的星载激光雷达的地面测高点作为控制点，根据距离-多普勒模型将其地理编码到sar影像上对应的像素坐标，并使用匹配窗口的斜率相关性分析进行精化校正；

9、步骤4，使用星载激光雷达控制点对insar基线参数进行精炼并计算相位偏移获取解缠后的绝对相位，再通过相高转换得到像素坐标下的数字高程模型dem，以及利用距离-多普勒模型将dem从像素坐标转换到地理坐标；

10、步骤5，以星载激光雷达为控制点，根据匹配窗口的斜率相关性分析均匀获取dem不同位置的地理坐标平面偏移量，然后使用多项式拟合偏移量对整景dem的平面位置进行校正；

11、步骤6，构建考虑轨道误差和因绝对相位偏移校正不准确导致的系统偏差的高程误差模型，以星载激光雷达地面控制点提供绝对的高程约束，选用多个不同轨道获取的dem提供相对高程约束，基于步骤5平面位置校正后的dem和高程误差模型进行区域网平差计算以获取高程误差模型系数，然后将高程误差从初始dem中移除获取高程校正后dem；

12、步骤7，循环迭代运行步骤5和步骤6，并以激光雷达地面控制点为参考计算均方根误差，当均方根误差变化量小于给定阈值，则迭代结束，并输出最终估计的dem。

13、进一步地，去平地相位时，平地相位的计算公式为：

14、

15、式中，λ为雷达信号波长，b是基线长度，θ是入射角，α为基线与水平方向的夹角，m在双站insar系统中为1，在单站insar系统中为2。

16、进一步地，所述距离-多普勒模型表示为：

17、

18、式中，(xs,ys,zs)表示sar卫星在t时刻地理坐标系中的位置，(vxs,vys,vzs)表示sar卫星t时刻的速度，(xp,yp,zp)表示星载激光控制点的地理位置，r表示斜距，a和b表示参考椭球的长短半轴，h为激光雷达控制点的高程；

19、所述根据距离-多普勒模型将其地理编码到sar影像上对应的像素坐标，具体为：

20、基于距离-多普勒模型中的r、a、b和h为已知值，卫星的位置(xs,ys,zs)和速度(vxs,vys,vzs)为待估值，通过牛顿迭代法不断迭代sar卫星的轨道时间t；当时间t所对应的卫星位置(xs,ys,zs)和速度(vxs,vys,vzs)满足公式(3)时，通过以下公式得到激光雷达控制点对应的sar像素坐标：

21、

22、式中，r为控制点在sar像素坐标的行号，c为控制点在sar像素坐标的列号，t为上述迭代所得到的时间，t0为sar影像的成像开始时间，δt为sar成像的时间间隔，r为公式(3)中所用斜距，r0为近距端斜距，δr为距离向分辨率。

23、进一步地，步骤3中使用匹配窗口的斜率相关性分析对每个控制点的像素坐标进行精化校正，具体包括：

24、首先，在像素坐标下以待校正的控制点为中心选取预设大小的匹配窗口；设匹配窗口内包括n个离散的星载激光雷达控制点，计算匹配窗口内每两个控制点之间的斜率：

25、

26、式中，hi和hj表示匹配窗口内编号为i和j的两个控制点的高程，(xi,yi)和(xj,yj)表示两个控制点的像素坐标；

27、然后，将步骤2得到的相对相位转换为基于insar的相对高程，并建立以待校正控制点为中心的搜索窗口；其中，搜素窗口的尺寸大于匹配窗口；在搜索窗口内滑动匹配窗口，并计算基于insar相对高程的斜率进而计算匹配窗口内星载激光雷达控制点斜率k和基于insar的相对高程的斜率的相关性；并选取最大相关性且大于设定阈值的位置，计算该位置相对于匹配窗口中心的行和列偏移量；

28、最后，将偏移量与控制点在校正前的像素坐标相加，即得到控制点校正后的像素坐标。

29、进一步地，步骤5中使用多项式拟合偏移量对整景dem的平面位置进行校正的表达式为：

30、

31、式中，aj和bj为多项式拟合系数，j＝0,1,2,3,4,5；x和y为dem在地理坐标系下的行列号，δx和δy为dem相对于星载激光雷达控制点的偏移量；

32、通过解算多项式拟合系数，将地理坐标系下的dem重采样到新的地理位置，从而实现平面位置校正。

33、进一步地，步骤6的具体过程包括：

34、步骤6.1，构建考虑轨道误差和因绝对相位偏移校正不准确导致的系统偏差的高程误差模型，表示为：

35、g(x,y)＝a0+a1x+a2y+a3xy+a4x2+a5x3# (7)

36、式中，ak表示轨道误差和相位偏移校正的系统偏差的多项式系数，k＝0,1,2,3,4,5；x和y为dem在sar影像像素坐标系下的行列号；

37、步骤6.2，对于每一景insar估计的dem，都使用步骤5中的星载激光雷达控制点为区域网平差提供绝对高程约束：

38、hi＝hi+vi# (8)

39、式中，hi表示第i个星载激光雷达控制点的高程；hi表示第i个控制点对应的初始dem的高程值；vi表示第i个控制点的高程误差值，使用高程误差模型表示，即vi＝g(xi,yi)，xi和yi为第i个控制点在sar影像像素坐标系下的行列号；假设每景影像有m个控制点，对应的矩阵表达式为：

40、

41、步骤6.3，在相邻的两景insardem的重叠区域内，均匀地选取控制点作为连接点；按照不同insardem在同一地理位置高程大小相等原理，为区域网平差提供相对高程约束，任一连接点在不同的insardem上满足条件：

42、hi,j+gj(xi,j,yi,j)＝hi,k+gk(xi,k,yi,k)# (10)

43、式中，j和k用于区分相邻的两景insar；hi,j表示第i个连接点在demj中的高程值，hi,k表示第i个连接点在demk中的高程值，gj(xi,j,yi,j)表示第i个连接点在demj中的高程误差，gk(xi,k,yi,k)表示第i个连接点在demk中的高程误差；对于n个连接点，矩阵表达式为：

44、hi,j-hi,k＝-gj(xi,j,yi,j)+gk(xi,k,yi,k)# (11)

45、

46、步骤6.4，根据步骤6.2中建立的m个控制点的绝对高程约束方程和步骤6.3中建立的n个连接点的相对高程约束方程，联立建立区域网平差的方程：

47、v＝bx-l# (13)

48、式中，v是残余误差，b是整体系数矩阵，x是待估参数矩阵，l是高程误差向量；对于j和k两景sar影像，矩阵b、x、l分别表示为：

49、

50、l＝[[h1,j … hm,j]-[h1,j … hm,j][h1,j … un,j]-[h1,k … hn,k][h1,k … hm,k]-[h1,k … hm,k]]t

51、在此基础上，使用最小二乘估计对待估参数x进行解算：

52、x＝(btpb)-1(btpl)# (14)

53、式中，p表示权重矩阵，其值为对角线元素为1的单位矩阵，表示所有控制点和连接点的权重相等；

54、步骤6.5，利用估计的模型系数矩阵x计算每景insar dem的高程误差并将其从原始dem中移除。

55、进一步地，步骤7循环迭代运行步骤5和步骤6的停止标准为：以研究区域的所有星载激光雷达数据为参考，计算当前得到的dem的均方根误差与上一次运行得到的dem的均方根误差的差异，当均方根误差差异的绝对值小于给定阈值时停止迭代并输出最终的dem结果。

56、一种基于星载双站insar和激光雷达数据估计dem的反演装置，包括：

57、预处理模块，用于：对获取的待研究森林区域的insar数据进行干涉处理，得到影像的相干性和相位；并对相位进行平地相位去除、滤波和解缠处理，得到未缠绕的相对相位；

58、地理编码模块，用于：将获取的星载激光雷达编码到sar影像像素坐标和将insardem编码到地理坐标，并且利用斜率相关性分析对初始编码结果进行精化校正；

59、初始dem获取模块，用于：使用星载激光雷达作为控制点对insar基线参数进行精炼并计算相位偏移获取解缠后的绝对相位，并进行相高转换获取初始insar dem；

60、区域网平差计算模块，用于：构建考虑轨道误差和因绝对相位偏移校正不准确导致的系统偏差的高程误差模型，并选用裸地区域的星载激光控制点为区域网平差提供绝对高程约束、选用连接点建立相对高程约束；基于初始的dem和高程误差模型进行平差运算，得到高程误差模型的系数，并将高程误差从初始dem中移除获取高程校正后的dem；

61、迭代模块，用于：再次调用地理编码模块进行平面位置校正和调用区域网平差计算模块进行高程校正，计算均方根误差差异判断是否满足要求和是否终止迭代，并得到最终估计的dem结果。

62、一种电子设备，包括存储器及处理器；其中，所述存储器用于存储计算机程序指令，所述处理器用于执行所述存储器的计算机程序指令，具体执行上述任一技术方案所述的方法。

63、一种计算机可读存储介质，其上存储有计算机程序，所述计算机程序被处理器执行时实现上述任一项技术方案所述方法。

64、有益效果

65、本发明基于星载双站insar和激光雷达数据估计dem的方法，实现了仅基于星载平台数据估计大范围的dem：首先利用常用常规的数据处理方法对双站星载双站insar数据进行处理得到解缠相位；之后将星载激光雷达数据编码到sar影像像素坐标系下作为控制点，对insar基线参数进行精化和对解缠相位进行绝对校正，利用常规的相高转换方法估计得到初始dem并将其编码到地理坐标系；然后使用地理坐标系下的星载激光雷达作为控制点，提出了一种基于斜率相关性分析的方法对其平面位置进行校正；最后，建立区域网平差模型对初始dem进行高程校正，并进行迭代估计最终的dem结果。

66、本方法的有益效果为：构建了仅基于星载平台数据估计大范围的dem的方法体系，克服了现有的insar dem估计方法严重依赖于现有地形产品和地面控制数据的难题；提出了一种基于斜率相关性分析的星载激光雷达数据与dem配准的方法；建立了一种考虑轨道误差和因绝对相位偏移校正不准确导致的系统偏差的高程误差模型，实现了dem高程误差的校正。在精度上，本发明在不依赖于现有地形产品和地面控制数据的前提下，达到了与现有方法相当的精度水平，是一种自主可控、大范围、高精度dem估计的有效方法。