CN101900546A - 对地球表面地形地貌离散表达的数字高程模型的构建方法 - Google Patents

Abstract

本发明为对地球表面地形地貌离散表达的数字高程模型的构建方法，首先对地球表面原始数据提取采样点的高程值及空间离散化处理，生成多个已知采样点高程值的地形曲面网格点；求算出未采样点的地形曲面网格点高程值；根据所有地形曲面网格点高程值计算每个地形曲面网格点的第一和第二类基本量，用高斯方程和高程采样点建立地形曲面方程组；利用高斯-赛德尔迭代算法解算地形曲面方程组，每次迭代结束得出数字高程模型的模拟结果；判断数字高程模型模拟结果是否满足精度要求，如果满足要求，停止迭代计算，输出一组有序数值阵列形式表示地面高程实体的数字高程模型用于对地球表面地形地貌进行离散表达；如果不满足要求，重复迭代计算。

Description

对地球表面地形地貌离散表达的数字高程模型的构建方法

技术领域

本发明属于地球观测与导航技术领域，涉及地理信息系统曲面生成方法，尤其是数字高程模型的构建方法。

背景技术

数字高程模型(Digital Elevation Model，DEM)是对地球表面地形地貌的离散表达，是GIS进行各种空间分析的数据基础，广泛应用于地貌定量分析、流域分析、土壤侵蚀分析、遥感图像辅助分类、道路设计、军事演练等领域，DEM作为数字地形模拟的重要成果已经成为国家空间数据基础设施的基本内容之一，并被纳入数字化空间数据框架进行规模化生产。

自1958年美国麻生理工学院Miller教授为公路设计首次提出DEM以来，DEM已经取得了长足的进展。由于DEM的广泛应用，对DEM精度的研究也受到有关专家的高度重视。在近30年的发展中，产生了许多关于DEM精度评价指标，如中误差、标准差、误差直方图、空间相关系数等，并提出了许多改进DEM精度的方法，这主要表现在以下几个方面：

(1)原始数据精度的评价分析

由于DEM是对地球表面地形地貌的一种离散的近似数字表达，因此不可避免的含有地形表达误差、人为采样误差等。原始数据误差决定着DEM精度，因此在进行各种插值之前，应分析采样点误差，剔除误差较大点，充分挖掘原始采样数据的有用信息，提高插值精度。如等高线数据用于DEM构建时，如果仅考虑等高线提供的离散数据点，忽略对应区域山脊点(线)、山谷点(线)等地形特征点(线)这些隐藏的信息，会导致构建的DEM与原始地形不符。

(2)各种内插算法精度的对比评价

各种内插算法都是为不同的应用目的而提出的，因此在进行插值前应根据采样点的空间分布详细对比分析各种算法的使用条件、最优参数、插值精度等，确定最佳插值方法或者对原始插值算法改进，提高插值精度。如反距离权重插值(IDW)没有考虑插值点周围采样点的空间分布、密度，而将模型参数设成定值导致插值结果出现“牛眼”等异常现象。

(3)专业数据的综合利用

DEM一直是水文分析的重要数据来源，高精度DEM可以提高各种水文分析的精度。如Hutchinson以粗糙度罚函数(roughness penalty)为基本方程^[1]Hutchinson M F.Optimising the degree of data smoothing for locally adaptive finite element bivariate smoothing splines.Australian & New Zealand Industrial and Applied Mathematics Journal，2000，42：744-796，引入了地形强化算法(drainage enforcement algorithm)生成的DEM极大提高了流域面积、水系位置、水系长度等计算结果精度。由于充分利用了流域水系矢量数据，并对计算数据进行了专业处理，因此Hutchinson研发的ANUDEM已成为世界上流行的建立水文、地貌关系正确的DEM专业化软件。但需要注意的是，尽管对DEM的改进提高了水文模型的精度，但是改进的DEM同时改变了原始DEM的高程和对应地形参数。

现有的空间插值模型获取内插点数值时往往忽略采样点的空间相关性，使得已知点之间的关系遭到破坏，最终导致较大的误差甚至出现各种错误。在以往的研究中，Evans虽然认识到坡面、坡向和曲率是反映局部地形规律的重要变量^[2]Evans I S.An Integrated System of Terrain Analysis & Slope Mapping.Dept.of Geography，University of Durham，1979，但他并没有将它们用于表达数字地面模型。事实上根据微分几何学原理，坡面、坡向和曲率只是地面剖面线的决定要素，数字地面模型的决定要素是曲面的第一类基本量和第二类基本量。

