CN101271596A - 基于曲面论基本定律对地形曲面建立数字模型的方法 - Google Patents

Abstract

本发明属于地球观测与导航技术，涉及地理信息系统处理，公开一种基于曲面论基本定律对地形曲面建立数字模型的方法，将曲面进行空间离散化，然后根据曲面论中的高斯方程建立曲面离散化格网方程，求解超定线性系统，即可获得曲面的数字模型。本发明提高了数字地面模型的精度，可应用于卫星影像校正、双星定位系统以及国防信息化建设。

Description

基于曲面论基本定律对地形曲面建立数字模型的方法

技术领域

本发明属于地球观测与导航技术，涉及地理信息系统处理，尤其涉及对曲面三维建模的方法。

背景技术

数字地面模型(DTM)可定义为诸如数字高程模型(DEM)、坡度模型和坡向模型等有关地形变量的数字表达。数字地面模型的精度问题是地理信息系统(Geographical Information System，GIS)领域近20年来的研究热点之一。数字高程模型的方法是影响数字地面模型精度的主要因素。

数字高程模型作为地理位置函数的地面高程表达，是所有其它数字地面模型的基础。一个数字高程模型可推演出许多数字地面模型。随着高空间分辨率卫星影像的发展，数字高程模型在遥感数据的获取中发挥着越来越重要的作用。数字高程模型通常被用于校正卫星影像，卫星影像的校正精度取决于数字高程模型的精度。因此，高精度数字高程模型研究对有高精度要求的数字地面模型应用问题显得尤其重要。

数字地面模型是现实世界中连续地表的逼近。数字地面模型的误差包括采样误差、采样仪器产生的误差、变换地面控制点产生的误差、数学模型带来的误差、数据源的传播误差、表达误差和栅格分辨率转换引起的误差等。数字地面模型的误差会通过建模过程被传播，并体现在最终成果中。地理信息系统中各种误差的累积使最终成果的可用性令人怀疑。尤其是地理信息系统最常用的叠合(overlay)操作对固有误差和运算误差都非常敏感，通过多次叠合获得的结果对有较高精度要求的应用毫无使用价值。地理信息系统的误差问题不是单纯的技术问题，而是地理信息系统的理论基础问题。

许多学者已对数字地面模型的误差问题进行了长期不懈的研究。例如，Goodchild(1982)将布朗分形过程引入地面模拟模型以提高数字地面模型的精度。Walsh等(1987)发现，通过识别输入数据的固有误差和运算误差，可以使总体误差达到最小。Hutchinson和Dowling(1991)为了构建反映流域自然结构的数字地面模型，引入了试图消除假深洼信息的流域强迫规则。Unwin(1995)在回顾了有关研究成果(Heuvelink et al.1989，Heuvelink and Burrough 1993)之后提出，检验地理信息系统在运算过程中误差传播的通用工具有助于提高数字地面模型的精度。Wise(2000)认为，为了提高数字地面模型的精度，当使用地理信息系统的时候，必须区分栅格模型和像元模型，存储在栅格中的信息只与网格的中心点有关，而存储在像元的值代表整个网格。美国地质调查局(1997)数字地面模型质量控制系统的主要内容包括精度统计检验、数据文件物理与逻辑格式检验和视觉检验。Shi等(2005)提出了减小数字地面模型误差的高次插值方法。Podobnikar(2005)认为，通过使用一切可用的数据源(甚至没有高度属性的低质量数据集)，可以提高数字地面模型的精度。

发明内容

现有技术这些方法都没能从根本上解决数字地面模型的误差问题，为了解决现有技术的问题，本发明的目的是建立高精度的数字地面模型(DTM)，为此，本发明提供一种基于曲面论基本定律的曲面建立数字模型的方法。

为了实现所述的目的，本发明基于曲面论基本定律对地形曲面建立数字模型的方法(A method of Surface Modeling based on Theory ofSurface---SMTS，以后简称SMTS)，步骤如下：

步骤S1：采集地形表面曲面数据；

步骤S2：对地形曲面的采样数据进行空间离散化处理，得到离散化的曲面数据Q(0)；

步骤S3：根据曲面论，计算地形曲面数据Q(K)，Q(K)基本要素为：地形曲面上每个点的第一类基本量(First Fundamental Form)的系数E、F、G和第二类基本量(Second Fundamental Form)的系数L、M、N；

