CN113222871B - 一种多视卫星影像数字表面模型融合方法 - Google Patents

Abstract

本发明属于多视卫星影像三维重建技术领域，具体涉及一种多视卫星影像数字表面模型融合方法。本发采用先分割后聚类的思想，先通过分割对高程落差严重的地物边缘进行定位，然后通过表面模型聚类为每个区域内DSM融合提供模型约束。本发明同时利用颜色信息和高程信息实现相同区域的聚集，克服因高曝光导致的边缘弱化问题，实现对地物边缘的准确定位；通过超像素分割方法定位DSM的地物边缘，为多视DSM融合提供准确的边缘依据，使得区域内的多视DSM融合能够沿着区域边缘进行。本发明获得的融合后的DSM能够解决因大气辐射变化和强曝光导致的地物边缘区域高程模糊、高程落差划分不规则的问题。

Description

一种多视卫星影像数字表面模型融合方法

技术领域

本发明属于多视卫星影像三维重建技术领域，具体涉及一种多视卫星影像数字表面模型融合方法。

背景技术

通过卫星影像进行三维重建能够突破机载航空影像三维重建的地域限制，在全球三维态势感知和大规模地面信息判读等任务中发挥关键作用。因此，在遥感工程、摄影测量和计算机视觉的研究领域，通过多视卫星影像获取高精度三维信息一直都是重点研究内容。

多视图三维重建是对两视图三维重建的扩展和延续。根据三维重建中密集匹配执行的顺序，可以将多视三维重建方法划分为多视匹配方法(multi-view matching，MVM)和多立体匹配方法(multi-stereo matching，MSM)两类。MVM方法需要同时考虑多个视图间的匹配关系，并利用光度一致性和可见性约束在不同的视图之间建立联系来过滤误匹配，因此MVM更善于处理从地面移动平台和机载平台拍摄的多视透视影像，这些影像具有大基线和大视差的特点。在三维模型生成阶段，MVM通过多种融合方法生成完整的点云图或数字表面模型(digital surface model，DSM)，包括基于体素的方法、基于表面演化的方法和基于插值增长的方法。在这些常用的方法中，基于体素的方法仅适用于紧密封闭的空间内的小型紧凑物体；基于表面演化的方法需要可靠的初始猜测，然而对于大规模场景则很难获得可靠的初始点；基于插值生长的方法需要预生成准确的三维起始点作为生长的指导，以生成准确的表面模型。但是，对于远离起始点的区域则可能会产生不合实际的插值。上述这些MVM方法都需要通过光束法区域网平差获得准确的高程，光束法区域网平差依赖数量充足的、准确的视图间匹配点作为地面控制点。然而对于大篇幅卫星影像来说，获得数量充足的准确地面控制点则需要耗费大量的人力成本。

MSM对两视图三维重建的直接扩展，这些方法通过融合每一对两视三维重建结果来获得完整、准确的DSM。得益于MSM方法能够直接使用最先进的密集匹配方法，因此MSM更加适用于大场景可见光遥感影像的三维重建任务。就准确性和完整性而言，航拍影像三维重建结果通常高于卫星影像三维重建结果。这是因为多视立体卫星影像的获取很容易受到卫星运行机制、采集时间间隔和大气辐射等方面的影响。最常见的MSM的优化方法是通过边缘检测或者区域分割方法对密集匹配进行边缘约束，以获得边缘更加准确的视差图，这样能通过两视三维重建计算出边缘精度更高的DSM。最后直接对多幅DSM执行融合获得最终的结果。但在多视DSM融合的过程中，多幅DSM之间并没有边缘约束，即使两视图重建结果被优化过，最终获得多视DSM在高程落差严重的区域依旧存在大量误差，并且弱纹理区域、阴影区域的高程空洞也没有得到有效的填充。

发明内容

本发明的目的在于提供一种多视卫星影像数字表面模型融合方法。

本发明的目的通过如下技术方案来实现：包括以下步骤：

步骤1：卫星影像与DSM的配准，通过卫星影像配置文件中的RPC系数构建RFM，建立图上坐标与大地坐标的直接映射关系，并根据大地坐标获得图上目标区域对应的DSM，在目标区域中实现卫星影像与DSM的配准；

