技术领域
[0001] 本
发明涉及一种SFS三维重建加密稀疏DEM方法，尤其是涉及一种顾及地表光谱信息的SFS三维重建加密稀疏DEM方法。
背景技术
[0002] 数字高程模型(DEM)作为地形
基础数据，在测绘遥感、军事等方面具有广泛应用，并直接受到DEM
精度和
分辨率的影响。此外地质灾害应急响应救援也需要高分辨率的地表三维
可视化模型，能够辅助突发事故及灾害的救援工作准确进行，因此建立高分辨率的DEM具有极为重要的意义。
[0003] 传统的DEM建立是根据同一
传感器对同一地区获取
立体像对，基于摄影测量的原理产生的。然而当在某些地区还无法快速有效的获取高分辨率的立体像对时，高分辨率的DEM建立显得十分困难。另一方面，经典插值方法是获取高分辨率DEM的常用策略，通过保持内插面的连续性，在原有粗分辨率数据的基础之上得到更高分辨率的DEM格网。目前比较有效的DEM数据内插方法如克里金插值等方法都只基于已知的高程数据，无法增加新的细节信息，一旦地面控制点或
采样点丢失了某个重要特征高程，插值方法就无法恢复这个区域地形起伏变换，因此传统插值方法对加密DEM有很大局限性。
[0004] 明暗恢复形状(SFS)是三维信息提取的一种经典方法，也是
计算机视觉中三维形状恢复问题的关键技术，其目的是利用单幅图像中物体表面光
密度明暗变化来恢复其表面各点的相对高度或表面法矢量方向等参数值，为进一步对物体进行三维重构奠定基础。自70年代以来国外学者陆续展开SFS的理论研究工作，horn于1979年提出了SFS思想，后续其他学者不断提出SFS问题的各类解法，Zhang(1999年)对经典
算法进行总结，基于此Jean-Denis Durou(2008年)对SFS解算体系进行了进一步完善。
[0005] 根据SFS算法的思路，本
专利利用地形
地貌造成的遥感影像光谱变化，应用明暗恢复形状的原理，结合原有粗分辨率DEM数据和稀疏控制点高程数据，依据数学模型反演得到加密DEM三维高程，既能满足立体像对无法获取地区构建DEM的要求，又顾及了DEM格网内部地表真实变化趋势，对于快速、有效的建立高分辨率网格DEM具有重要意义。同时利用遥感影像地物的光谱信息，修复图像阴影区域，恢复图像上各点的真实
亮度值，对不同地物进行识别，精确得到
各向异性地表下不同地物类别的反射率值；基于各个地物内部进行同质区域内三维重建，提高单幅影像三维重建的精度，反演地表真实起伏，减少SFS引入的误差。
[0006] Rajabi(2002年)首先将SFS算法引入DEM格网内插，基于单幅图像较精确的恢复地表三维表面形状，配合稀疏地面控制点，能克服立体像对与传统插值方法的缺点，生成客观反映地形的高分辨率DEM。然而Rajabi(2003，2004年)均是基于灰度图像，在已知
光源条件下解算SFS反射图方程，并用双线性内插结果作为SFS算法初值，结果误差较大，反映地形的结构信息容易产生
混叠模糊。
发明内容
[0007] 本发明主要是解决已有SFS技术中单幅灰度图像光源估算方法精度低，也没有考虑图像阴影和植被
覆盖对于单幅遥感影像三维重建造成的误差等科学问题。本专利将光源估算方法扩展到多光谱图像领域，提出了一种利用多光谱影像代替全色图像，运用地物光谱信息辅助提高重建精度；同时通过抑制地表植被覆盖类型和地物阴影造成建立地表起伏引起与影像光谱亮度变换关系的影响；顾及地表光谱信息的SFS三维重建加密稀疏DEM方法。
[0008] 本发明另一目的是解决
现有技术所存在的常用DEM插值算法精度低，并不能完成基于稀疏控制点的大范围DEM快速建立技术问题；提供了一种基于稀疏控制点的大范围DEM快速建立的一种顾及地表光谱信息的SFS三维重建加密稀疏DEM方法。
[0009] 本发明的上述技术问题主要是通过下述技术方案得以解决的：
[0010] 一种顾及地表光谱信息的SFS三维重建加密稀疏DEM方法，其特征在于，包括以下步骤：
[0011] 步骤1，给出在复杂地表覆盖情况下三维表面上的点集高程值在x方向、y方向上的梯度与该坐标点(x，y)处光谱亮度值的关系，构建不同假设条件下的反射图模型，进而建立面向不同地物覆盖类型的反射图方程，所述数据包括：粗分辨率DEM格网，同一地区高分辨率遥感影像；
[0012] 步骤2，对步骤1所建立的反射图模型进行预处理；
[0013] 步骤3，对已完成步骤3的反射图模型反演地表法向量，针对不同的地物对象，利用基于对象的最小化方法反演地物内的地表法向量；
[0014] 步骤4，利用粗分辨率DEM格网点，恢复影像与真实地表的空间缩放尺度，运用滑动最小二乘方法来求解局部点；
[0015] 步骤5，结合步骤3获取的地表法向量以及步骤4恢复影像与真实地表的空间缩放尺度进行多约束条件下对DEM数据内插：在优化内插过程中，引入滑动最小二乘算法进行内插的
迭代计算，并用保持结构纹理的ICBI算法处理不同地物对象间的接边内插，产生高分辨率DEM格网。
[0016] 在上述的一种顾及地表光谱信息的SFS三维重建加密稀疏DEM方法，所述的反射图方程基于公式：
[0017] I(x，y)＝R(p，q)
[0018] 并假设：
[0020] 2)太阳光源可看作是无穷远处点光源平行投影；
[0022] 其中，原始光学影像数据记为I(x，y)，根据光源参数(Sx，Sy，Sz)和反射率albedo恢复得到的反射图记为R，p，q是表面方向，是z在x，y方向的导数：
[0023] 朗伯体反射模型R的具体表达式如下：
[0024] R(p，q)＝albedo*(Sz-p*Sx-q*Sy)/sqrt(p2+q2+1)
[0025] 设表面点(x，y，z(x，y))处的法线为N＝(cosαsinβ，sinαsinβ，cosβ)，光源方向为L＝(cosτsinγ，sinτsinγ，cosγ)，τ为光源偏
角，γ为光源倾角，α为法线偏角，β为法线倾角，上式即为：
[0026]
[0027]
[0028] 在上述的一种顾及地表光谱信息的SFS三维重建加密稀疏DEM方法，所述的步骤2中，对步骤1所建立的反射图模型进行预处理包括以下步骤：
[0029] 步骤3.1，对多光谱遥感影像与DEM数据进行裁剪与配准，使阴影信息对应的
位置匹配；
[0030] 步骤3.2，检查多光谱遥感影像上是否有存在光谱亮度，即DN值为零的区域，由于亮度为零的格网缺少二阶属性，被视为阴影，在步骤3.4进行该区域影像修复。
[0031] 步骤3.3，对完成步骤3.2的遥感影像进行的基于图像光谱与纹理的地物分割分类；对于分割后图像，估算光源方向以及各个对象的表面反射率；
[0032] 步骤3.4，对已完成步骤3.3的遥感影像利用
颜色恒常理论模型进行图像修复，滤除阴影区域；
[0033] 步骤3.5，对已完成步骤3.4的遥感影像进行各向异性条件下明暗恢复形状：对象内加入灰度梯度约束，光滑性约束和可积性约束，进行变分法和快速行进法解算，得到地表相对高程。
[0034] 本发明创造性的提出利用多光谱遥感影像结合SFS三维重建结果的粗分辨率DEM格网内插方法；提出利用地物光谱信息估算不同地物类型的反射率，去除地表植被覆盖类型对SFS重建精度的影响；提出基于阴影去除的单幅多光谱遥感影像三维重建方法，减少图像阴影对地物光谱亮度的影响；建立各向异性复杂地表的光源估算，引入多光谱信息，添加地面辅助控制点，提高光源方向估算精度；基于三维重建的结果，利用滑动最小二乘进行高程值内插，配合ICBI内插的地物接边方法，避免了DEM内插方法的过平滑现象，提高DEM加密精度。专利提出的方法主要应用于应急情况下三维数字地表模型(DSM)快速建立、复杂区域的高分辨率
三维建模、无人区或数据匮乏地区的DEM格网建立。
[0035] 在上述一种顾及地表光谱信息的SFS三维重建加密稀疏DEM方法，所述的步骤3.1中，对遥感影像与DEM数据进行裁剪与配准的具体步骤如下：
[0036] 1)读取DEM数据以及遥感影像数据文件头，得到其数据的覆盖范围；
[0037] 2)建立仿射变换模型
[0038]
[0039] 3)根据影像上下左右四个角点范围值，带入仿射变换模型，得到八个方程，求解6个变换模型参数；
[0040] 4)计算每一点的变换后的地理坐标。
[0041] 在上述的一种顾及地表光谱信息的SFS三维重建加密稀疏DEM方法，所述的步骤3.3中，地物分割分类的具体步骤如下：
[0042] 1)依照光谱信息结合地物纹理特性，定义类别个数c，价值函数
阈值ε，最大迭代次数t；在此步骤中本专利给出的实例(见图2)，定义类别个数为3，价值函数阈值为0.1，最大迭代次数为50；
[0043] 2)用值在0-1间的随机数初始化隶属矩阵U，使其满足下式中的约束条件：
[0044]
[0045] 3)用下式计算c个聚类中心ci，i＝1，...，c：
[0046]
[0047] 其中，xi(i＝1，2，...，n)为n个样本点光谱向量；
[0048] 4)根据下式计算价值函数：如果价值函数小于确定的阈值ε，或价值函数相对上次价值函数值的改变量小于阈值，则算法停止：
[0049]
[0050] 这里uij介于0到1间；ci为模糊组I的聚类中心，dij＝||ci-xj||为第i个聚类中心与第j个数据点间的欧几里德距离；且m∈[1，∞)是一个加权指数，这里m取值为2；
[0051] 若价值函数大于阈值，则用下式计算新的U矩阵；返回步骤2)；
[0052]
[0053] 5)若迭代次数超过阈值t，则计算停止，最终得到多光谱遥感影像的地物分类图。
[0054] 在上述的一种顾及地表光谱信息的SFS三维重建加密稀疏DEM方法，所述的步骤3.3中，对分割后图像，估算光源方向的具体步骤如下：
[0055] 1)利用线性函数转换对图像亮度值进行归一化：
[0056] Ii′＝Ii/Imax
[0057] 其中Ii为i点处的某一波段原始亮度值，Imax为该波段全局影像亮度最大值，Ii′为归一化后的亮度值；
[0058] 2)取同一类别地物的原始DEM格网高程值及其所在位置对应遥感影像上的亮度值，作为控制点，建立多个反射图方程：
[0059] Ii＝albedo*(Sz-pi*Sx-qi*Sy)/sqrt(pi2+qi2+1)；
[0060] 其中i∈N为DEM格网点个数，Ii为该点的亮度值，(p，q)为该点的x，y方向的梯度，对于同一类地物，反射率albedo为常量；即利用N个方程(N＞3)解算3个光源参数(Sx，Sy，Sz)；利用最小二乘法得到参数估计值。
[0061] 在上述的一种顾及地表光谱信息的SFS三维重建加密稀疏DEM方法，所述的步骤3.3中，获取各个对象的表面反射率的具体步骤如下：其中，各个对象是指不同的地物类别：
[0062]
[0063] 其中E为总体均值，
[0064]
[0065] 令：
[0066]
[0067]
[0068]
[0069] 若 则γ＝0；否则，γ为 的解；
[0070] 对于反射率： 其中，I为图像上同一类型地物的归一化亮度值。
[0071] 在上述的一种顾及地表光谱信息的SFS三维重建加密稀疏DEM方法，所述的步骤3.4中，进行图像修复，滤除阴影区域的具体步骤如下：
[0072] 1)根据多光谱遥感影像各波段DN值及结构信息进行阴影区域检测：认为在所有波段上的反射率均接近于0的区域为阴影区域；
[0073] 2)读入原图像I和检测出的阴影区域模板，根据模板将原图像分成两部分：阴影区域图像Is和非阴影区域图像Ins；
[0074] 3)分别对阴影区域图像Is和非阴影区域图像Ins用Shade of Gray算法进行颜色恒常性计算，得到阴影区域光照颜色估计[s_r，s_g，s_b]和非阴影区光照颜色估计[n_r，n_g，n_b]，计算公式为：
[0075]
[0076] 式中， 表示当前图像中的光源颜色，f＝{R，G，B}表示图像各波段的DN值，这里假设为红，绿，蓝三个波段，p＝6；
[0077] 各个波段的校正增益为：
[0078]
[0079]
[0080]
[0081] 4)阴影区域变换到非阴影区域的比率常向量为
[0082] 5)将比率常向量k转换到对数空间
[0083] 6)把阴影区域图像Is转换到对数空间Ls＝log(Is)；
[0084] 7)得到去除阴影后的图像，I＝exp(Ls+c)+Ins。
[0085] 在上述的一种顾及地表光谱信息的SFS三维重建加密稀疏DEM方法，进行明暗恢复形状的具体步骤如下：
[0086] 1)将SFS反射图方程改写为地表梯度的函数，并加入亮度约束，光滑性约束等条件来求解方程；用图像和反射图的整体误差构造目标泛函，并使用变分法求解三维表面函数；
[0087] 已知图像灰度与反射模型所确定的物体表面图像亮度之间可能存在误差，将亮度方程变为误差函数的形式，即亮度约束：
[0088]
[0089] 2)引入灰度梯度约束和可积性约束，建立误差方程，得：2 2 2 2 2
[0090] ∫ ∫ F(p，q，z)dxdy＝ ∫ ∫((I-R)+(Rx-Ix)+(Ry-Iy)+μ((zx-p)+(zy-q)))dxdy；
[0091] 3)利用欧拉方程求泛函数极小值，通过逐次迭代，求出地表真实法向量(p，q)和高程z。
[0092] 在上述的一种顾及地表光谱信息的SFS三维重建加密稀疏DEM方法，所述的步骤5中，执行下列两个步骤：
[0093] 步骤5.1，对于同一个地物对象，即保持了对象内地表变化的连续性和光滑性，利用遥感影像三维重建的结果，根据下列公式得到地表坡度坡向图：
[0094]
[0095]
[0096] 其中，(fx，fy)为三维重建的相对高度在x，y方向上的梯度值；
[0097] 建立代价函数，使内插后高程格网地表坡度坡向与三维重建的地表坡度坡向图的差最小，用滑动最小二乘法迭代解算DEM格网；
[0098] 步骤5.2，若相邻两次迭代结果之差小于阈值0.1，则用ICBI算法处理不同地物对象间的接边内插，产生高分辨率DEM格网，即：
[0099] 设格网内的偶数位点为高程已知点，对于位置为(2i+1，2j+1)点，首先计算两个方向上的值：
[0100] I11(2i+1，2j+1)＝I(2i-2，2j+2)+I(2i，2j)+I(2i+2，2j-2)
[0101] -3I(2i，2j+2)-3I(2i+2，2j)+I(2i，2j+4)+I(2i+2，2j+2)+I(2i+4，2j)；
[0102] I22(2i+1，2j+1)＝I(2i，2j-2)+I(2i+2，2j)+I(2i+4，2j+2)
[0103] -3I(2i，2j)-3I(2i+2，2j+2)+I(2i-2，2j)+I(2i，2j+2)+I(2i+2，2j+4)；
[0104] 最终其内插值按下式确定：
[0105]
[0106] 因此，本发明具有如下优点：
[0107] 1.利用多光谱影像代替
全色影像，运用地物光谱信息辅助提高三维重建精度；能有效抑制地表植被覆盖类型、地物阴影对于建立地表三维模型与影像光谱亮度变换关系的影响；
[0108] 2.相对于常用DEM插值算法内插精度有显著提高，实例验证能将粗网格DEM空间分辨率提高2倍，并能完成基于稀疏控制点的大范围DEM快速建立。
附图说明
[0110] 图2a1为试验区域多光谱遥感图像；
[0111] 图2a2为试验区域原始DEM影像(2.5米)；
[0112] 图2a3为分类结果，分别为
岩石、植被、
水体；
[0113] 图2b1为本专利SFS算法结果；
[0114] 图2b2为传统SFS算法结果；
[0115] 图2c1为本发明DEM加密结果；
[0116] 图2c2为图图2c1的三维显示；
[0117] 图2d1为本发明加权最小二乘曲面拟合内插DEM结果的DEM加密格网；
[0118] 图2d2为图2d1的三维显示；
[0119] 图2e1为倒距离加权内插DEM结果的DEM加密格网；
[0120] 图2e2为图图2e1的三维显示；
[0121] 图2f1为克里金插值DEM加密；
[0122] 图2f2为图图2f1的三维显示DEM格网；
[0123] 图3a为各插值算法生成DEM高程比较；
[0124] 图3b1为本发明算法的等高线；
[0125] 图3b2为克里金算法等高线；
具体实施方式
[0126] 下面通过
实施例，并结合附图，对本发明的技术方案作进一步具体的说明。
[0127] 实施例：
[0128] 如图1所示，一种顾及地表光谱信息的SFS三维重建加密稀疏DEM方法，包括以下步骤：
[0129] 步骤1，给出在复杂地表覆盖情况下三维表面上的点集高程值在x方向、y方向上的梯度与坐标点(x，y)处光谱亮度值的关系，构建不同假设条件下的反射图模型，进而建立面向不同地物覆盖类型的反射图方程，所述数据包括：粗分辨率DEM格网，同一地区高分辨率遥感影像；反射图方程基于公式：
[0130] I(x，y)＝R(p，q)
[0131] 并假设：
[0132] 1)地形表面为朗伯体模型；
[0133] 2)太阳光源可看作是无穷远处点光源平行投影；
[0134] 3)成像几何关系为正交投影；
[0135] 其中，原始光学影像数据记为I(x，y)，根据光源参数(Sx，Sy，Sz)和反射率albedo恢复得到的反射图记为R，p，q是表面方向，是z在x，y方向的导数：
[0136] 朗伯体反射模型R的具体表达式如下：
[0137] R(p，q)＝albedo*(Sz-p*Sx-q*Sy)/sqrt(p2+q2+1)
[0138] 设表面点(x，y，z(x，y))处的法线为N＝(cosαsinβ，sinαsinβ，cosβ)，光源方向为L＝(cosτsinγ，sinτsinγ，cosγ)，τ为光源偏角，γ为光源倾角，α为法线偏角，β为法线倾角，上式即为：
[0139]
[0140]
[0141] 步骤2，对步骤1所建立的反射图模型进行预处理；预处理包括以下步骤：
[0142] 步骤3.1，对多光谱遥感影像与DEM数据进行裁剪与配准，使阴影信息对应的位置匹配；具体步骤如下：
[0143] 1)读取DEM数据以及遥感影像数据文件头，得到其数据的覆盖范围；
[0144] 2)建立仿射变换模型
[0145]
[0146] 3)根据影像上下左右四个角点范围值，带入仿射变换模型，得到八个方程，求解6个变换模型参数；
[0147] 4)计算每一点的变换后的地理坐标。
[0148] 步骤3.2，检查多光谱遥感影像上是否存在光谱亮度(即DN值)为零的区域，由于亮度为零的格网缺少二阶属性，被视为阴影，在步骤3.4进行该区域影像修复；
[0149] 步骤3.3，对完成步骤3.2的遥感影像进行的基于图像光谱与纹理的地物分割分类；对于分割后图像，估算光源方向以及各个对象的表面反射率；地物分割分类的具体步骤如下：
[0150] 1)依照光谱信息和地物纹理特性定义类别个数c，价值函数阈值ε，最大迭代次数t；在此步骤中本专利给出的实例定义类别个数为3，价值函数阈值为0.1，最大迭代次数为50；
[0151] 2)用值在0，1间的随机数初始化隶属矩阵U，使其满足下式中的约束条件：
[0152]
[0153] 3)用下式计算c个聚类中心ci，i＝1，...，c：
[0154]
[0155] 其中，xi(i＝1，2，...，n)为n个样本点光谱向量；
[0156] 4)根据下式计算价值函数：如果价值函数小于确定的阈值ε，或价值函数相对上次价值函数值的改变量小于阈值，则算法停止：
[0157]
[0158] 这里uij介于0到1间；ci为模糊组I的聚类中心，dij＝||ci-xj||为第i个聚类中心与第j个数据点间的欧几里德距离；且m∈[1，∞)是一个加权指数，这里m取值为2；
[0159] 若价值函数大于阈值，则用下式计算新的U矩阵；返回步骤2)；
[0160]
[0161] 5)若迭代次数超过阈值t，则计算停止，最终得到多光谱遥感影像的地物分类图。
[0162] 然后，对分割后图像，估算光源方向的具体步骤如下：
[0163] 1)利用线性函数转换对图像亮度值进行归一化：
[0164] Ii′＝Ii/Imax
[0165] 其中Ii为i点处的某一波段原始亮度值，Imax为该波段全局影像亮度最大值，Ii′为归一化后的亮度值；
[0166] 2)取同一类别地物的原始DEM格网高程值及其所在位置对应遥感影像上的亮度值，作为控制点，建立多个反射图方程：
[0167] Ii＝albedo*(Sz-pi*Sx-qi*Sy)/sqrt(pi2+qi2+1)；
[0168] 其中i∈N为DEM格网点个数，Ii为该点的亮度值，(p，q)为该点的x，y方向的梯度，对于同一类地物，反射率albedo为常量；即利用N个方程(N＞3)解算3个光源参数(Sx，Sy，Sz)；利用最小二乘法得到参数估计值。
[0169] 最后，获取各个对象的表面反射率的具体步骤如下：其中，各个对象是指不同的地物类别：
[0170]
[0171] 其中E为总体均值，
[0172]
[0173] 令：
[0174]
[0175]
[0176]
[0177] 若 则γ＝0；否则，γ为 的解；
[0178] 对于反射率： 其中，I为图像上同一类型地物的归一化亮度值。
[0179] 步骤3.4，对已完成步骤3.3的遥感影像利用颜色恒常理论模型进行图像修复，滤除阴影区域；具体步骤如下：
[0180] 1)根据多光谱遥感影像各波段DN值及结构信息进行阴影区域检测：认为在所有波段上的反射率均接近于0的区域为阴影区域；
[0181] 2)读入原图像I和检测出的阴影区域模板，根据模板将原图像分成两部分：阴影区域图像Is和非阴影区域图像Ins；
[0182] 3)分别对阴影区域图像Is和非阴影区域图像Ins用Shade of Gray算法进行颜色恒常性计算，得到阴影区域光照颜色估计[s_r，s_g，s_b]和非阴影区光照颜色估计[n_r，n_g，n_b]，计算公式为：
[0183]
[0184] 式中， 表示当前图像中的光源颜色，f＝{R，G，B}表示图像各波段的光谱亮度DN值，这里假设为红，绿，蓝三个波段，p＝6；
[0185] 各个波段的校正增益为：
[0186]
[0187]
[0188]
[0189] 4)阴影区域变换到非阴影区域的比率常向量为
[0190] 5)将比率常向量k转换到对数空间
[0191] 6)把阴影区域图像Is转换到对数空间Ls＝log(Is)；
[0192] 7)得到去除阴影后的图像，I＝exp(Ls+c)+Ins。
[0193] 步骤3.5，对已完成步骤3.4的遥感影像进行各向异性条件下明暗恢复形状：对象内加入灰度梯度约束，光滑性约束和可积性约束，进行变分法和快速行进法解算，得到地表相对高程，明暗恢复形状具体步骤如下：
[0194] 1)将SFS反射图方程改写为地表梯度的函数，并加入亮度约束，光滑性约束等条件来求解方程；用图像和反射图的整体误差构造目标泛函，并使用变分法求解三维表面函数；
[0195] 已知图像灰度与反射模型所确定的物体表面图像亮度之间可能存在误差，将亮度方程变为误差函数的形式，即亮度约束：
[0196]
[0197] 2)引入灰度梯度约束和可积性约束，建立误差方程，得：
[0198] ∫ ∫ F(p，q，z)dxdy＝ ∫ ∫((I-R)2+(Rx-Ix)2+(Ry-Iy)2+μ((zx-p)2+(zy-q)2))dxdy；
[0199] 3)利用欧拉方程求泛函数极小值，通过逐次迭代，求出地表真实法向量(p，q)和高程z。
[0200] 步骤3，对已完成步骤3的反射图模型反演地表法向量，针对不同的地物对象，利用基于对象的最小化方法反演地物内的地表法向量；
[0201] 步骤4，利用粗分辨率DEM格网点，恢复影像与真实地表的空间缩放尺度，运用滑动最小二乘方法来求解局部点；
[0202] 步骤5，结合步骤3获取的地表法向量以及步骤4恢复影像与真实地表的空间缩放尺度进行多约束条件下对DEM数据内插：在优化内插过程中，引入滑动最小二乘算法进行内插的迭代计算，并用保持结构纹理的ICBI算法处理不同地物对象间的接边内插，产生高分辨率DEM格网，本步骤具体执行下列两个步骤：
[0203] 步骤5.1，对于同一个地物对象，即保持了对象内地表变化的连续性和光滑性，利用遥感影像三维重建的结果，根据下列公式得到地表坡度坡向图：
[0204]
[0205]
[0206] 其中，(fx，fy)为三维重建的相对高度在x，y方向上的梯度值；
[0207] 建立代价函数，使内插后高程格网地表坡度坡向与三维重建的地表坡度坡向图的差最小，用滑动最小二乘法迭代解算DEM格网；
[0208] 步骤5.2，若相邻两次迭代结果之差小于阈值0.1，则用ICBI算法处理不同地物对象间的接边内插，产生高分辨率DEM格网，即：
[0209] 设格网内的偶数位点为高程已知点，对于位置为(2i+1，2j+1)点，首先计算两个方向上的值：
[0210] I11(2i+1，2j+1)＝I(2i-2，2j+2)+I(2i，2j)+I(2i+2，2j-2)
[0211] -3I(2i，2j+2)-3I(2i+2，2j)+I(2i，2j+4)+I(2i+2，2j+2)+I(2i+4，2j)；
[0212] I22(2i+1，2j+1)＝I(2i，2j-2)+I(2i+2，2j)+I(2i+4，2j+2)
[0213] -3I(2i，2j)-3I(2i+2，2j+2)+I(2i-2，2j)+I(2i，2j+2)+I(2i+2，2j+4)；
[0214] 最终其内插值按下式确定：
[0215]
[0216] 如图2所示，采用已配准的SPOT遥感图像与DEM格网影像。遥感
图像空间分辨率为5m，原始DEM格网分辨率为2.5m。为验证算法精度，将DEM格网
下采样至空间分辨率为10m，并用遥感影像三维重建的方法以及DEM常用插值的方法恢复到5m分辨率，比较分析重建结果，验证方法的有效性。首先对遥感影像进行监督分类，将地表覆盖分为植被、水体、岩石和农田区域。对于不同的地表覆盖类型，分别用类间统计信息估算地物反射率和光源位置估算。对于不同覆盖类型区域，进行面向对象的SFS恢复，反演地表相对起伏状态。再根据地表坡度坡向信息，利用最小二乘结合ICBI方法内插高分辨率DEM格网。并用最小二乘拟合、倒距离加权插值，样条插值的结果进行比较分析，下表为本发明的单幅影像三维重建方法与传统内插算法精度对比：
[0217]均误差 最大误差 最小误差 绝对均误差
SFS重建 -0.1399 13 -38 0.4673
最小二乘拟合 0.2574 66 -69 1.7914
距离加权反比插值 -0.2638 12 -68 0.7487
克里金插值 0.6373 22 -21 1.3490
[0218] 图2b中左图为本发明SFS算法恢复的图像相对高度，右图为传统SFS算法恢复的相对高度。首先对于图(a)中遥感影像上存在的三类地物：岩石，植被，水体，本文分别估算三者的地表反射率值，在SFS解算中用地表反射参数滤除地物覆盖类型不同带来的图像灰度差异。因此本发明的SFS算法可以将植被和岩石处明显分开，计算出植被区域对于图像上部的岩石区域具有相对较低的高度值。而传统的SFS算法由于没有考虑地物覆盖类型的影响，导致结果中植被区域和岩石区域高度平滑相连，捕捉不到类别间变化的高程。
[0219] 在图2c为基于本发明算法重建的DEM格网，从三维显示中可以看到重建后的DEM细节突出，原始DEM数据中的地形结构，包括山脊，山谷，水体池塘，都能清晰看见，结构特征保持较好。
[0220] 图2d和e的结果显示，距离加权反比插值存在过平滑的问题，丢失很多地表起伏细节。最小二乘拟合则存在振荡，造成伪信息。同时由误差表-1看出，两者与本发明算法相比，具有较大的最大最小误差，与真实高分辨率DEM格网差异较大。
[0221] 克里金插值结果如图2f所示，可以看到DEM形态结构较清晰，比最小二乘拟合与距离加权反比插值过平滑现象较轻，但结合图3a和3b中可以看出克里金DEM加密的结果的细节信息仍不如本发明方法的结果丰富，这是由于本专利算法增加了高分辨率遥感影像光谱信息。
[0222] 从图3a中可以看出，最小二乘拟合存在振荡和伪信息，尤其当有噪声存在时干扰非常严重。距离加权反比插值存在过平滑的问题，丢失很多地表起伏细节，同时表-1看出，两者与本发明算法相比，具有较大的最大和最小误差，与真实高分辨率DEM格网差异较大。结合克里金插值结果，如图(f所示，可以看到DEM形态结构较清晰，过平滑现象较最小二乘拟合与距离加权反比插值较轻，但也存在一定程度的平滑现象。从图3、图2以及表-1可以看出，本发明增加了单幅遥感影像和地表光谱信息约束，利用SFS三维重建加密稀疏DEM方法，能够恢复更丰富的细节信息，与克里金插值算法等方法比较有显著优势。
[0223] 本发明能够快速建立复杂地表三维模型，主要用于高分辨率DEM格网快速获取，在应急情况下三维地表模型快速建立，卫星遥感影像目标识别和遥感影像测量方面的应用，复杂区域的高分辨率三维建模，无人区的DEM格网建立，
飞行器在复杂地形条件下的应急着陆，地表三维模型
数据压缩方面具有重要意义。
[0224] 本发明是一种顾及地表光谱信息的SFS三维重建加密稀疏DEM方法，具体包括：基于各向异性复杂地表的单幅遥感影像三维重建，利用异质地表
光谱特性去除地表植被覆盖类型对地表起伏引起的灰度变化的影响。对于不同地物类型局部估算地表反射率值，实现了多光谱影像的光源位置估算方法。去除地物遮挡阴影，避免图像光谱亮度值不能真实反映地表起伏的情况。明确建立复杂地表图像灰度与光源方向，地表法向量以及地表反射率之间的关系，利用滑动最小二乘方法来求解局部点和ICBI插值方法，利用粗分辨率DEM格网点，恢复影像与真实地表的缩放尺度，解算高密度网格DEM高程。
[0225] 本发明用于应急情况下三维地表模型快速建立，卫星遥感影像军事目标识别和遥感影像测量等方面，在复杂地形区域的高分辨率三维建模，数据匮乏区域的大范围DEM格网建立、DEM格网数据压缩、DEM与影像的三维一体化显示，飞行器在复杂地形下的应急着陆等方面也有广泛的应用前景。
[0226] 本专利所描述的具体实施例仅仅是对本发明思想作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的
修改或补充，或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附
权利要求书所定义的范围。