发明内容

为了克服现有技术的不足，且从本质上解决曲面建模的误差问题，本发明的目的是提供一种能对地球表面地形地貌离散表达的数字高程模型的构建方法。

为达成所述目的，本发明提供对地球表面地形地貌离散表达的数字高程模型的构建方法，该方法的实现步骤包括：

步骤S1：利用数据采集设备对地球表面原始数据进行高程采样，获得地球表面多个采样点的高程值；对多个采样点的高程值进行空间离散化处理，得到多个离散化的地形曲面网格点；将每个采样点的高程值赋给离其最近的离散化的地形曲面网格点，生成多个已知采样点高程值的地形曲面网格点；

步骤S2：利用插值方法，根据已知采样点高程值的地形曲面网格点对未采样点的地形曲面网格点进行空间预插值，求算出未采样点的地形曲面网格点高程值；

步骤S3：根据所有地形曲面网格点高程值计算每个地形曲面网格点的第一类基本量E，F，G和第二类基本量L，M，N，其中第一类基本量表示地形曲面上曲线的弧长、地形曲面上两个方向的夹角和地形曲面域的面积，第二类基本量刻画了地形曲面空间中的弯曲性；用高斯方程作为基本方程和高程采样点作为约束方程建立地形曲面方程组SUⁿ⁺¹＝Tⁿ，其中，S为地形曲面方程组的系数矩阵，U为待求的数字高程模型数据，T为计算出的地形曲面常数项，n为地形曲面方程组迭代次数；

步骤S4：利用高斯-赛德尔迭代算法解算地形曲面方程组SUⁿ⁺¹＝Tⁿ，整个解算过程为迭代计算，每次迭代结束得出数字高程模型的模拟结果；

步骤S5：判断数字高程模型模拟结果是否满足精度要求，如果满足要求，执行步骤6；如果不满足要求，执行步骤3；

步骤S6：停止迭代计算，输出一组有序数值阵列形式表示地面高程实体的数字高程模型，该数字高程模型用于对地球表面地形地貌进行离散表达。

优选地，设所述地形曲面z的表达式为：z＝u(x，y)，其中x，y分别表示地形曲面点的横、纵坐标；函数u表示地形曲面点横、纵坐标(x，y)和地形曲面高程点z的映射关系。

优选地，所述第一类基本量E，F，G和第二类基本量L，M，N的表达式为：

$E=1+u_x^2,F=u_x{u}_y,G=1+u_y^2$

$L=\frac{u_{\mathrm{xx}}}{\sqrt{1+u_x^2+u_y^2}},M=\frac{u_{\mathrm{xy}}}{\sqrt{1+u_x^2+u_y^2}},N=\frac{u_{\mathrm{yy}}}{\sqrt{1+u_x^2+u_y^2}},$ 其中：u_x，u_y分别表示函数u对地形曲面点的横、纵坐标x，y求一阶偏导数；u_xy，u_xx，u_yy分别表示函数u对地形曲面点的横、纵坐标x，y求二阶偏导数。

优选地，所述地形曲面方程组SUⁿ⁺¹＝Tⁿ为大型稀疏线性方程组，且该方程组系数矩阵S每行最多有九个非零元素。

优选地，所述解算地形曲面方程组所需计算的时间与地形曲面网格点总数呈线性关系。

为达成所述目的，本发明的第二方面，是提供一种地球表面地形地貌离散表达的数字高程模型用作水文模型、地貌模型的数据基础。

本发明的有益效果：本发明以GS解算地形曲面方程组(该过程简称为：HASM-GS)的计算时间与地形曲面网格点总数呈线性相关，比采用直接法计算时间降低两个数量级。本发明生成的数字高程模型精度要比传统的方法(SPLINE，IDW和KRIGING)高多倍。本发明可作为水文模型、地貌模型的数据基础。

选择高斯合成曲面作为研究对象如图2示出，将HASM-GS与MATLAB提供的经典方程组解算方法比较表明(表1)，各种方法在达到相同的精度2.65e-006时，HASM-GS耗时远远小于MATLAB提供的各种函数。除去HASM-GS，耗时最少的为二次共轭梯度法(CGS)，该法耗时约为HASM-GS的12.409倍；耗时最多的为最小二乘QR分解法(LSQR)，该法耗时约为HASM-GS的470.912倍。根据不同的模拟区域验证HASM-GS的计算时间表明如表2和图3示出，HASM-GS计算时间与地形曲面网格点总数成线性相关，HASM-GS极大的提高了数字高程模型的生成速度。由表2和图3可见，HASM-GS计算时间与网格数呈非常好的线性相关关系，其回归方程为：