步骤2：卫星影像及其对应的DSM联合数据超像素分割；将LSC超像素使用的十维映射函数扩展为十一维映射函数，然后对联合数据进行映射，并对十一维的像素进行利用谱聚类，获得卫星影像和DSM的联合数据的超像素分割结果，实现具有相同高程结构区域的聚集；

步骤3：在每个超像素区域内，利用DSM的3-D坐标进行T-Linkage表面模型聚类，在超像素区域内获得集群标签及其对应的表面模型；

步骤3.1：通过均匀采样与局部采样相结合的随机抽样获取估计表面模型所需的最小数据3-D点集，生成m个表面模型H＝{h_j}_j＝1...m；

步骤3.2：拟合表面模型h；

通过奇异值分解SVD的方法获得集群3-D空间中的表面模型ax+by+cz+d＝0；根据模型性质，系数a,b,c构成的向量N＝[a,b,c]^T为表面模型的法向量；首先求最小样本集中三维点{(x_i,y_i,z_i)}_i＝1...n的质心

n≥3，根据法向量的性质构建方程组：

令方程组左边矩阵为A，右边矩阵为法向量N，则拟合的目标函数为min||AN||，并且保证||N||＝1，则对A做奇异值分解A＝UDV^T；D为对角阵，U和V均为酉矩阵，根据酉矩阵的性质能得到：

||AN||＝||UDV^TN||＝||DV^TN||

当满足条件||V^TN||＝||N||＝1时，D的对角元素为奇异值，最后一个对角元素为最小奇异值；当且仅当V^TN＝[0 0 1]^T时，||AN||能取到最小值，最小奇异值对应的奇异向量就是平面的法方向，也就求出了系数a,b,c：最终可以得出

步骤3.3：在每个样本集中计算数据点x的偏好函数，每个数据点x都用偏好函数φ来描述；偏好函数在取值区间为闭合区间[0,1]，最喜欢被表示成1，最不喜欢被表示成0；通过阈值的截断效果来计算数据x的偏好函数；

τ为时间常数能够实现截断效果；d(x,h)表示数据点x与表面模型h之间的残差；

步骤3.4：通过数据点x与表面模型h之间的残差公式d(x,h)计算点到平面之间的距离；

d(x,h)＝|ax_i+by_i+cz_i+d|

步骤3.5：如果数据子集S＝{x}，则S的偏好函数就是x的偏好函数表示为q_x；否则S的偏好函数是向量p∈[0,1]^m，使其第i个分量为x的最小偏好函数

然后通过Tanimoto距离计算每个点偏好函数[0,1]^m之间的相似性：

其中，符号<·,·>表示

的标准内积；Tanimoto距离范围是[0,1]，对于

有如果p＝q，则d_T(p,q)＝0；如果p⊥q，则d_T(p,q)＝1；

步骤3.6：如果任意p和q都正交，那么输出集群标签及其对应的表面模型，聚类结束；否则通过公式d_T(p,q)＝min_r,sd_T(r,s)找到p和q，并用二者的并集替换这两个集群，重新标记新集群，返回步骤3.2；

步骤4：将T-Linkage聚类获得的N个集群作为DSM高程值融合的基本单位，利用表面模型的约束作用在每个集群区域中执行DSM的融合以获得最终的DSM；

步骤4.1：初始化集群标签i＝1，融合结果

步骤4.2：对于第i个集群的3-D点集合

中的点p，若p到虚拟表面h_i的距离小于阈值τ_i，则将像素点p收入符合模型的点集合

步骤4.3：通过中值滤波的方法计算最终的DSM，将点集合

作为中值滤波的作用范围，中值滤波计算如下：

其中，将

作为第i个集群的高程值，得到第i个集群的DSM：

步骤4.4：合并所有集群中3-D点的中值滤波结果，获得最终融合的DSM：

步骤4.5：令i＝i+1，如果当前集群标签i＝N+1，那么输出融合结果

否则返回步骤4.2。

本发明还可以包括：

所述的步骤1中卫星影像与DSM的配准的方法具体为：

步骤1.1：通过RFM建立卫星影像坐标与大地坐标的直接映射关系；

获取构建RFM的80个RPC系数，通过这80个RPC系数构建图上坐标(l,s)与大地坐标