步骤S4：判断迭代次数K，第一次迭代时，K等于0；当K等于0时，执行步骤S6；当K不等于0时，执行步骤S5；

步骤S5：判断迭代过程中相临两次曲面数据的差异Q(K)-Q(K-1)，当差异大于某一个阈值时，执行步骤S6；当Q(K)-Q(K-1)等于小于或等于某一个阈值时，执行步骤S7；

步骤S6：根据曲面论中的高斯方程(Gauss Equations)，结合地形曲面上点的数据参数，建立地形曲面的差分微分方程为：S×f＝T，求解差分微分方程，得到地形曲面的模拟值为：f＝S^-1×T，完成一次地形曲面的迭代过程，并且K＝K+1；其中S是曲面方程的系数矩阵，f是待求的曲面未知数，T为计算出的曲面常数项；

步骤S7：输出地形曲面数据Q(K)，得到高精度的数字地形模型。

采用离散化曲面的方法建立高斯方程(Gauss Equations)，并通过迭代计算的方法进行求解。

所述曲面由第一类基本量(First Fundamental Form)和第二类基本量决定(Second Fundamental Form)。

所述地形曲面为蒙格面片(Monge Patch)形式：

z＝f(x，y)，其中：x，y分别表示x坐标和y坐标。

所述地形曲面的第一类基本量(First Fundamental Form)的系数为：

$E=1+f_x^2$ F＝f_xf_y $G=1+f_y^2$ ，第一类基本量的系数E表示当变化x时，变化将离曲面的多远，G表示当变化y时，变化将离曲面多远，F表示当变化x，y同时变化时，他们之间的夹角是多大；

所述地形曲面的第二类基本量(Second Fundamental Form)的系数为：

$L=\frac{f_{\mathrm{xx}}}{\sqrt{1+f_x^2+f_y^2}},$

$M=\frac{f_{\mathrm{xy}}}{\sqrt{1+f_x^2+f_y^2}},$

$N=\frac{f_{\mathrm{yy}}}{\sqrt{1+f_x^2+f_y^2}},$

第二类基本量的系数L、M、N都与曲面的法向量曲率紧密相关的。

选择地形数学曲面f(x，y)＝2sin(πx)sin(πy)+1作为标准曲面，对标准曲面进行数值模拟分析，当标准曲面顶部没有采样点时，经典模型的模拟结果有峰值削平和边界振荡，曲面整体扭曲；将曲面整体扭曲经过HASM迭代后，获得被削平的峰值，消除边界振荡，使模拟曲面与标准曲面相同。

对地形数学曲面f进行数值模拟时，先根据采样点的值求解出基本量E、F、G、L、M、N，再利用基本量E、F、G、L、M、N求出曲面的第二类克里斯托弗尔变量Γ_i ^j，i＝11，12，22；j＝1，2，在给定的边界条件下，模拟得到曲面的数值解。

所述的方法用于空间图像数据处理。

本发明产生的技术效果：

本发明的SMTS目前主要用于建立高精度的数字地面模型(DTM)。是基于微分几何的曲面论的曲面建模方法，本发明这种空间数据的处理方式，采用经典的曲面论是本发明HASM-DTM的出发点，曲面论也是使SMTS与以往生成DTM的系统的不同之处。实验结果表明：本发明基于曲面论的曲面建模方法(SMTS)的精度远远高于GIS中常用的空间插值方法，如三次样条(Cubic Spline)插值方法，反距离权重(Inverse DistanceWeight，IDW)方法，克里金法(Kriging)，TIN(Triangular Irregular Network)基础上的线性插值(TLI)方法。当采样间距为2h、迭代次数为64时，高精度曲面建模迭代模拟方法HASM的平均绝对误差较插值模拟方法TIN减小了47469倍、较Cubic减小了4940倍、较Spline减小了2745倍、较IDW减小449519倍、较Kriging减小450529倍；HASM的平均相对误差较插值模拟方法TIN减小了52169倍、较Cubic减小了6205倍、较Spline减小了3706倍、较IDW减小573251倍、较Kriging插值模拟方法减小575048倍。

本发明应用于国防信息化建设：数字高程模型是国防信息化建设的重要基础数据平台之一，在战场准备、战场仿真、作战指挥、后勤保障等方面处于举足轻重的地位。在数字高程模型建设中，精度是其核心。如果精度得不到保障，数字高程模型将毫无意义。如果用来指挥决策，将会造成难以弥补的损失。高精度数字地面模型快速生成系统的研制完成，将为我国国防信息化建设做出重要贡献。