t＝2.036×10^-5gn+0.010(R²＝1)     (1)

式(1)中，t为计算时间，gn为网格数目。

表1均匀采样条件下各种方法计算效率比较

表2HASM-GS计算时间与模拟区域网格数比较

附图说明

图1本发明流程图；

图2高斯合成曲面；

图3HASM-GS计算时间和地形曲面网格点总数的回归曲线；

图4(a)-(d)各种方法的检核点离散图及对应的相关系数；

具体实施方式

下面结合附图详细说明本发明技术方案中所涉及的各个细节问题。

本发明以曲面论为理论基础建立了对地球表面地形地貌离散表达的高精度数字高程模型(HASM)。经典的曲面论是本发明的出发点，也是HASM区别于以往数字高程模型生成方法的不同之处。基于曲面论基本定律建立的数字高程模型理论更加成熟，使其精度远远高于其它方法。本发明是以高斯(Gauss)方程作为基本方程，以高程采样点作为约束条件，将数字高程模型模拟转换为解算地形曲面方程组SUⁿ⁺¹＝Tⁿ，该地形曲面方程组为大型稀疏线性方程组，且该方程组系数矩阵S每行最多有九个非零元素。求解该大型稀疏线性方程组即可得到高精度的数字高程模型。请参见图1本发明方法的流程图，具体实施如下：

利用数据采集设备(如全站仪、水准仪、全球定位系统(GPS)等)对地球表面原始数据进行高程采样，采集地球表面多个采样点的高程值，其中，高程数据源主要有地形图、航天航空遥感图像、雷达、野外实测等；将多个采样点的高程值进行空间离散化处理，得到离散化的地形曲面网格点；将每个采样点的高程值赋给离其最近的地形曲面网格点，这些网格点即为已知采样点高程值的地形曲面网格点；利用插值算法(如线性插值、三次样条插值等)借助已知知采样点高程值的地形曲面网格点对未知采样点的地形曲面网格点空间预插值，求算未采样点的地形曲面网格点的高程值。利用所有地形曲面网格点的高程值计算所得地形曲面网格点的第一类基本量E，F，G和第二类基本量L，M，N，其中第一类基本量可以表示地形曲面上曲线的弧长、地形曲面上两个方向的夹角和地形曲面(域)的面积，第二类基本量刻画了地形曲面空间中的弯曲性。以高斯(Gauss)方程作为基本方程(请详见公式2)，以高程采样点作为约束方程，建立对地球表面地形地貌离散表达的高精度数字高程模型地形曲面方程组：SUⁿ⁺¹＝Tⁿ，该方程组为大型稀疏线性方程组，且方程组系数矩阵S每行最多有九个非零元素，U为待求数字高程模型数据，T为计算出的地形曲面常数项，n为方程组迭代次数；解算方程组称为内迭代，更新T称为外迭代。利用高斯-赛德尔(Gauss-Seidel)解算地形曲面方程组，SUⁿ⁺¹＝Tⁿ，整个解算过程为迭代计算，每次迭代结束即可得出高精度数字高程模型的模拟结果。

设地形曲面表达式为z＝u(x，y)，其中x，y分别表示地形曲面点的横、纵坐标；函数u表示地形曲面点横、纵坐标(x，y)和地形曲面高程点z的映射关系，则Gauss方程的有限差分形式为：