之间的关系；令(B_n,L_n,H_n)为一个点的正则化的经度、纬度和高程坐标；

L_n＝(λ-L_OFF)/L_SCALE

H_n＝(h-H_OFF)/H_SCALE

其中，B_OFF、L_OFF和H_OFF为卫星配置文件中的正则化偏移参数；B_SCALE、L_SCALE和H_SCALE为卫星配置文件中的正则化比例参数；

有理函数的多项式方程为：

其中，L_num(·)、L_den(·)、S_num(·)和S_den(·)是三阶多项式函数并且拥有20个RPC系数；

步骤1.2：根据RFM将图像上坐标转换为经纬度坐标，再根据经纬度坐标在DSM上获取对应目标区域的DSM，实现卫星影像与多幅DSM之间的配准。

所述的步骤2中卫星影像及其对应的DSM联合数据超像素分割的方法具体为：

步骤2.1：将卫星影像的颜色转化为CIELAB颜色空间中的颜色分量，构建每个像素点p的六维坐标p＝(l_p，α_p，β_p，x_p，y_p，e_p)；其中l_p、α_p和β_p是像素点p在CIELAB颜色空间中的颜色分量值；x_p和y_p是图像平面中像素点p的水平和垂直坐标；e_p是像素点p的高程值；

步骤2.2：根据映射函数φ_e将每个六维像素点p映射到十一维；

其中，ω为权重；C_c和C_s用来控制颜色和空间信息的相对重要性，C_s/C_c的值越大空间上越接近的像素越倾向于聚成一类；C_e用来平衡感官信息与高程信息之间的重要性，C_e/C_c的值越大高程相近的像素越倾向于聚成一类；

步骤2.3：初始化加权K-means聚类；赋予每个像素点的权重ω(p)，令K表示聚类数，并在整个图像上分别以水平间隔v_x和垂直间隔v_y对K个种子像素进行均匀采样；π_k表示第k个聚类，k＝1,2,...,K，c_k是该类的聚类中心；初始化每个像素p的标签L(p)＝0，初始化每个像素p的距离d(p)＝∞；

步骤2.4：加权K-means聚类，通过欧几里得距离来计算每个像素点φ_e(p)与聚类中心c_k之间的距离d_e(p)；

其中，聚类中心c_k的计算如下式：

步骤2.5：如果d_e(p)＜d(p)，则更新距离d(p)＝d_e(p)和标签L(p)＝k；

步骤2.6：更新所有的像素的标签和聚类中心，如果聚类中心不再发生变化，出现收敛，则输出当前聚类标记结果，否则返回步骤2.4。

本发明的有益效果在于：

本发明针对传统多视深度图融合方法应用在多视卫星影像DSM融合上出现的地物边缘误差大、细节丢失以及高程漏洞多的问题，采用先分割后聚类的思想，先通过分割对高程落差严重的地物边缘进行定位，然后通过表面模型聚类为每个区域内DSM融合提供模型约束。本发明同时利用颜色信息和高程信息能够实现相同区域的聚集，并且能够克服因高曝光导致的边缘弱化问题，实现对地物边缘的准确定位。本发明通过超像素分割方法定位DSM的地物边缘，为多视DSM融合提供准确的边缘依据，使得区域内的多视DSM融合能够沿着区域边缘进行。本发明获得的融合后的DSM能够解决因大气辐射变化和强曝光导致的地物边缘区域高程模糊、高程落差划分不规则的问题。本发明在每个分割区域中都执行表面模型聚类，实现相同高程结构的3-D点聚集及其对应的模型拟合，既能够通过聚类挖掘出“隐藏”的表面细节，恢复地物的精细表面结构，又能为因强曝光和建筑物阴影导致的高程漏洞区域提供表面模型的约束，进而提高DSM的精细度和完整度。

附图说明

图1是本发明的流程图。

图2是本发明中卫星影像及其对应的DSM的联合数据超像素分割结果图。

图3是本发明中表面模型聚类过程图。

图4是本发明中多视DSM三维点融合过程图。

图5是本发明获得的DSM对比结果图。

具体实施方式

下面结合附图对本发明做进一步描述。