本发明的方法可用于遥感图象处理的空间数据的处理、医学图象处理等各个领域。

附图说明

图1是本发明HASM的基本原理图

图2标准曲面f(x，y)＝2sin(πx)sin(πy)+1及其等高线图

图3Kriging插值模拟结果与HASM迭代模拟结果比较

图4IDW插值模拟结果与HASM迭代模拟结果比较

图5TIN插值模拟结果与HASM迭代模拟结果比较

图6空间点数据离散化示意图

具体实施方式

下面将结合附图对本发明加以详细说明，应指出的是，所描述的实施例仅旨在便于对本发明的理解，而对其不起任何限定作用。

本发明应用于卫星影像校正：数字高程模型作为地理位置函数的地面高程表达，是所有其它数字地面模型的基础。一个数字高程模型可推演出许多数字地面模型。随着高空间分辨率卫星影像的发展，数字高程模型在遥感数据的获取中发挥着越来越重要的作用。数字高程模型通常被用于校正卫星影像，卫星影像的校正精度取决于数字高程模型的精度。因此，高精度数字高程模型研究对有高精度要求的数字地面模型应用问题显得尤其重要。

本发明应用于双星定位系统：由于双星定位导航系统定位是要依靠地球椭球面，所以对于实际应用来说，存在高度误差，必须依靠数字高程模型进行矫正。现用的数字高程模型是由经典模型构造的，因此，如果改用较经典模型精度高多个数量级的高精度数字高程模型，双星定位系统的精度将会大幅度提高。

本发明基于曲面论基本定律的曲面建模方法(SMTS)目前主要用于建立高精度的数字地面模型(DTM)。是基于微分几何的曲面论的曲面建模方法。经典的曲面论是本发明系统HASM-DTM的出发点，也是使SMTS与以往生成DTM的系统的不同之处。基于曲面论的曲面建模方法使得SMTS的精度远远高于其他方法。

本发明的方法，将曲面进行空间离散化，然后根据曲面论中的高斯方程(Gauss Equations)建立曲面离散化格网方程，求解超定线性系统，即可获得曲面的数字模型。具体如下：

采集地形表面曲面数据后，对地形曲面的采样数据进行空间离散化处理，得到离散化的曲面数据Q(0)；空间数据的离散化是HASM的重要第一步，也是极易产生误差的一步。离散化的示意图如上图6所示，其中三角形的点是随机的采样数据点，圆形的点是每个格子的中心点，理论上，采样点应该都落在每个格子的中心点也即圆形点上是最佳的。但是实际应用中，常常不能恰好落在中心点上，而且会偏离中心点很远。于是需要进行近似的离散化，将三角形的点当作是其最临近的圆形点的值，然后再进行模拟计算。

根据微分几何学中的曲面理论，计算地形曲面数据Q(K)(Q(K)表示曲面的离散化结果)，Q(K)基本要素为：地形曲面上每个点的第一类基本量的系数E、F、G和第二类基本量的系数L、M、N；空间曲线由曲率和挠率决定；判断迭代次数K，第一次迭代时，K等于0；当K等于0时，根据曲面论中的高斯方程(Gauss Equations)，结合地形曲面上点的数据参数，建立地形曲面的差分微分方程为：S×f＝T，求解差分微分方程，得到地形曲面的模拟值为：f＝S^-1×T，完成一次地形曲面的迭代过程，并且K＝K+1；其中S是曲面方程的系数矩阵，f是待求的曲面未知数，T为计算出的曲面常数项；当K不等于0时，判断迭代过程中相临两次曲面数据的差异Q(K)-Q(K-1)，当差异大于某一个阈值时，执行K等于0的步骤；当Q(K)-Q(K-1)等于小于或等于某一个阈值时，输出地形曲面数据Q(K)，得到高精度的数字地形模型。

曲面由第一类基本量和第二类基本量决定，如果曲面表达为z＝f(x，y)，其中：x，y分别表示x坐标和y坐标，则此曲面的第一类基本量可表达为：

$E=1+f_x^2$ F＝f_xf_y $G=1+f_y^2$ ，第一类基本量的系数E表示当变化x时，变化将离曲面的多远，G表示当变化y时，变化将离曲面多远，F表示当变化x，y同时变化时，他们之间的夹角是多大；

第二类基本量可表达为，