$\left\{\begin{array}{c}(u_{i+1,j}-2u_{i,j}+u_{i-1,j})/h^2={\left({\mathrm{\Γ}}_{11}^1\right)}_{i,j}(u_{i+1,j}-u_{i-1,j})/\left(2h\right)+{\left({\mathrm{\Γ}}_{11}^2\right)}_{i,j}(u_{i,j+1}-u_{i,j-1})/\left(2h\right)+\\ L_{i,j}/{(E_{i,j}+G_{i,j}-1)}^{1/2},(x_i,y_i)\∈\mathrm{\Ω}\backslash \∂\mathrm{\Ω}\\ (u_{i,j+1}-2u_{i,j}+u_{i,j-1})/h^2={\left({\mathrm{\Γ}}_{22}^1\right)}_{i,j}(u_{i+1,j}-u_{i-1,j})/\left(2h\right)+{\left({\mathrm{\Γ}}_{22}^2\right)}_{i,j}(u_{i,j+1}-u_{i,j-1})/\left(2h\right)+\\ N_{i,j}/{(E_{i,j}+G_{i,j}-1)}^{1/2},(x_i,y_j)\∈\mathrm{\Ω}\backslash \∂\mathrm{\Ω}\end{array}\right.---\left(2\right)$

式中，u_i，j表示地形曲面第i行j列点高程值；i、j为自然数，

为地形曲面的第二类克里斯托费尔符号，x_i，y_j表示地形曲面点的平面坐标；

表示所有地形曲面网格点的集合；h表示网格分辨率；E_i，j、F_i，j、G_i，j、L_i，j、N_i，j表示地形曲面网格点的第一类基本量和第二类基本量，它们的有限差分形式为：

$E_{i,j}\≈1+{[({\tilde{u}}_{i+1,j}-{\tilde{u}}_{i-1,j})/\left(2h\right)]}^2---\left(3\right)$

$F_{i,j}\≈({\tilde{u}}_{i+1,j}-{\tilde{u}}_{i-1,j})({\tilde{u}}_{i,j+1}-{\tilde{u}}_{i,j-1})/\left(4h^2\right)---\left(4\right)$

$G_{i,j}\≈1+{[({\tilde{u}}_{i,j+1}-{\tilde{u}}_{i,j-1})/2h]}^2---\left(5\right)$

$L_{i,j}\≈({\tilde{u}}_{i+1,j}-2{\tilde{u}}_{i,j}+{\tilde{u}}_{i-1,j})/\left\{h^2{\{1+{[({\tilde{u}}_{i+1,j}-{\tilde{u}}_{i-1,j})/\left(2h\right)]}^2+[{({\tilde{u}}_{i,j+1}-{\tilde{u}}_{i,j-1})/\left(2h\right)]}^2\}}^{1/2}\right\}---\left(6\right)$

$N_{i,j}\≈({\tilde{u}}_{i,j+1}-2{\tilde{u}}_{i,j}+{\tilde{u}}_{i,j-1})/\left\{h^2{\{1+{[({\tilde{u}}_{i+1,j}-{\tilde{u}}_{i-1,j})/\left(2h\right)]}^2+[{({\tilde{u}}_{i,j+1}-{\tilde{u}}_{i,j-1})/\left(2h\right)]}^2\}}^{1/2}\right\}---\left(7\right)$

式中，

表示地形曲面网格点的高程值。地形曲面的第二类克里斯托费尔符号Γ_k ^l(k＝11，22；l＝1，2)可表示为：

${\left({\mathrm{\Γ}}_{11}^1\right)}_{i,j}\≈{[G}_{i,j}(E_{i+1,j}-E_{i-1,j})-2U_{i,j}(U_{i+1,j}-U_{i-1,j})+U_{i,j}(E_{i,j+1}-E_{i,j-1})\left]\right[4h(E_{i,j}{G}_{i,j}-U_{i,j}^2)]---\left(8\right)$

${\left({\mathrm{\Γ}}_{11}^2\right)}_{i,j}\≈[2E_{i,j}(U_{i+1,j}-U_{i-1,j})-E_{i,j}(E_{i,j+1}-E_{i,j-1})-U_{i,j}(E_{i,j+1}-E_{i,j-1})]/\left[4h(E_{i,j}{G}_{i,j}-U_{i,j}^2)\right]---\left(9\right)$

${\left({\mathrm{\Γ}}_{22}^1\right)}_{i,j}\≈[2G_{i,j}(U_{i,j+1}-U_{i,j-1})-G_{i,j}(G_{i+1,j}-G_{i-1,j})-U_{i,j}(G_{i,j+1}-G_{i,j-1})]/\left[4h(E_{i,j}{G}_{i,j}-U_{i,j}^2)\right]---\left(10\right)$

${\left({\mathrm{\Γ}}_{22}^2\right)}_{i,j}\≈[E_{i,j}(G_{i,j+1}-G_{i,j-1})-2U_{i,j}(U_{i,j+1}-U_{i,j-1})+U_{i,j}(G_{i+1,j}-G_{i-1,j})]/\left[4h(E_{i,j}{G}_{i,j}-U_{i,j}^2)\right]---\left(11\right)$

高斯(Gauss)方程的矩阵表达式为：

$\left\{\begin{array}{c}A{U}^{n+1}=U^n\\ B{U}^{n+1}={\mathrm{EE}}^n\end{array}---\left(12\right)\right.$

式中，

L，C分别表示矩阵的行列数；A，B分别表示方程(2)左端系数矩阵；Dⁿ，EEⁿ分别表示方程(2)右端向量。

在高程采样点处的模拟值等于已知高程采样值，因此高精度数字高程模型(HASM)数值模拟转换为等式约束的最小二乘问题，

$\left\{\begin{array}{c}{\min \left|\right|P{U}^{n+1}-Q^n\left|\right|}_2\\ s.t.M{U}^{n+1}=H\end{array}---\left(13\right)\right.$

式中，

P，Qⁿ均为矩阵集合；s.t.表示等式约束条件；M表示高程采样点矩阵，H表示采样点高程值矩阵；当高程采样点在(i，j)时，M(i，(i-1)×C+j)＝1，H(i)＝u_i，j。

对于充分大的采样点权重因子λ，数字高程模型求算转换为解算无约束的最小二乘问题，即HASM表达式为：

$\underset{F}{\min }{||\begin{array}{c}\end{array}\left[\begin{array}{c}P\\ \mathrm{\λM}\end{array}\right]U^{n+1}-\left[\begin{array}{c}{Q}^n\\ \mathrm{\λH}\end{array}\right]||}_2---\left(14\right)$

也就是求解

SUⁿ⁺¹＝Tⁿ           (15)

式中，其中解算方程组(15)称为内迭代，更新右端向量Tⁿ称为外迭代。

利用Gauss-Seidel(GS)分量形式解算方程组(15)可表示为：

for k＝0，1，2，…

  y₁＝u_k

  for i＝1，2，…，n

   y_i+1＝y_i+α_ie_i    (16)

 end for i

 u_k+1＝y_n+1

 stopping criteria

end for k

式中：u_k表示数字高程模型第k次模拟值，k为0时，u_k用u₀来表示数字高程模型初始模拟值；其余各个符号的表达式为：e₁＝(1，0，0，…，0)，e₂＝(0，1，0，…，0)，…，e_n＝(0，0，0，…，1)；α_i＝-p_i/s_ii，P_i＝(S_i，y_i)-T_i；e_i表示单位向量；α_i表示数字高程模型每次迭代更新项；p_i为临时变量；S_ii表示HASM系数矩阵第i行，i列元素；S_i表示HASM的系数矩阵第i行向量；T_i表示计算出的地形曲面常数项。公式(16)本发明称为高精度数字高程模型-高斯赛德尔迭代(HASM-GS)。

利用计算机和MATLAB语言，根据曲面论的思想，先计算地形曲面上每个网格点的第一类基本量E，F，G和第二类基本量L，M，N，然后根据曲面的差分形式建立曲面方程，利用GS迭代法即(16)求解该方程即可得到高精度数字高程模型。

实例选自山东省济南市某测区，该测区面积为3.7km²，测区平均高程为180.5m，高程标准差为60.5m。为了提高高程采样点的精度，首先在该测区均匀布设GPS控制点321个，然后利用全站仪采集9756个高程点，高程采样点高程精度为5cm，平面精度为1cm。高程采样点的平均密度为2637points/km²，对应的高程采样点的间隔(ES)为19.5m。随机选择9260个点用于模拟DEM，其余496个点用于精度检验。除了HASM-GS外，传统的插值方法(SPLINE，KRIGING和IDW)也用于构建DEM，并将其精度与HASM-GS进行比较。精度指标除了采用中误差(RMSE)外，平均绝对误差(MAE)也用于精度评价，平均绝对误差(MAE)的计算公式为：

$\mathrm{MAE}=\frac{\underset{i=1}{\overset{N}{\mathrm{\Σ}}}|E{V}_i-S{V}_i|}{N}---\left(17\right)$

式中，EV_i表示高程观测值；SV_i表示数字高程模型模拟值；N为检核的高程点数目。各种方法的模拟结果如表3所示。

表3各种方法的插值精度比较(比值：经典方法模拟误差与HASM-GS比值)

由表3可见，无论是MAE还是RMSE作为精度指标，HASM-GS插值精度要远高于传统的方法，SPLINE，KRIGING，IDW的MAE分别是HASM-GS的9.0125，8.5890和8.3383倍；RMSE分别是HASM-GS的7.0611，6.4628和8.6510倍。

为了更准确的描述每个检核点的模拟值与观测值之间的差异，本发明提供了检核点的离散图请见图4(a)-图4(d)，其中，横轴表示观测值(EV)，纵轴表示模拟值(SV)，因此当观测值等于模拟值时，检核点落在y＝x这条直线上。由图可见，除了HASM-GS外，每种方法均存在模拟值与观测值相差较大的检核点。利用SPSS软件计算各种方法的模拟值和观测值的拟合曲线相关系数得，HASM-GS的相关系数最大，表明精度最高(图4a)；而IDW相关系数最小，表明精度最低(图4b)。图4(a)HASM-GS；相关系数：0.9999943；图4(b)IDW；相关系数：0.9995738；图4(c)KRIGING；相关系数：0.9997633；图4(d)SPLINE；相关系数：0.9997172。

应指出的是，所描述的实施例仅旨在便于对本发明的理解，而对其不起任何限定作用。但本发明的保护范围并不局限于此，任何熟悉该技术的人在本发明所揭露的技术范围内，可理解想到的变换或替换，都应属于本发明的范围。

Claims (6)

1.一种对地球表面地形地貌离散表达的数字高程模型的构建方法，其特征在于，该方法的实现步骤如下：

步骤S1：利用数据采集设备对地球表面原始数据进行高程采样，获得地球表面多个采样点的高程值；对多个采样点的高程值进行空间离散化处理，得到多个离散化的地形曲面网格点；将每个采样点的高程值赋给离其最近的离散化的地形曲面网格点，生成多个已知采样点高程值的地形曲面网格点；

步骤S2：利用插值方法，根据已知采样点高程值的地形曲面网格点对未采样点的地形曲面网格点进行空间预插值，求算出未采样点的地形曲面网格点高程值；

步骤S3：根据所有地形曲面网格点高程值计算每个地形曲面网格点的第一类基本量E，F，G和第二类基本量L，M，N，其中第一类基本量表示地形曲面上曲线的弧长、地形曲面上两个方向的夹角和地形曲面域的面积，第二类基本量刻画了地形曲面空间中的弯曲性；用高斯方程作为基本方程和高程采样点作为约束方程建立地形曲面方程组SUⁿ⁺¹＝Tⁿ，其中，S为地形曲面方程组的系数矩阵，U为待求的数字高程模型数据，T为计算出的地形曲面常数项，n为地形曲面方程组迭代次数；

步骤S4：利用高斯-赛德尔迭代算法解算地形曲面方程组SUⁿ⁺¹＝Tⁿ，整个解算过程为迭代计算，每次迭代结束得出数字高程模型的模拟结果；

步骤S5：判断数字高程模型模拟结果是否满足精度要求，如果满足要求，执行步骤6；如果不满足要求，执行步骤3；

步骤S6：停止迭代计算，输出一组有序数值阵列形式表示地面高程实体的数字高程模型，该数字高程模型用于对地球表面地形地貌进行离散表达。

2.根据权利要求1所述的数字高程模型的构建方法，其特征在于，设所述地形曲面z的表达式为：z＝u(x，y)，其中x，y分别表示地形曲面点的横、纵坐标；函数u表示地形曲面点横、纵坐标(x，y)和地形曲面高程点z的映射关系。

3.根据权利要求1所述的数字高程模型的构建方法，其特征在于，所述第一类基本量E，F，G和第二类基本量L，M，N的表达式为：

$E=1+u_x^2,F=u_x{u}_y,G=1+u_y^2$

$L=\frac{u_{\mathrm{xx}}}{\sqrt{1+u_x^2+u_y^2}},M=\frac{u_{\mathrm{xy}}}{\sqrt{1+u_x^2+u_y^2}},N=\frac{u_{\mathrm{yy}}}{\sqrt{1+u_x^2+u_y^2}},$ 其中：n_x，u_y分别表示函数u对地形曲面点的横、纵坐标x，y求一阶偏导数；u_xy，u_xx，u_yy分别表示函数u对地形曲面点的横、纵坐标x，y求二阶偏导数。

4.根据权利要求1所述的数字高程模型的构建方法，其特征在于，所述地形曲面方程组SUⁿ⁺¹＝Tⁿ为大型稀疏线性方程组，且该方程组系数矩阵S每行最多有九个非零元素。

5.根据权利要求1所述的数字高程模型的构建方法，其特征在于，所述解算地形曲面方程组所需计算时间与地形曲面网格点总数呈线性关系。

6.根据权利要求1所述的对地球表面地形地貌离散表达的数字高程模型用作水文模型、地貌模型的数据基础。

Images