本发明提供的是一种多视卫星影像数字表面模型(digital surface model，DSM)融合方法，它属于多视卫星影像三维重建技术领域。本发明解决了现有多视卫星影像三维重建方法在高程落差严重的建筑物边缘的高程误差、在平滑的地面目标表面出现的高程漏洞以及地面目标表面细节丢失的问题。

本发明首先通过卫星影像配置文件中的有理多项式系数(rational polynomialcoefficient，RPC)构建有理函数模型(rational function model，RFM)建立图上坐标与大地坐标的直接映射关系，并根据大地坐标获得目标区域的DSM，在卫星影像与DSM之间实现配准。然后实现卫星影像及其对应的DSM联合数据的超像素分割。在线性谱聚类(LinearSpectral Clustering，LSC)超像素中添加高程信息，将原线性谱聚类LSC超像素方法的映射函数的维数扩展到十一维，将其应用于卫星影像和DSM的联合数据的超像素分割。之后在每个超像素区域内利用DSM的三维点坐标进行T-Linkage表面模型聚类，在超像素区域内获得集群标签及其对应的表面模型。最后以聚类获得的集群为基本单位，依次在其覆盖范围内执行基于表面模型的约束的DSM融合，获得最终的DSM融合结果。本发明的流程图如图1所示。与现有的方法相比，本发明能够根据分割结果准确地定位到高程落差严重的地物边缘，使该处的多视DSM融合始终沿着分割边缘进行，因此获得的DSM能够准确划分高程落差处的边界，实现不同高度区域的准确划分。另外本发明通过模型聚类方法对不同高程区域进行聚类，使表面平滑的区域、强曝光区域和阴影区域也能在聚类获得的表面模型的约束下完成高程数据的准确填充。本发明能够应用于多视卫星影像三维重建领域。

本发明是这样实现的：

步骤一、卫星影像与DSM的配准。通过卫星影像配置文件中的RPC系数构建RFM，建立图上坐标与大地坐标的直接映射关系，并根据大地坐标获得图上目标区域对应的DSM，在目标区域中实现卫星影像与DSM的配准。

步骤二、卫星影像及其对应的DSM联合数据超像素分割。提出新的应用于2-D卫星影像与3-D DSM联合数据的超像素方法。将LSC超像素使用的十维映射函数扩展为十一维映射函数，然后对联合数据进行映射，并对十一维的像素进行利用谱聚类，获得卫星影像和DSM的联合数据的超像素分割结果，实现具有相同高程结构区域的聚集。

步骤三、T-Linkage表面模型聚类。在每个超像素区域内，利用DSM的3-D坐标进行T-Linkage表面模型聚类，在超像素区域内获得集群标签及其对应的表面模型。

步骤四、多视DSM融合。将T-Linkage表面模型聚类获得的集群作为融合的基本单位，在其覆盖范围内执行基于表面模型的约束的DSM融合。

本发明的有益效果是：

本发明提出了一种多视卫星影像数字表面模型融合方法，针对传统多视深度图融合方法应用在多视卫星影像DSM融合上出现的地物边缘误差大、细节丢失以及高程漏洞多的问题，本发明采用先分割后聚类的思想，先通过分割对高程落差严重的地物边缘进行定位，然后通过表面模型聚类为每个区域内DSM融合提供模型约束。本发明具有以下优点：

第一，本发明提出一种适用于2-D卫星影像与3-D DSM联合数据的超像素分割方法获得联合数据的超像素结果。该方法同时利用颜色信息和高程信息能够实现相同区域的聚集，并且能够克服因高曝光导致的边缘弱化问题，实现对地物边缘的准确定位。因此该方法能够作为基础方法作用于基于DSM的目标识别、高程信息判读等应用。另外影像与DSM联合数据是一种包含颜色信息和高度信息的数据，与包含颜色信息和深度信息的RGB-D数据的数据组成一样。这种超像素分割方法能够用于RGB-D数据的超像素分割。

第二，本发明通过超像素分割方法定位DSM的地物边缘，为多视DSM融合提供准确的边缘依据，使得区域内的多视DSM融合能够沿着区域边缘进行。因此本发明获得的融合后的DSM能够解决因大气辐射变化和强曝光导致的地物边缘区域高程模糊、高程落差划分不规则的问题。

