一种利用卫星立体像对进行网格优化生成DSM的方法

coefficients)参数，并通过构网方法获取粗略的三维网格(以经纬度、高程为坐标)；
12.(2)依据初始的三维网格和线阵卫星立体像对，计算能量函数的光度一致性能量项和梯度平滑能量项，构建能量函数；
13.e(s)＝λ
photoephoto
(s) λ
smoothesmooth
(s)
14.s为三维网格表面；
15.e(s)为三维格网的能量函数；
16.e
photo
(s)为光度一致性能量项；
17.e
smooth
(s)为三维格网梯度平滑的能量项；
18.λ
photo
为光度一致性能量项权重；λ
smooth
为梯度平滑能量项权重；
19.λ
photo
和λ
smooth
的权重值默认为0.5，可根据迭代优化的最终结果来更改权重值，当该值变小时，对应λ
smooth
大，而梯度平滑项控制建立的三维网格的平滑程度，因此会使三维网格更趋于平滑，可控制该权重值来调整最终网格的平滑程度。以zncc(zero-normalized cross-correlation)为主要指标计算光度一致性能量项，以拉普拉斯一阶和二阶梯度为主要指标计算梯度平滑能量项；
20.(3)以自适应阈值对整个三维网格进行细分，使网格细节更丰富，更接近真实地表；
21.以三维网格的三角网投影到二维图像上的区域面积为判断值(实际为投影区域包含的像素和)，默认设置阈值为24，当投影值大于该阈值时，对三角网进行细分，每次迭代都进行细分判断，当达到迭代次数时，细分同时停止。
22.(4)梯度下降法求解，迭代更新，驱动三角网顶点同步移动，得到最优解即优化三角网，生成dsm。
23.本发明提出利用线阵卫星影像直接对粗略的原始格网进行优化，通过rpc模型来转换物点像点间的对应关系，之后通过计算zncc来表征能量函数中的光度一致性能量项，该项用以衡量点由一影像经三维网格重投影至另一影像时，点的偏移程度；同时计算三维网格表面的梯度，该项用以衡量三维网格表面的平滑程度；之后加以正则化约束，使得前两项(光度项和平滑项)平衡；最终对构建的能量函数进行优化求解(梯度下降法调整构成三位网格的三角网顶点位置)。
24.本发明方法对线阵卫星影像构建的粗略网格直接网优化，无需使用摄影测量常用的密集匹配方法。直接对初始粗略网格进行优化，可以从一个近乎空白的网格上优化出接近地表真实情况的三维网格(dsm)；为dsm生产提供了一种新的方法和思路，避免了传统生成方式中密集匹配的窗口的大小适应性问题以及隐含的窗口内视差一致性限制问题，精度更高。
附图说明
25.图1为输入的足够粗略的原始网格。1为某建筑物轮廓，2为原始网格存在的空洞。
26.图2网格优化驱动三角网顶点移动示意图。1为构成格网的三角形的顶点，2为构成网格的基本单位，3为顶点的移动。
27.图3为重投影误差示意图。
28.图4为三维网格细分示意图。1、2、3为三边中点，4、5、6为三角网细分后生成的四个
三角形。
29.图5为三维网格细分后示意图。相较于图2的原始网格，三角网数量更多。
30.图6为图1经网格优化后生成的dsm图。
31.图7为基于网格优化的卫星立体像对dsm生成流程图。
具体实施方式
32.本发明涉及一种利用卫星立体像对进行网格优化生成dsm的方法。本发明基于线阵卫星立体像对的rpc模型，将基于一般图像的三维格网顶点驱动优化问题转换为直接利用线阵卫星立体像对进行三维格网顶点驱动的优化问题，并通过构建能量函数进行优化求解得到优化后三维格网(dsm)。首先根据线阵卫星立体像对生成粗略的原始三维格网；之后依据rpc模型确定物点像点对应关系，计算图像的重投影误差和梯度，以两者表征能量函数的组成项；以自适应阈值对三维网格不断进行细分；采用先离散后优化方法进行优化求解过程。
33.基于以上原理，参照附图来详细说明本发明的实施方式，本实施例是基于网格优化的卫星立体像对dsm的生成方法。
34.1.获取初始网格
35.1.1获取某区域的多景线阵卫星影像，通过摄影测量一般方法构网，得到十分粗糙的三维网格，该网格可足够粗糙，目的是作为网格优化的初始平面。在构建初始三维格网时，可通过摄影测量常用的密集匹配方法(如sgm)或者colmap、visualsfm、cloudcompare等开源项目生成。只需要一个足够粗略的三维网格作为处理平台来进行。
36.本实施例如如图1所示，是通过摄影测量密集匹配算法中的sgm算法利用高分辨率卫星立体像对生成的一个初始平面(足够粗糙)。其中，该初始网格的顶点总数为4443个，构成网格的面个数为8671个。整个网格十分粗糙，依稀可见建筑轮廓，存在明显的空洞情况。
37.1.2初始网格的预处理。对初始的粗略三维网格进行预处理，投影转换等。根据构建的初始网格的坐标系统进行预处理，如输入的初始网格是经纬度高程坐标系统时，就需要对其进行投影转换，转为投影坐标。处理前偏移坐标，处理后复原坐标(偏移时依据rpc参数来进行加减)。其目的是使三维坐标单位一致，便于显示，并保证输入后的精度。
38.本实施例中输入为高分辨率卫星立体像对生成的初始三维网格，先进行x，y坐标偏移，其偏移参数为：-81.65295044719999850713，30.32576204769999961286(根据rpc文件中的latoffset及longoffset参数确定)。初始网格单位为经纬度、高程，通过投影变换，将utm 17n的经纬度转换为utm 17n的投影坐标，使得x，y，高程z单位都变为米。
39.2.构建能量函数
40.2.1构建光度一致性能量项。以zncc为主要指标计算重投影误差，目的是衡量一景影像的点通过三维网格重投影至另一景影像时，其偏离程度。如图3，参考图像(reference image)的顶点经过格网表面s重投影至另一张影像(reprojected image)时，存在重投影误差。由此构建光度一致能量项。其公式如下：
[0041][0042]ephoto
(s)为光度一致性能量项；
[0043]
h(i，j)(x)是衡量像片i和像片j之间光度一致性的递减函数(与归一化互相关zncc相反)；
[0044]
表示像片ij通过网格s重投影至ii上。
[0045]
表示重投影的有效区域；
[0046]
按照上式构建光度一致性能量项，主要以zncc表示，如下计算：
[0047][0048]
其中，f(x，y)是原图像，t(x，y)为模板图像，n是模板中像素(元素)的个数，σ为样本的标准差，μ为样本的均值。
[0049]
在本实施例中，通过计算右影像重投影通过三维网格重投影至左影像的重投影误差，计算得zncc，以zncc表征h(i，j)(x)，构成能量函数的光度一致性能量项。
[0050]
2.2构建梯度平滑能量项。具体实现为计算拉普拉斯梯度，目的是保证三维网格较平滑，惩罚弯曲。其构建公式如下：
[0051][0052]
k1，k2为点所在表面的主曲率。
[0053]
2.3合理分配λ
photo
和λ
smooth
的权重，至此能量函数构建完成。
[0054]
本实施例中，λ
photo
和λ
smooth
的权重分别设置为0.5，0.5，其和为1。λ
photo
和λ
smooth
分别控制三维网格的真实程度和平滑程度，当权重λ
smooth
增大时，最终三维网格棱角更为平滑，其值根据最终的结果来进行调整，不断调整权重值以达到真实性和平滑程度的平衡。
[0055]
3.自适应网格细分
[0056]
依据判断条件对三维网格细分。当影像重投影到另一景影像对应的像素数量超过阈值和纹理复杂程度超过阈值时(投影区域包含的像素数之和以及投影区域的纹理复杂度)，对三角形进行划分。细分阈值的确定，是由构成三维格网的每个三角面与影像的投影关系决定的，计算每个三角面在影像上的投影区域内的像素数及纹理复杂度之和，设置为阈值，大于该阈值则进行细分。划分原理如图4所示，取每个边中点作为新的顶点，将原三角形划分为新的四个三角形，目的是实现网格的细分，使细节更加完善。
[0057]
本实施例中，如图5所示，细分阈值的参数设置为24，该阈值可手动设置。随参数设置的大小，最终网格的精细程度也不同，阈值越大，则网格细节越粗糙；反之，网格细节越丰富。
[0058]
4.能量函数最优化
[0059]
4.1能量函数最优化求解。迭代初始值为初始的三维网格，实际使用梯度下降法迭代求得最优解，每次迭代计算三角形各顶点的偏移量，不断迭代直至三角网顶点移动至接近地表真实情况的最优位置。如图5示，整个三维网格由成千上万个三角网构成，在进行梯度下降法求解时，初始梯度为网格中所有顶点的集合，每个顶点都由x，y，z三方向的坐标构成。求解过程中，x，y，z三方向的值迭代变化，使得顶点坐标不断移动。当达到设置的迭代次数时停止，停止时结果即为最优化结果。其顶点运动示意图如图2。
[0060]
本实施例中，迭代参数设置为255.05，步长设置为1.5，移动后优化结果如图6。与原始输入网格图1相比，细节层次明显更加丰富，建筑物轮廓明显可见，原始网格存在的空
洞也填补完毕。
[0061]
4.2限制步长，同步移动。迭代求解过程中，为避免三角网的某个顶点移动量过大，必须对每个顶点的移动量加以限制。取三维网格平均边长为阈值，对驱动顶点的移动量进行约束，使每个顶点同步移动。计算构成整个三维网格的每个三角网的边长长度作为限制步长的判断依据，当x，y，z方向上的移动量大于每个三角网的平均边长时，将x，y，z方向上的移动量(梯度变化)限制为三维网格平均边长的一半。这样使得三角网的三个顶点同步移动，避免出现某个顶点移动量过大而使得梯度下降法求解结果不佳的情况。
[0062]
本实施例中，当某个顶点的移动量过大时，在迭代次数一定的情况下，可能出现无法优化到最佳位置的情况，如图5所示，圈定区域的凸起，考虑原因可能是某几个三角网的顶点移动量相对于其他顶点过大，因此出现“凸起”、“凹陷”情况。因此需要对顶点的移动加以限制。