$L=\frac{f_{\mathrm{xx}}}{\sqrt{1+f_x^2+f_y^2}},$ $M=\frac{f_{\mathrm{xy}}}{\sqrt{1+f_x^2+f_y^2}},$ $N=\frac{f_{\mathrm{yy}}}{\sqrt{1+f_x^2+f_y^2}},$

第二类基本量的系数L、M、N与曲面的法向量曲率紧密相关。

设：

$\begin{array}{ccc}{\mathrm{\Γ}}_{11}^1=\frac{{\mathrm{GE}}_x-2{\mathrm{FF}}_x+{\mathrm{FE}}_y}{2(\mathrm{EG}-F^2)}& {\mathrm{\Γ}}_{12}^1=\frac{{\mathrm{GE}}_y-{\mathrm{FG}}_x}{2(\mathrm{EG}-F^2)}& {\mathrm{\Γ}}_{22}^1=\frac{{2\mathrm{GF}}_y-{\mathrm{GG}}_x-{\mathrm{FG}}_y}{2(\mathrm{EG}-F^2)}\\ {\mathrm{\Γ}}_{11}^2=\frac{{2\mathrm{EF}}_x-{\mathrm{EE}}_y-{\mathrm{FE}}_x}{2(\mathrm{EG}-F^2)}& {\mathrm{\Γ}}_{12}^2=\frac{{\mathrm{EG}}_x-{\mathrm{FE}}_y}{2(\mathrm{EG}-F^2)}& {\mathrm{\Γ}}_{22}^2=\frac{{\mathrm{EG}}_y-{2\mathrm{FF}}_y+{\mathrm{FG}}_x}{2(\mathrm{EG}-F^2)}\end{array}$

则Gauss方程组(Gauss Equations)可表达为：