第三，本发明在每个分割区域中都执行表面模型聚类，实现相同高程结构的3-D点聚集及其对应的模型拟合。这样既能够通过聚类挖掘出“隐藏”的表面细节，恢复地物的精细表面结构，又能为因强曝光和建筑物阴影导致的高程漏洞区域提供表面模型的约束，进而提高DSM的精细度和完整度。

实施例1：

本发明直接对多视DSM进行融合获得准确性和完整性更高的DSM。通过对卫星影像及其对应的DSM的联合数据进行区域分割，获得DSM的区域分割结果，这样在高程落差严重的区域就有了边缘约束。另外，为了解决超像素分割方法存在的表面细节区域丢失的问题，本发明通过基于表面模型的3-D点聚类的方式检索区域中的表面细节，在每个区域中执行T-Linkage表面模型聚类获得隐藏在区域内的细节区域，包括局部凸起、凹陷等区域，这样既能在表面模型的约束下填充平滑区域的高程漏洞，也能利用聚类对物体表面的其他细节区域进行发掘和保护，提高了融合DSM的完整度和准确度。

图2是本发明中卫星影像及其对应的DSM的联合数据超像素分割结果。其中(a)是卫星影像的超像素分割结果，(b)是DSM的超像素分割结果。红色线条表示超像素边缘。

图3是本发明中表面模型聚类过程。图3中(a)为部分卫星影像；(b)为对应于区域(a)的真实LiDAR点云；(c)为T-Linkage的表面模型聚类结果，不同的颜色表示不同的类别；(d)为表面模型结果。(c)与(d)的点云对应于(a)中的红色区域和(b)中的蓝色方框。黑色虚线箭头表示点云模型之间的对应关系。

图4是本发明中多视DSM三维点融合过程。其中左图表示通过T-Linkage获得的其中一个类别，并计算出属于该类别的表面模型，而且能够根据点到表面模型之间的距离d筛选出类内点(绿点)和离群点(红点)；右图表示根据获得的类内点(绿点)计算中值，并将其作用于整个聚类区域。

图5是本发明获得的DSM对比结果。其中，(a)是卫星影像；(b)对应于(a)中红色框的真实LiDAR点云；(c)对应于(a)中红色框的融合前的DSM；(d)对应于(a)中红色框的本发明获得的融合后的DSM。

步骤一、卫星影像与DSM配准

本步骤在多幅卫星影像与多幅DSM之间实现配准，即获取能够用于融合的公共区域，本步骤作为多视DSM的基础步骤，详细步骤如下：

Step1：通过RFM建立卫星影像坐标与大地坐标的直接映射关系。卫星影像供应商都会提供构建RFM的80个RPC系数。通过这80个RPC系数构建图上坐标(l,s)与大地坐标

之间的关系。令(B_n,L_n,H_n)为一个点的正则化的经度、纬度和高程坐标。

其中B_OFF，L_OFF和H_OFF为卫星配置文件中的正则化偏移参数，B_SCALE，L_SCALE和H_SCALE为卫星配置文件中的正则化比例参数，那么有理函数的多项式方程为：

其中L_num(·)，L_den(·)，S_num(·)和S_den(·)是三阶多项式函数并且拥有20个RPC系数。

Step2：根据RFM将图像上坐标转换为经纬度坐标，再根据经纬度坐标在DSM上获取对应目标区域的DSM，实现卫星影像与多幅DSM之间的配准。

步骤二、卫星影像与DSM的联合数据超像素分割

本步骤利用配准的卫星影像与DSM的联合数据实现DSM的超像素分割，目的是聚集具有相同颜色和高程结构的区域和定位高程落差严重的目标边缘，分割结果如图2所示，详细步骤如下：

Step1：将卫星影像的颜色转化为CIELAB颜色空间中的颜色分量，构建每个像素点p的六维坐标p＝(l_p，α_p,β_p,x_p,y_p,e_p)，其中l_p，α_p和β_p是像素点p在CIELAB颜色空间中的颜色分量值，x_p和y_p是图像平面中像素点p的水平和垂直坐标，e_p是像素点p的高程值。

Step2：根据映射函数φ_e将每个六维像素点p映射到十一维，

其中，ω为权重，C_c和C_s被用来控制颜色和空间信息的相对重要性，C_s/C_c的值越大空间上越接近的像素越倾向于聚成一类。C_e被用来平衡感官信息与高程信息之间的重要性，C_e/C_c的值越大高程相近的像素越倾向于聚成一类。

Step3：初始化加权K-means聚类。赋予每个像素点的权重ω(p)，令K表示聚类数，并在整个图像上分别以水平间隔v_x和垂直间隔v_y对K个种子像素进行均匀采样。π_k表示第k个聚类(k＝1,2,…,K)，c_k是该类的聚类中心。初始化每个像素p的标签L(p)＝0，初始化每个像素p的距离d(p)＝∞。

Step4：加权K-means聚类。通过欧几里得距离来计算每个像素点φe_e(p)与聚类中心c_k之间的距离d_e(p)，如下：

其中聚类中心c_k的计算如下式：

Step5：如果d_e(p)＜d(p)，那么更新距离d(p)＝d_e(p)和标签L(p)＝k；

Step6：更新所有的像素的标签和聚类中心，如果聚类中心不再发生变化，出现收敛，那么返回当前聚类标记结果，否则返回Step4；

步骤三、T-Linkage表面模型聚类

在获得区域分割结果之后，在每一个超像素区域中执行T-Linkage聚类以获得每个超像素区域中的一个或多个3-D点集群及其对应的表面模型，如图3所示，详细步骤如下：

Step1：生成最小样本集。通过均匀采样与局部采样相结合的随机抽样获取估计表面模型所需的最小数据3-D点集(称为最小样本集)，然后生成m个表面模型如下:

H＝{h_j}_j＝1...m (6)

Step2：拟合表面模型h。通过奇异值分解(SVD)的方法获得集群3-D空间中的表面模型ax+by+cz+d＝0。根据模型性质，系数a，b,c构成的向量N＝[a,b,c]^T为表面模型的法向量。首先求最小样本集中三维点{(x_i,y_i,z_i)}_i＝1...n(n≥3)的质心

根据法向量的性质构建方程组：

令左边矩阵为A，右边矩阵为法向量N，则拟合的目标函数为min||AN||，并且保证||N||＝1，则对A做奇异值分解A＝UDV^T。那么D为对角阵，U和V均为酉矩阵，根据酉矩阵的性质能得到：

||AN||＝||UDV^TN||＝||DV^TN|| (8)

当满足条件||VTN||＝||N||＝1时，D的对角元素为奇异值，最后一个对角元素为最小奇异值，当且仅当V^TN＝[0 0 1]^T时，||AN||能取到最小值，最小奇异值对应的奇异向量就是平面的法方向，也就求出了系数a,b，c：

最终可以得出

Step3：在每个样本集中计算数据点x的偏好函数，每个数据点x都用偏好函数φ来描述。偏好函数在取值区间为闭合区间[0,1]，最喜欢被表示成1，最不喜欢被表示成0。通过阈值的截断效果来计算数据x的偏好函数，如下：

τ为时间常数能够实现截断效果，d(x，h)表示数据点x与表面模型h之间的残差。

Step4：通过数据点x与表面模型h之间的残差公式d(x,h)计算点到平面之间的距离，又因为

因此距离计算如下：

d(x，h)＝|ax_i+by_i+cz_i+d| (11)

Step5：如果数据子集S＝{x}，则S的偏好函数就是x的偏好函数表示为q_x，否则S的偏好函数是向量p∈[0,1]^m，使其第i个分量为x的最小偏好函数

然后通过Tanimoto距离计算每个点偏好函数[0,1]^m之间的相似性：

其中，符号<·,·>表示

的标准内积。Tanimoto距离范围是[0,1]，对于

有如果p＝q，那么d_T(p，q)＝0；如果p⊥q，那么d_T(p,q)＝1。

Step6：如果任意p和q都正交，那么输出集群标签及其对应的表面模型，聚类结束；否则通过公式d_T(p,q)＝min_r,sd_T(r,s)找到p和q，并用二者的并集替换这两个集群，重新标记新集群，返回Step2。

步骤四、多视DSM融合