$\left\{\begin{array}{c}{f}_{\mathrm{xx}}={\mathrm{\Γ}}_{11}^1{f}_x+{\mathrm{\Γ}}_{11}^2{f}_y+\frac{L}{\sqrt{E+G-1}}---\left(1\right)\\ f_{\mathrm{xy}}={\mathrm{\Γ}}_{12}^1{f}_x+{\mathrm{\Γ}}_{12}^2{f}_y+\frac{M}{\sqrt{E+G-1}}---\left(2\right)\\ f_{\mathrm{yy}}={\mathrm{\Γ}}_{22}^1{f}_x+{\mathrm{\Γ}}_{22}^2{f}_y+\frac{N}{\sqrt{E+G-1}}---\left(3\right)\end{array}\right.$

根据曲面论的原理，以上三个方程在地形曲面为蒙格面片形式：z＝f(x，y)上微小的“小范围”内才成立，但是为了模型实际应用的简单起见，我们假设在曲面z＝f(x，y)在足够密集的离散化程度时也成立。设{(x_i，y_j)}是对计算区域Ω进行均匀正交剖分产生的网格点，则高斯方程(Gauss Equations)组(1)，(3)方程的有限差分形式可简化为，

$\frac{f_{i+1,j}-{2f}_{i,j}+f_{i-1,j}}{h^2}={\left({\mathrm{\Γ}}_{11}^1\right)}_{i,j}\frac{f_{i+1,j}-f_{i-1,j}}{2h}+{\left({\mathrm{\Γ}}_{11}^2\right)}_{i,j}\frac{f_{i,j+1}-f_{i,j-1}}{2h}+\frac{L_{i,j}}{\sqrt{E_{i,j}+G_{i,j}-1}},(x_i,y_j)\∈\mathrm{\Ω}\backslash \∂\mathrm{\Ω}---\left(4\right)$

$\frac{f_{i,j+1}-{2f}_{i,j}+f_{i,j-1}}{h^2}={\left({\mathrm{\Γ}}_{22}^1\right)}_{i,j}\frac{f_{i,j+1}-f_{i,j-1}}{2h}+{\left({\mathrm{\Γ}}_{22}^2\right)}_{i,j}\frac{f_{i,j+1}-f_{i,j-1}}{2h}+\frac{N_{i,j}}{\sqrt{E_{i,j}+G_{i,j}-1}},(x_i,y_j)\∈\mathrm{\Ω}\backslash \∂\mathrm{\Ω}---\left(5\right)$

其中，

${\left({\mathrm{\Γ}}_{11}^1\right)}_{i,j}=\frac{G_{i,j}(E_{i+1,j}-E_{i-1,j})-{2F}_{i,j}(F_{i+1,j}-F_{i-1,j})+F_{i,j}(E_{i,j+1}-E_{i,j-1})}{4(E_{i,j}{G}_{i,j}-F_{i,j}^2)h};$

${\left({\mathrm{\Γ}}_{12}^1\right)}_{i,j}=\frac{G_{i,j}(E_{i,j+1}-E_{i,j-1})-F_{i,j}(G_{i+1,j}-G_{i-1,j})}{4(E_{i,j}{G}_{i,j}-F_{i,j}^2)h};$

${\left({\mathrm{\Γ}}_{22}^1\right)}_{i,j}=\frac{{2G}_{i,j}(F_{i,j+1}-F_{i,j-1})-G_{i,j}(G_{i+1,j}-G_{i-1,j})-F_{i,j}(G_{i,j+1}-G_{i,j-1})}{4(E_{i,j}{G}_{i,j}-F_{i,j}^2)h};$

${\left({\mathrm{\Γ}}_{11}^2\right)}_{i,j}=\frac{{2E}_{i,j}(F_{i+1,j}-F_{i-1,j})-E_{i,j}(E_{i,j+1}-E_{i,j-1})-F_{i,j}(E_{i,j+1}-E_{i,j-1})}{4(E_{i,j}{G}_{i,j}-F_{i,j}^2)h};$

${\left({\mathrm{\Γ}}_{12}^2\right)}_{i,j}=\frac{E_{i,j}(G_{i+1,j}-G_{i-1,j})-F_{i,j}(E_{i,j+1}-E_{i,j-1})}{4(E_{i,j}{G}_{i,j}-F_{i,j}^2)h};$

${\left({\mathrm{\Γ}}_{22}^2\right)}_{i,j}=\frac{E_{i,j}(G_{i,j+1}-G_{i,j-1})-{2F}_{i,j}(F_{i,j+1}-F_{i,j-1})+F_{i,j}(G_{i+1,j}-G_{i-1,j})}{4(E_{i,j}{G}_{i,j}-F_{i,j}^2)h};$

$F_{i,j}=\left(\frac{{\tilde{f}}_{i+1,j}-{\tilde{f}}_{i-1,j}}{2h}\right)\left(\frac{{\tilde{f}}_{i,j+1}-{\tilde{f}}_{i,j-1}}{2h}\right);$

$G_{i,j}=1+{\left(\frac{{\tilde{f}}_{i,j+1}-{\tilde{f}}_{i,j-1}}{2h}\right)}^2;$

$L_{i,j}=\frac{\frac{{\tilde{f}}_{i+1,j}-2{\tilde{f}}_{i,j}+{\tilde{f}}_{i-1,j}}{h}}{\sqrt{1+{\left(\frac{{\tilde{f}}_{i+1,j}-{\tilde{f}}_{i-1,j}}{2h}\right)}^2+{\left(\frac{{\tilde{f}}_{i,j+1}-{\tilde{f}}_{i,j-1}}{2h}\right)}^2}};$

$M_{i,j}=\frac{\left(\frac{{\tilde{f}}_{i+1,j+1}-{\tilde{f}}_{i+1,j-1}}{4h^2}\right)-\left(\frac{{\tilde{f}}_{i-1,j+1}-{\tilde{f}}_{j-1,j-1}}{4h^2}\right)}{\sqrt{1+{\left(\frac{{\tilde{f}}_{i+1,j}-{\tilde{f}}_{i-1,j}}{2h}\right)}^2+{\left(\frac{{\tilde{f}}_{i,j+1}-{\tilde{f}}_{i,j-1}}{2h}\right)}^2}};$

$N_{i,j}=\frac{\frac{{\tilde{f}}_{i,j+1}-2{\tilde{f}}_{i,j}+{\tilde{f}}_{i,j-1}}{h^2}}{\sqrt{1+{\left(\frac{{\tilde{f}}_{i+1,j}-{\tilde{f}}_{i-1,j}}{2h}\right)}^2{\left(\frac{{\tilde{f}}_{i,j+1}-{\tilde{f}}_{i,j-1}}{2h}\right)}^2}};$

其中，

h为计算步长(像元分辨率)；

$(x_i,y_j)\∈\mathrm{\Ω}\backslash \∂\mathrm{\Ω}$ 表示(x_i，y_j)为计算区域Ω的内点；

设 $f_{i,j}=g_{i,j},(x_i,y_j)\∈\∂\mathrm{\Ω}$ 为边界条件，简写成矩阵形式为CF＝D；

$f_{i,j}={\stackrel{\‾}{f}}_{i,j},(x_i,y_j)\∈\mathrm{\Φ}\⋐\mathrm{\Ω},$ Φ表示采样点集；

表示在点(x_i，y_j)上，利用采样点的数据(x_i，y_j，f_i，j)∈Φ进行插值而得到的f的近似值。

矩阵表达形式为：

$\left\{\begin{array}{c}{A}_{1}F=B_1\\ A_{2}F=B_2\end{array}\right.---\left(6\right)$

其中A₁和B₁分别为方程组(1)中前三个方程的系数矩阵和右端常数项矩阵；A₂和B₂分别为方程组(3)中最后一个方程的系数矩阵和右端常数项矩阵；假定计算区域Ω在x与y方向的网格点数相等，F＝[f_1，1，f_1，2，f_1，3，…，f_1，N，f_2，1，f_2，2，f_2，3…，f_2，N，…，f_N，1，f_N，2，f_N，3，…，f_N，N]^T。

这样，HASM的数值模拟问题转化为等式约束的最小二乘问题，

$\left\{\begin{array}{c}\min \left|\right|\mathrm{AF}-B{\left|\right|}_2\\ s.t.\mathrm{CF}=D\end{array}\right.---\left(7\right)$

对充分大的λ，HASM可以转化为解无约束的最小二乘问题，

$\underset{F}{\min }{||\begin{array}{c}\left[\begin{array}{c}A\\ C\end{array}\right]F-\left[\begin{array}{c}B\\ D\end{array}\right]\end{array}||}_2$

也就是求解，

$\left[\begin{array}{cc}{A}^T& \mathrm{\λ}{C}^T\end{array}\right]\left[\begin{array}{c}A\\ \mathrm{\λC}\end{array}\right]F=\left[\begin{array}{cc}{A}^T& {\mathrm{\λC}}^T\end{array}\right]\left[\begin{array}{c}B\\ \mathrm{\λD}\end{array}\right]---\left(9\right)$

设 $S=\left[\begin{array}{cc}{A}^T& {\mathrm{\λC}}^T\end{array}\right]\left[\begin{array}{c}A\\ \mathrm{\λC}\end{array}\right],$ $T=\left[\begin{array}{cc}{A}^T& {\mathrm{\λC}}^T\end{array}\right]\left[\begin{array}{c}B\\ \mathrm{\λD}\end{array}\right],$ 则HASM最终转化为

求解线性代数方程组，

Sf＝T(10)

对f进行数值模拟时，先根据采样点的值求解出基本量E、F、G、L、M、N，再利用基本量E、F、G、L、M、N求出曲面的第二类克里斯托弗尔变量(Christoffel)Γ_i ^j(i＝11，12，22；j＝1，2)。在给定的边界条件下，就可以模拟得到曲面的数值解，整个SMTS的实现流程图如图1。

根据以上原理我们已完成了高精度曲面建模(High Accuracy SurfaceModeling，HASM)方法的雏形，而且数值实验结果表明，高精度曲面建模HASM较经典模型的精度提高了多个数量级，但是HASM的发展尚处于理论实验的初级阶段，存在着计算量大、CPU时间长和运算速度慢等理论问题。因此，在理论方面还需要对HASM做进一步深入研究。同时，对运用HASM建立高精度数字地面模型过程中出现的理论技术难题，需要开展有针对性的研究。

利用Visual C++6.0，从底层进行源码开发，根据曲面论的思想，先计算曲面上每个点的第一、二类基本式的系数E、F、G和L、M、N，将它们分别写成类函数的形式，建立曲面的差分微分方程，求解该方程即可得到该曲面的模拟值。这是一次迭代过程，通过多次这样的迭代就可以得到高精度的数字地面模型。

为了避免数据源误差对比较分析的影响，我们选择数学曲面f(x，y)＝2sin(πx)sin(πy)+1作为标准曲面(如图2)，进行数值模拟分析实验。模拟结果表明，当标准曲面顶部没有采样点时，经典模型(Kriging，IDW)的模拟结果都会出现程度不同的峰值削平和边界振荡现象，导致曲面整体扭曲如图3，图4所示。经过HASM迭代后，捕捉到了被削平的峰值，消除了边界振荡，使模拟曲面与标准曲面相差无几如图5所示。

图3中：图3(a)是Kriging的插值曲面图，图3(b)是Kriging的等高线图，图3(c)是HASM的曲面图，图3(d)是HASM的等高线图。

图4中：图4(a)是IDW的插值曲面图，图4(b)是IDW的等高线图，图4(c)是HASM的曲面图，图4(d)是HASM的等高线图。

图5中：图5(a)是TIN的插值曲面图，图5(b)是TIN的等高线图，图5(c)是HASM的曲面图，图5(d)是HASM的等高线图。

以上所述，仅为本发明中的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应涵盖在本发明的包含范围之内，因此，本发明的保护范围应该以权利要求书的保护范围为准。

Claims (9)

1. 一种基于曲面论基本定律对地形曲面建立数字模型的方法，步骤如下：

步骤S1：采集地形表面曲面数据；

步骤S2：对地形曲面的采样数据进行空间离散化处理，得到离散化的曲面数据Q(0)；

步骤S3：根据曲面论，计算地形曲面数据Q(K)，Q(K)基本要素为：地形曲面上每个点的第一类基本量的系数E、F、G和第二类基本量的系数L、M、N；

步骤S4：判断迭代次数K，第一次迭代时，K等于0；当K等于0时，执行步骤S6；当K不等于0时，执行步骤S5；

步骤S5：判断迭代过程中相临两次曲面数据的差异Q(K)-Q(K-1)，当差异大于某一个阈值时，执行步骤S6；当Q(K)-Q(K-1)等于小于或等于某一个阈值时，执行步骤S7；

步骤S6：根据曲面论中的高斯方程，结合地形曲面上点的数据参数，建立地形曲面的差分微分方程为：S×f＝T，求解差分微分方程，得到地形曲面的模拟值为：f＝S^-1×T，完成一次地形曲面的迭代过程，并且K＝K+1；其中S是曲面方程的系数矩阵，f是待求的曲面未知数，T为计算出的曲面常数项；

步骤S7：输出地形曲面数据Q(K)，得到高精度的数字地形模型。

2. 根据权利要求1所述的方法，其特征在于还包括，采用离散化曲面的方法建立高斯方程，并通过迭代计算的方法进行求解。

3. 根据权利要求1所述的方法，其特征在于，所述曲面由第一类基本量和第二类基本量决定。

4. 根据权利要求1所述的方法，其特征在于，所述地形曲面为蒙格面片形式：

z＝f(x，y)，其中：x，y分别表示x坐标和y坐标。

5. 根据权利要求1所述的方法，其特征在于，所述地形曲面的第一类基本量的系数为：

$E=1+f_x^2$ F＝f_xf_y $G=1+f_y^2$ ，第一类基本量的系数E表示当变化x时，变化将离曲面的多远，G表示当变化y时，变化将离曲面多远，F表示当变化x，y同时变化时，他们之间的夹角是多大。

6. 根据权利要求1、4所述的方法，其特征在于，所述地形曲面的第二类基本量的系数为：

$L=\frac{f_{\mathrm{xx}}}{\sqrt{1+f_x^2+f_y^2}},$

$M=\frac{f_{\mathrm{xy}}}{\sqrt{1+f_x^2+f_y^2}},$

$N=\frac{f_{\mathrm{yy}}}{\sqrt{1+f_x^2+f_y^2}},$

第二类基本量的系数L、M、N与曲面的法向量曲率紧密相关。

7. 根据权利要求4所述的方法，其特征在于还包括，选择地形数学曲面f(x，y)＝2sin(πx)sin(πy)+1作为标准曲面，对标准曲面进行数值模拟分析，当标准曲面顶部没有采样点时，经典模型的模拟结果有峰值削平和边界振荡，曲面整体扭曲；将曲面整体扭曲经过HASM迭代后，获得被削平的峰值，消除边界振荡，使模拟曲面与标准曲面相同。

8. 根据权利要求1所述的方法，其特征在于还包括，对地形数学曲面f进行数值模拟时，先根据采样点的值求解出基本量E、F、G、L、M、N，再利用基本量E、F、G、L、M、N求出曲面的第二类克里斯托弗尔变量Γ_i ^j，i＝11，12，22；j＝1，2，在给定的边界条件下，模拟得到曲面的数值解。

9. 权利要求1所述的方法用于空间图像数据处理。

Images