在获得区域分割结果及其对应的表面模型之后，将T-Linkage聚类获得的N个集群作为DSM高程值融合的基本单位，利用表面模型的约束作用在每个集群区域中执行DSM的融合以获得最终的DSM，融合过程如图4所示，融合结果如图5所示，详细步骤如下：

Step1：初始化集群标签i＝1，融合结果

Step2：对于第i个集群的3-D点集合

中的点p，若p到虚拟表面h_i的距离小于阈值τ_i，则将像素点p收入符合模型的点集合

计算如下：

Step3：通过中值滤波的方法计算最终的DSM，将点集合

作为中值滤波的作用范围，中值滤波计算如下：

其中，将

作为第i个集群的高程值，得到第i个集群的DSM：

Step4：合并所有集群中3-D点的中值滤波结果，获得最终融合的DSM：

i＝i+1 (16)

Step5：如果当前集群标签i＝N+1，那么输出融合结果

否则返回Step2。

以上所述仅为本发明的优选实施例而已，并不用于限制本发明，对于本领域的技术人员来说，本发明可以有各种更改和变化。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (3)

1.一种多视卫星影像数字表面模型融合方法，其特征在于，包括以下步骤：

步骤1：卫星影像与DSM的配准，通过卫星影像配置文件中的RPC系数构建RFM，建立图上坐标与大地坐标的直接映射关系，并根据大地坐标获得图上目标区域对应的DSM，在目标区域中实现卫星影像与DSM的配准；

步骤2：卫星影像及其对应的DSM联合数据超像素分割；将LSC超像素使用的十维映射函数扩展为十一维映射函数，然后对联合数据进行映射，并对十一维的像素进行利用谱聚类，获得卫星影像和DSM的联合数据的超像素分割结果，实现具有相同高程结构区域的聚集；

步骤3：在每个超像素区域内，利用DSM的3-D坐标进行T-Linkage表面模型聚类，在超像素区域内获得集群标签及其对应的表面模型；

步骤3.1：通过均匀采样与局部采样相结合的随机抽样获取估计表面模型所需的最小数据3-D点集，生成m个表面模型H＝{h_j}_j＝1...m；

步骤3.2：拟合表面模型h；

通过奇异值分解SVD的方法获得集群3-D空间中的表面模型ax+by+cz+d＝0；根据模型性质，系数a,b,c构成的向量N＝[a,b,c]^T为表面模型的法向量；首先求最小样本集中三维点{(x_i,y_i,z_i)}_i＝1...n的质心

n≥3，根据法向量的性质构建方程组：

令方程组左边矩阵为A，右边矩阵为法向量N，则拟合的目标函数为min||AN||，并且保证||N||＝1，则对A做奇异值分解A＝UDV^T；D为对角阵，U和V均为酉矩阵，根据酉矩阵的性质能得到：

||AN||＝||UDV^TN||＝||DV^TN||

当满足条件||V^TN||＝||N||＝1时，D的对角元素为奇异值，最后一个对角元素为最小奇异值；当且仅当V^TN＝[0 0 1]^T时，||AN||能取到最小值，最小奇异值对应的奇异向量就是平面的法方向，也就求出了系数a,b,c：最终可以得出

步骤3.3：在每个样本集中计算数据点x的偏好函数，每个数据点x都用偏好函数φ来描述；偏好函数在取值区间为闭合区间[0,1]，最喜欢被表示成1，最不喜欢被表示成0；通过阈值的截断效果来计算数据x的偏好函数；

τ为时间常数能够实现截断效果；d(x,h)表示数据点x与表面模型h之间的残差；

步骤3.4：通过数据点x与表面模型h之间的残差公式d(x,h)计算点到平面之间的距离；

d(x,h)＝|ax_i+by_i+cz_i+d|

步骤3.5：如果数据子集S＝{x}，则S的偏好函数就是x的偏好函数表示为q_x；否则S的偏好函数是向量p∈[0,1]^m，使其第i个分量为x的最小偏好函数

然后通过Tanimoto距离计算每个点偏好函数[0,1]^m之间的相似性：

其中，符号<·,·>表示

的标准内积；Tanimoto距离范围是[0,1]，对于

有如果p＝q，则d_T(p,q)＝0；如果p⊥q，则d_T(p,q)＝1；

步骤3.6：如果任意p和q都正交，那么输出集群标签及其对应的表面模型，聚类结束；否则通过公式d_T(p,q)＝min_r,sd_T(r,s)找到p和q，并用二者的并集替换这两个集群，重新标记新集群，返回步骤3.2；

步骤4：将T-Linkage聚类获得的N个集群作为DSM高程值融合的基本单位，利用表面模型的约束作用在每个集群区域中执行DSM的融合以获得最终的DSM；

步骤4.1：初始化集群标签i＝1，融合结果

步骤4.2：对于第i个集群的3-D点集合

中的点p，若p到虚拟表面h_i的距离小于阈值τ_i，则将像素点p收入符合模型的点集合

步骤4.3：通过中值滤波的方法计算最终的DSM，将点集合

作为中值滤波的作用范围，中值滤波计算如下：

其中，将

作为第i个集群的高程值，得到第i个集群的DSM：

步骤4.4：合并所有集群中3-D点的中值滤波结果，获得最终融合的DSM：

步骤4.5：令i＝i+1，如果当前集群标签i＝N+1，那么输出融合结果

否则返回步骤4.2。

2.根据权利要求1所述的一种多视卫星影像数字表面模型融合方法，其特征在于：所述的步骤1中卫星影像与DSM的配准的方法具体为：

步骤1.1：通过RFM建立卫星影像坐标与大地坐标的直接映射关系；

获取构建RFM的80个RPC系数，通过这80个RPC系数构建图上坐标(l,s)与大地坐标

之间的关系；令(B_n,L_n,H_n)为一个点的正则化的经度、纬度和高程坐标；

L_n＝(λ-L_OFF)/L_SCALE

H_n＝(h-H_OFF)/H_SCALE

其中，B_OFF、L_OFF和H_OFF为卫星配置文件中的正则化偏移参数；B_SCALE、L_SCALE和H_SCALE为卫星配置文件中的正则化比例参数；

有理函数的多项式方程为：

其中，L_num(·)、L_den(·)、S_num(·)和S_den(·)是三阶多项式函数并且拥有20个RPC系数；

步骤1.2：根据RFM将图像上坐标转换为经纬度坐标，再根据经纬度坐标在DSM上获取对应目标区域的DSM，实现卫星影像与多幅DSM之间的配准。

3.根据权利要求1或2所述的一种多视卫星影像数字表面模型融合方法，其特征在于：所述的步骤2中卫星影像及其对应的DSM联合数据超像素分割的方法具体为：

步骤2.1：将卫星影像的颜色转化为CIELAB颜色空间中的颜色分量，构建每个像素点p的六维坐标p＝(l_p,α_p,β_p,x_p,y_p,e_p)；其中l_p、α_p和β_p是像素点p在CIELAB颜色空间中的颜色分量值；x_p和y_p是图像平面中像素点p的水平和垂直坐标；e_p是像素点p的高程值；

步骤2.2：根据映射函数φ_e将每个六维像素点p映射到十一维；

其中，ω为权重；C_c和C_s用来控制颜色和空间信息的相对重要性，C_s/C_c的值越大空间上越接近的像素越倾向于聚成一类；C_e用来平衡感官信息与高程信息之间的重要性，C_e/C_c的值越大高程相近的像素越倾向于聚成一类；

步骤2.3：初始化加权K-means聚类；赋予每个像素点的权重ω(p)，令K表示聚类数，并在整个图像上分别以水平间隔v_x和垂直间隔v_y对K个种子像素进行均匀采样；π_k表示第k个聚类，k＝1,2,...,K，c_k是该类的聚类中心；初始化每个像素p的标签L(p)＝0，初始化每个像素p的距离d(p)＝∞；

步骤2.4：加权K-means聚类，通过欧几里得距离来计算每个像素点φ_e(p)与聚类中心c_k之间的距离d_e(p)；

其中，聚类中心c_k的计算如下式：

步骤2.5：如果d_e(p)＜d(p)，则更新距离d(p)＝d_e(p)和标签L(p)＝k；

步骤2.6：更新所有的像素的标签和聚类中心，如果聚类中心不再发生变化，出现收敛，则输出当前聚类标记结果，否则返回步骤2.4。

Images