CN102054294A - 基于曲面论和优化控制理论的曲面建模方法 - Google Patents

Abstract

本发明提供一种基于曲面论和优化控制理论的曲面建模方法，包括：采集等高线高程值和已知离散点高程值；根据所述等高线高程值和已知离散点高程值计算原始模拟区域中未采样网格点的高程值上下界，得到所述原始模拟区域中所有网格点的高程值上下界；用所述原始模拟区域中所有网格点的高程值上下界作为约束方程建立第一不等式约束方程组；根据得到的网格点的高程值计算所述原始模拟区域中每个网格点的第一类基本量和第二类基本量，并结合高斯方程以及第一不等式约束方程组建立地形曲面的差分微分方程；利用迭代算法解算该差分微分方程，当相邻两次曲面数据的差值满足精度要求时，停止迭代，输出地形曲面数据，从而得到高精度的数字地形模型。

Description

基于曲面论和优化控制理论的曲面建模方法

技术领域

本发明涉及地理信息系统中所使用的一种曲面建模方法，特别涉及一种基于曲面论和优化控制理论的曲面建模方法。

背景技术

DTM(Digital Terrain Model，数字地面模型)是地形表面形态属性信息的数字表达，是带有空间位置特征和地形属性特征的数字描述。DTM应用广泛，在测绘领域中被用于绘制等高线、坡度坡向图、立体透视图等。在遥感应用中可作为分类的辅助数据。另外，DTM还是地理信息系统的基础数据，可用于土地利用现状的分析、合理规划及洪水险情预报等。

由于DTM的误差会通过建模过程被传播，并体现在最终成果中。地理信息系统中各种误差的累积使最终成果的可用性令人怀疑，尤其是地理信息系统最常用的叠合操作对固有误差和运算误差都非常敏感。因此，对DTM的研究主要为如何有效减少DTM的误差问题。常见的DTM的误差包括：采样误差、采样仪器产生的误差、变换地面控制点产生的误差、数学模型带来的误差、数据源的传播误差、表达误差和栅格分辨率转换引起的误差等。

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

DEM(Digital Elevation Model，数字高程模型)是建立DTM的基础数据，其它的地形要素(例如：坡度、坡向等)可以通过DEM直接或间接导出。因此，DEM的精度和构建DEM的方法是影响DTM精度的主要因素。

因此，寻找一种有效的建立DEM的方法是提高DTM精度的关键任务。

发明内容

本发明提供一种基于曲面论和优化控制理论的曲面建模方法，主要用于建立高精度的DTM。

本发明提供的技术方案如下：

本发明提供一种基于曲面论和优化控制理论的曲面建模方法，包括以下步骤：

(1)采集等高线高程值和已知离散点高程值，并设置迭代参数初始值，所述迭代参数包括：迭代次数k和容差delta；

(2)将等高线所在的原始模拟区域空间离散化为网格点形式，然后根据所述等高线高程值和已知离散点高程值计算所述原始模拟区域中未采样网格点的高程值上下界，得到所述原始模拟区域中所有网格点的高程值上下界；

(3)用所述原始模拟区域中所有网格点的高程值上下界作为约束方程建立第一不等式约束方程组；

(4)根据得到的网格点的高程值计算所述原始模拟区域中每个网格点的第一类基本量E、F、G和第二类基本量L、M、N；其中，所述第一类基本量用于表示地形曲面上曲线的弧长、地形曲面上两个方向的夹角和地形曲面域的面积；所述第二类基本量用于刻画地形曲面空间中的弯曲性；并结合曲面论中的高斯方程以及所述第一不等式约束方程组建立地形曲面的差分微分方程：S*f^k+1＝t^k；其中，S是曲面方程的系数矩阵，f是待求的数字高程模型数据，上标k表示第k次迭代，f^k+1表示第k次迭代得到的地形曲面，t为计算出的地形曲面常数项；

(5)利用迭代算法解算所述S*f^k+1＝t^k，每次迭代结束均得出数字高程模型的模拟结果；其中，迭代运算的初始值S(0)(k＝0)为预先设定值；

(6)在迭代过程中，需要判断相临两次曲面数据的差值S(k)-S(k-1)(k＞＝1)是否小于所述容差，如果判断结果为否，，则重复执行步骤(4)(5)中的迭代过程；如果判断结果为是，则终止迭代运算，执行步骤(7)；

(7)输出地形曲面数据S(k)，得到高精度的数字地形模型。

优选的，所述根据所述等高线高程值和已知离散点高程值计算所述原始模拟区域中未采样网格点的高程值上下界，得到所述原始模拟区域中所有网格点的高程值上下界包括：

直接在网格点形式的基础上并根据所述等高线高程值和已知离散点高程值计算所述原始模拟区域中未采样网格点的高程值上下界，得到所述原始模拟区域中所有网格点的高程值上下界。

优选的，所述计算所述原始模拟区域中未采样网格点的高程值上下界包括以下步骤：

判断位于所述网格点中心位置处的格网中心点与最接近等高线间距离是否小于预设阈值；

如果判断结果为是，则所述网格点高程值上下界分别为所述最接近等高线的高程值加减预设微小值所得到的值；

如果判断结果为否，则包括以下情况：

当所述格网中心点被两条父子关系的等高线包围时，所述格网中心点所在的网格点的高程值上下界分别为所述两条父子关系中父等高线的高程值和子等高线的高程值；

当所述格网中心点为马鞍处的点时，所述格网中心点所在的网格点的高程值上下界分别为包围该格网中心点的所有等高线高程值的极大值和极小值；

当所述格网中心点为山峰点时，所述格网中心点所在的网格点的高程值下界为包围该格网中心点的最接近等高线的高程值；并进一步判断所述格网中心点所在的山峰是否还存在采样点，如果判断结果为是，则所述格网中心点所在的网格点的高程值上界为该采样点的高程值；如果判断结果为否，则所述格网中心点所在的网格点的高程值上界为包围该格网中心点的最接近等高线的高程值加上1倍的等高距的值；

当所述格网中心点为山谷点时，所述格网中心点所在的网格点的高程值上界为包围该格网中心点的最接近等高线的高程值；并进一步判断所述格网中心点所在的山谷是否还存在采样点，如果判断结果为是，则所述格网中心点所在的网格点的高程值下界为该采样点的高程值；如果判断结果为否，则所述格网中心点所在的网格点的高程值下界为包围该格网中心点的最接近等高线的高程值减去1倍的等高距的值。

优选的，所述根据所述等高线高程值和已知离散点高程值计算所述原始模拟区域中未采样网格点的高程值上下界，得到所述原始模拟区域中所有网格点的高程值上下界包括以下步骤：

建立所述原始模拟区域中各等高线的等高线树；

根据建立的等高线树以及所述等高线高程值和已知离散点高程值计算所述原始模拟区域中未采样网格点的高程值上下界。

优选的，在步骤(1)之后还包括：向外围扩展所述原始模拟区域，获得扩展区域，所述扩展区域包含内部区域和新增区域两部分；其中，所述内部区域范围与所述原始模拟区域范围相同；

利用外推的方法获得所述新增区域中各网格点的高程值上下界；

根据所述原始模拟区域和所述扩展区域的扩展关系，获得所述内部区域中各网格点的高程值上下界；

根据获得的新增区域中各网格点的高程值上下界和内部区域中各网格点的高程值上下界，在所述扩展区域上建立第二不等式约束方程组。

优选的，所述步骤(5)中对所述地形曲面的差分微分方程进行迭代时，还包括下述步骤：

首先进行前迭代，所述前迭代为根据所述第二不等式约束方程组在所述扩展区域上的迭代；当得到的相临两次曲面数据的差值S(k)-S(k-1)(k＞＝1)小于预设值时，转入后迭代；所述后迭代为根据所述第一不等式约束方程组在所述原始模拟区域上的迭代。

优选的，所述地形曲面为蒙格面片的形式：z＝f(x，y)，其中，x，y分别表示x坐标和y坐标。

优选的，所述地形曲面的第一类基本量的系数为：

F＝f_xf_y

其中，所述f_x、f_y分别表示函数f对地形曲面点的横、纵坐标求一阶偏导数。优选的，所述地形曲面的第二类基本量的系数为：

$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}}$

其中，f_x、f_y分别表示函数f对地形曲面点的横、纵坐标求一阶偏导数；f_xx、f_xy、f_yy分别表示函数f对地形曲面点的横、纵坐标求二阶偏导数。

本发明还提供一种空间图像数据处理方法，应用本发明所提供的基于曲面论和优化控制理论的曲面建模方法。

本发明的有益效果如下：

本发明提供的基于曲面论和优化控制理论的曲面建模方法，结合曲面论中的高斯方程以及不等式约束方程组建立地形曲面的差分微分方程，通过迭代算法解算该方程，从而得到了高精度的数字地形模型。

附图说明

图1为根据等高线确定离散网格点上下界约束示意图；

图2为本发明提供的一种基于曲面论和优化控制理论的曲面建模方法的流程示意图；

图3为本发明提供的另一种基于曲面论和优化控制理论的曲面建模方法的流程示意图；

图4为通过使用本发明所提供的HASMOC方法得到的DEM；

图5为使用传统的TIN方法得到的DEM；

图6为千烟洲扫描矢量化原始等高线图和离散点数据；

图7为通过使用本发明所提供的HASMOC方法对图6所示的等高线反演DEM得到的阴影图；

图8为对图7所得到的DEM阴影图进行等高线回放所得到的等高线回放图。

具体实施方式

以下结合附图详细介绍本发明。

一、等式约束的HASM算法

设{(x_i，y_i)|x_i＝i*h，y_j＝j*h，0≤i≤I+1，0≤j≤J+1}是计算区域Ω进行均匀正交剖分的网格，则HASM主方程Dirichlet边界问题的有限差分迭代形式可表达为，

$\left\{\begin{array}{c}\frac{f_{i+1,j}^{n+1}-2f_{i,j}^{n+1}+f_{i-1,j}^{n+1}}{h^2}={\left({\mathrm{\Γ}}_{11}^1\right)}_{i,j}^n\frac{f_{i+1,j}^n-f_{i-1,j}^n}{2h}+{\left({\mathrm{\Γ}}_{11}^2\right)}_{i,j}^n\frac{f_{i,j+1}^n-f_{i,j-1}^n}{2h}+\frac{L_{i,j}^n}{\sqrt{E_{i,j}^n+G_{i,j}^n-1}}\\ \frac{f_{i,j+1}^{n+1}-2f_{i,j}^{n+1}+f_{i,j-1}^{n+1}}{h^2}={\left({\mathrm{\Γ}}_{22}^1\right)}_{i,j}^n\frac{f_{i+1,j}^n-f_{i-1,j}^n}{2h}+{\left({\mathrm{\Γ}}_{22}^2\right)}_{i,j}^n\frac{f_{i,j+1}^n-f_{i,j-1}^n}{2h}+\frac{N_{i,j}^n}{\sqrt{E_{i,j}^n+G_{i,j}^n-1}}\end{array}\right.---\left(1\right)$

其中，表示迭代过程中网格点(x_i，y_j)上第n次迭代的模拟值，h表示网格分辨率，

(i＝11，12，22；j＝1，2)表示地形曲面的第二类克里斯托费尔符号。

设

$F^{n+1}=(f_{\mathrm{1,1}}^{n+1},f_{\mathrm{1,2}}^{n+1},...,f_{1,J-1}^{n+1},f_{1,J}^{n+1},f_{\mathrm{2,1}}^{n+1},f_{\mathrm{2,2}}^{n+1},...,f_{2,J-1}^{n+1},f_{2,J}^{n+1},.........,f_{I-\mathrm{1,1}}^{n+1},f_{I-\mathrm{1,2}}^{n+1},...,$

$f_{I-1,J-1}^{n+1},f_{I-1,J}^{n+1},f_{I,1}^{n+1},f_{I,2}^{n+1},...,f_{I,J-1}^{n+1},f_{I,J}^{n+1}{)}^T,$ n≥0

则方程组(1)可以用矩阵形式表达为：

$\left\{\begin{array}{c}{\mathrm{AF}}^{n+1}=b^n\\ {\mathrm{BF}}^{n+1}=c^n\end{array}\right.---\left(2\right)$

其中，A和bⁿ分别为HASM主方程组中第1个方程的系数矩阵和右端项向量；B和cⁿ分别为HASM主方程组中第2个方程的系数矩阵和右端项向量。设：

$P=\begin{array}{c}\left[\begin{array}{c}A\\ B\end{array}\right]\end{array},q^n=\left[\begin{array}{c}{b}^n\\ c^n\end{array}\right],$

则HASM方法可以表示为如下等式约束最小二乘问题：

$\left\{\begin{array}{c}\min {\left|\right|{\mathrm{PF}}^{n+1}-q^n\left|\right|}_2\\ s.t.S\×F^{n+1}=t\end{array}\right.$

其中S∈R^T×(I×J)和t∈R^T×1分别为采样矩阵和采样向量，T为采样点个数。建立等式约束最小二乘方法的目的是为了在保证采样点处模拟值接近采样值的条件下，曲面的整体模拟误差最小。

由于有限差分方法的计算特点，对模拟区域的边界，迭代过程中其高程不变，均采用初始插值，即

(n＞0，i＝0或j＝0)。

以上是离散点曲面建模的方法，对等高线反演DEM，为了充分利用等高线信息，保证回放等高线的质量，我们采用区域控制，构建不等式约束的方法。

二、不等式约束的HASM算法

在上述得到的等式约束的HASM算法的基础上，本发明提出了一种不等式约束的HASM算法。其中，不等式约束条件是通过下述方法获得的：

具体的，如图1所示，为根据等高线确定离散网格点上下界约束示意图。其中，等高线间距为2.5m，在每个等高线附近均标注出该等高线的高程值，比如：高程值为117.5、122.5等。将等高线所在的模拟区域离散化为网格点形式，即离散化为X轴方向有10个网格点，Y轴方向有7个网格点的网格点形式，每一个网格点的中心点处被标注有该网格点的标志号，即：20-89。

计算每一个网格点的高程值范围包括以下步骤：

(一)判断位于所述网格点中心位置处的格网中心点与最接近等高线间距离是否小于预设阈值；如果判断结果为是，则所述网格点高程值为所述最接近等高线的高程值加减预设微小值所得到的值。本步骤中，为了方便求解，通过加减预设微小值将网格点高程值转换为不等式约束问题。

例如：对于图1中标志号为43的网格点，可以看出，等高距为117.5的等高线正好穿过该网格点的格网中心点，因此，标志号为43的网格点的高程值即为117.5。类似的，标志号为22、27、43、62、56和57等网格点的高程值均为穿过其网格点中心点的等高线的高程值。

另外，可以将预设阈值定义为参数ε，当格网中心点与最接近等高线间距离小于ε时，同样可以将该格网中心点所在的网格点定义为采样点，其高程值仍为穿过其网格点中心点的等高线的高程值。其中，ε值的大小根据实际精度需要的程度确定。

(二)如果判断结果为否，则包括以下情况：

(1)当所述格网中心点被两条父子关系的等高线包围时，所述格网中心点所在的网格点的高程值上下界分别为所述两条父子关系中父等高线的高程值和子等高线的高程值。

例如：对于标志号为51的网格点，位于等高值为117.5和等高值为120这两条父子等高线之间，所以网格点51的高程值，应该满足下述不等式l(51)＜f(51)＜u(51)，其中l(51)＝117.5和u(51)＝120。

(2)当所述格网中心点为马鞍处的点时，所述格网中心点所在的网格点的高程值上下界分别为包围该格网中心点的所有等高线高程值的极大值和极小值。

(3)当所述格网中心点为山峰点时，所述格网中心点所在的网格点的高程值下界为包围该格网中心点的最接近等高线的高程值；并进一步判断所述格网中心点所在的山峰是否还存在采样点，如果判断结果为是，则所述格网中心点所在的网格点的高程值上界为该采样点的高程值；如果判断结果为否，则所述格网中心点所在的网格点的高程值上界为包围该格网中心点的最接近等高线的高程值加上1倍的等高距的值。

例如：对于标志号为46的网格点，其所在区域为山峰，，其中还存在一个离散采样点A，其高程值为134，因此，可以近似认为，采样点A是该山峰的高程最大值点。因此，对网格点46，其高程值应该满足不等式l(46)＜f(46)＜u(46)，其中l(46)＝132.5，u(46)＝134。

若网格点46所在的山峰没有采样点，则其高程值所满足的不等式中，高程值上界u(46)可以取为.135，即为围成这片区域等高线的等高值加上1倍的等高距。

(4)当所述格网中心点为山谷点时，所述格网中心点所在的网格点的高程值上界为包围该格网中心点的最接近等高线的高程值；并进一步判断所述格网中心点所在的山谷是否还存在采样点，如果判断结果为是，则所述格网中心点所在的网格点的高程值下界为该采样点的高程值；如果判断结果为否，则所述格网中心点所在的网格点的高程值下界为包围该格网中心点的最接近等高线的高程值减去1倍的等高距的值。

通过上述方法，得到了网格点上下界约束方程组，即：在HASM迭代过程中，曲面应该一直满足下述不等式组：

l＜Fⁿ⁺¹＜u

其中，l＝(l_1，1，l_1，2，…，…，l_1，J-1，l_1，J，l_2，1，l_2，2，…，…，l_2，J-1，l_2，J，…，…l_I-1，1，l_I-1，2，…，…，

l_I-1，J-1，l_I-1，J，l_I，1，l_I，2，…，…，l_I，J-1，l_I，J)^T

u＝(u_1，1，u_1，2，…，…，u_1，J-1，u_1，J，u_2，1，u_2，2，…，…，u_2，J-1，u_2，J，…，…u_I-1，1，u_I-1，2，…，…，

u_I-1，J-1，u_I-1，J，u_I，1，u_I，2，…，…，u_I，J-1，u_I，J)^T

这里l_i，j、u_i，j分别表示网格点(x_i，y_j)(0＜i＜I+1，0＜j＜J+1)高程的两个控制值。当等高线本身很密，等高线闭合的面积很小而结果DEM分辨率很粗时，等高线不包含任何网格的中心点，HASMOC方法(即：本发明所提出的基于曲面论和优化控制理论的曲面建模方法，以下简称HASMOC)模拟结果的回放等高线就会遗漏该条等高线，如高程值为132.5这条等高线，虽然跨越了网格点53、54、63、64这4个网格点，但是都没有包含任何一个格网的中心点，故在模拟过程中就会被忽略，高程值为135这条等高线也将被忽略，所以，欲使该等高线不被忽略，一种有效的方法是增加模拟分辨率，将格网尺寸减小。

另外，通过构建等高线树的方式同样可以依照上述方法获得各采样点的高程值的上下界，本发明在此不再赘述。

综上所述，针对等高线反演DEM的HASMOC方法得到如下等式和不等式约束优化方程组，即下述第一不等式约束方程组：

$\left\{\begin{array}{c}\min {\left|\right|{\mathrm{PF}}^{n+1}-q^n\left|\right|}_2\\ s.t.\\ s\&times;F^{n+1}=t\\ l<F^{n+1}<u\end{array}\right.$

如图2所示，为本发明提供的一种基于曲面论和优化控制理论的曲面建模方法的流程示意图：包括：

(1)采样等高线高程值和已知离散点高程值，并设置迭代参数初始值，所述迭代参数包括：迭代次数K和容差delta；

(2)将等高线所在的原始模拟区域空间离散化为网格点形式，然后根据所述等高线高程值和已知离散点高程值计算所述原始模拟区域中未采样网格点的高程值上下界，得到所述原始模拟区域中所有网格点的高程值上下界；

(3)用所述原始模拟区域中所有网格点的高程值上下界作为约束方程建立第一不等式约束方程组；

(4)根据得到的网格点的高程值计算所述原始模拟区域中每个网格点的第一类基本量E、F、G和第二类基本量L、M、N；其中，所述第一类基本量用于表示地形曲面上曲线的弧长、地形曲面上两个方向的夹角和地形曲面域的面积；所述第二类基本量用于刻画地形曲面空间中的弯曲性；并结合曲面论中的高斯方程以及所述第一不等式约束方程组建立地形曲面的差分微分方程：S*f^k+1＝t^k；其中，S是曲面方程的系数矩阵，f是待求的数字高程模型数据，上标k表示第k次迭代，f^k+1表示第k次迭代得到的地形曲面，t为计算出的地形曲面常数项；

(5)利用迭代算法解算所述S*f^k+1＝t^k，每次迭代结束均得出数字高程模型的模拟结果；其中，迭代运算的初始值S(0)(k＝0)为预先设定值；

(6)在迭代过程中，需要判断相临两次曲面数据的差值S(k)-S(k-1)是否小于所述容差，如果判断结果为是，则重复执行步骤(4)(5)中的迭代过程；如果判断结果为否，则终止迭代运算，执行步骤(7)；

(7)输出地形曲面数据S(k)，得到高精度的数字地形模型。

由于微分方程有限差分迭代过程中，边界不变，

(n＞0，i＝0或j＝0)，这在边界误差比较大的情况下会引起计算的极大困难，因为初始边界，有可能是不能满足等高线的区域控制的。但实际上，对边界，根据其网格点所在的区域，与内部区域一样处理，可以获得其高程的双边不等式。为了对模拟区域边界进行控制优化，我们采取模拟区域扩展的方法。

将原始模拟区域Ω分别向外围扩展M个栅格，获得扩展区域

所述扩展区域包含内部区域和新增区域两部分；其中，所述内部区域范围与所述原始模拟区域范围相同；

扩展区域

表达为下述形式：

$\left\{({\stackrel{\&OverBar;}{x}}_i,{\stackrel{\&OverBar;}{y}}_j)\right|0\&le;i\&le;I+1+2M,0\&le;j\&le;J+1+2M\},$

显然，原始模拟区域对应的网格与扩展区域对应的网格，存在如下对应关系：

$(x_i,y_j)=({\stackrel{\&OverBar;}{x}}_{i+M},{\stackrel{\&OverBar;}{y}}_{j+M}),(0\&le;i\&le;I+\mathrm{1,0}\&le;j\&le;J+1),$

在扩展区域

上，类似地离散HASM主方程，获得相应的代数方程组，设

表示迭代过程中网格点

上第n次迭代的模拟值，又：

${\stackrel{\&OverBar;}{F}}^{n+1}=({\stackrel{\&OverBar;}{f}}_{\mathrm{1,1}}^{n+1},{\stackrel{\&OverBar;}{f}}_{\mathrm{1,2}}^{n+1},...,...,{\stackrel{\&OverBar;}{f}}_{1,J+2M}^{n+1},{\stackrel{\&OverBar;}{f}}_{\mathrm{2,1}}^{n+1},{\stackrel{\&OverBar;}{f}}_{\mathrm{2,2}}^{n+1},...,...,{\stackrel{\&OverBar;}{f}}_{2,J+2M}^{n+1},...,...,{\stackrel{\&OverBar;}{f}}_{I+2M-\mathrm{1,1}}^{n+1},{\stackrel{\&OverBar;}{f}}_{I+2M-\mathrm{1,2}}^{n+1},...,$

${\stackrel{\&OverBar;}{f}}_{I+2M-1,J+2M}^{n+1},{\stackrel{\&OverBar;}{f}}_{I+2M,1}^{n+1},{\stackrel{\&OverBar;}{f}}_{I+2M,2}^{n+1},...,...,{\stackrel{\&OverBar;}{f}}_{I+2M,J+2M}^{n+1}{)}^T,$ n≥0

在扩展区域

上，类似地建立约束优化问题，得到下述第二不等式约束方程组：

$\left\{\begin{array}{c}\min {\left|\right|\stackrel{\&OverBar;}{P}{\stackrel{\&OverBar;}{F}}^{n+1}-{\stackrel{\&OverBar;}{q}}^n\left|\right|}_2\\ s.t.\\ \stackrel{\&OverBar;}{S}\&times;{\stackrel{\&OverBar;}{F}}^{n+1}\text{=}\stackrel{\&OverBar;}{t}\\ \stackrel{\&OverBar;}{l}<{\stackrel{\&OverBar;}{F}}^{n+1}<\stackrel{\&OverBar;}{u}\end{array}\right.$

值得注意的是，对于扩展区域

中的新增区域，各网格点的高程值上下界是通过外推方法获得的，因此精度很可能不高，因此，其内部区域(原始模拟区域Ω)的模拟值在扩展区域

上迭代到一定次数K后，误差可能不再容易衰减，因此，这时候，需要重新回到原始模拟区域Ω，利用已经优化好的Ω的边界，再次进行迭代，

我们把在扩展区域

上的迭代，即：根据第二不等式约束方程组在扩展区域上的迭代称为前迭代；当得到的相临两次曲面数据的差值S(k)-S(k-1)小于预设值时，转入后迭代；所述后迭代为根据第一不等式约束方程组在所述原始模拟区域上的迭代。由于前迭代，其已经满足等高线的区域控制，因此，在此基础上进行计算，会提高计算的精度。如图3所示，为本发明提供的另一种基于曲面论和优化控制理论的曲面建模方法的流程示意图。

使用Matlab 7.0中的lsqlin函数来解决不等式约束的最小二乘的曲线拟合问题：

$\underset{x}{\min }{\left|\right|C\&CenterDot;x-d\left|\right|}_2^2,$ 满足 $\left\{\begin{array}{c}A\&CenterDot;x\&le;b,\\ \mathrm{Aeq}\&CenterDot;x=\mathrm{beq},\\ \mathrm{lb}\&le;x\&le;\mathrm{ub}.\end{array}\right.$

x＝lsqlin(C，d，A，b，Aeq，beq，lb，ub，x0，options)，这里的C是一个m×n的稀疏矩阵，它由HASM的第1个与第3个方程离散化后的系数组成，1b和ub表示的从等高线中求得的变量的上下界，x0是待求量的初始值，options是关于优化求解参数的选项，本文的优化求解的参数皆为默认参数。约束条件方程A·x≤b与Aeq·x＝beq可以用来表示地形特征线及湖泊等约束条件。由于本文只研究了最简单的情形，没有A·x≤b与Aeq·x＝beq这些约束条件，故都为空。

通过使用本发明所提供的基于曲面论和优化控制理论的曲面建模方法，可以获得高精度的DEM。

如图4所示，为通过使用本发明所提供的HASMOC方法得到的DEM，如图5所示，为使用传统的TIN方法得到的DEM。从图中可以明显的看出，图4中的DEM精度远高于图5所示的DEM的精度。

另外，江西省千烟洲生态试验站位于江西省中部，隶属于吉安市泰和县，处于115°04’13”E，26°44’48”N，是一个由约80个小山丘，9条沟谷组成的小山村，总面积为204.16hm2，海拔高程60-150m，属典型的红壤丘陵区。地貌类型主要有4种：①低丘，分布面积最广，地形起伏，丘顶浑圆，海拔多在100m以下，最高达147m，相对高度20-50m，丘坡坡度以10-30°居多。②阶地，地表平坦，坡度为1-2°，相对高度2m左右。③高河漫滩，分布较零散，地表坡度为1-3°，相对高度0.6-1m。④低河漫滩，地表略有起伏，相对高度0.5-0.6m。如图6所示，为千烟洲扫描矢量化原始等高线图和离散点数据。如图7所示，为通过使用本发明所提供的HASMOC方法对图6所示的等高线反演DEM得到的阴影图。如图8所示，为对图7所得到的DEM阴影图进行等高线回放所得到的等高线回放图。对比图6和图8，可以看出，HASMOC模拟结果的千烟洲回放等高线与原始等高线基本吻合，从而有效的证明了本发明所提供的HASMOC方法建立DEM的高保真性。

本发明所提供的基于曲面论和优化控制理论的曲面建模方法应用广泛，可应用于空间图像数据处理中，例如：双星定位导航系统和国防信息化建设中。

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

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

综上所述，通过使用本发明所提供的曲面建模方法，可以建立高精度的DEM，从而能够有效提高双星定位系统的精度以及更好的服务于国防信息化建设中。

Claims (10)

1.一种基于曲面论和优化控制理论的曲面建模方法，其特征在于，包括以下步骤：

(1)采集等高线高程值和已知离散点高程值，并设置迭代参数初始值，所述迭代参数包括：迭代次数k和容差delta；

(2)将等高线所在的原始模拟区域空间离散化为网格点形式，然后根据所述等高线高程值和已知离散点高程值计算所述原始模拟区域中未采样网格点的高程值上下界，得到所述原始模拟区域中所有网格点的高程值上下界；

(3)用所述原始模拟区域中所有网格点的高程值上下界作为约束方程建立第一不等式约束方程组；

(4)根据得到的网格点的高程值计算所述原始模拟区域中每个网格点的第一类基本量E、F、G和第二类基本量L、M、N；其中，所述第一类基本量用于表示地形曲面上曲线的弧长、地形曲面上两个方向的夹角和地形曲面域的面积；所述第二类基本量用于刻画地形曲面空间中的弯曲性；并结合曲面论中的高斯方程以及所述第一不等式约束方程组建立地形曲面的差分微分方程：S*f^k+1＝t^k；其中，S是曲面方程的系数矩阵，f是待求的数字高程模型数据，上标k表示第k次迭代，f^k+1表示第k次迭代得到的地形曲面，t为计算出的地形曲面常数项；

(5)利用迭代算法解算所述S*f^k+1＝t^k，每次迭代结束均得出数字高程模型的模拟结果；其中，迭代运算的初始值S(0)(k＝0)为预先设定值；

(6)在迭代过程中，需要判断相临两次曲面数据的差值S(k)-S(k-1)(k＞＝1)是否小于所述容差，如果判断结果为否，则重复执行步骤(4)(5)中的迭代过程；如果判断结果为是，则终止迭代运算，执行步骤(7)；

(7)输出地形曲面数据S(k)，得到高精度的数字地形模型。

2.根据权利要求1所述的方法，其特征在于，所述根据所述等高线高程值和已知离散点高程值计算所述原始模拟区域中未采样网格点的高程值上下界，得到所述原始模拟区域中所有网格点的高程值上下界包括：

直接在网格点形式的基础上并根据所述等高线高程值和已知离散点高程值计算所述原始模拟区域中未采样网格点的高程值上下界，得到所述原始模拟区域中所有网格点的高程值上下界。

3.根据权利要求2所述的方法，其特征在于，所述计算所述原始模拟区域中未采样网格点的高程值上下界包括以下步骤：

判断位于所述网格点中心位置处的格网中心点与最接近等高线间距离是否小于预设阈值；

如果判断结果为是，则所述网格点高程值上下界分别为所述最接近等高线的高程值加减预设微小值所得到的值；

如果判断结果为否，则包括以下情况：

当所述格网中心点被两条父子关系的等高线包围时，所述格网中心点所在的网格点的高程值上下界分别为所述两条父子关系中父等高线的高程值和子等高线的高程值；

当所述格网中心点为马鞍处的点时，所述格网中心点所在的网格点的高程值上下界分别为包围该格网中心点的所有等高线高程值的极大值和极小值；

当所述格网中心点为山峰点时，所述格网中心点所在的网格点的高程值下界为包围该格网中心点的最接近等高线的高程值；并进一步判断所述格网中心点所在的山峰是否还存在采样点，如果判断结果为是，则所述格网中心点所在的网格点的高程值上界为该采样点的高程值；如果判断结果为否，则所述格网中心点所在的网格点的高程值上界为包围该格网中心点的最接近等高线的高程值加上1倍的等高距的值；

当所述格网中心点为山谷点时，所述格网中心点所在的网格点的高程值上界为包围该格网中心点的最接近等高线的高程值；并进一步判断所述格网中心点所在的山谷是否还存在采样点，如果判断结果为是，则所述格网中心点所在的网格点的高程值下界为该采样点的高程值；如果判断结果为否，则所述格网中心点所在的网格点的高程值下界为包围该格网中心点的最接近等高线的高程值减去1倍的等高距的值。

4.根据权利要求1所述的方法，其特征在于，所述根据所述等高线高程值和已知离散点高程值计算所述原始模拟区域中未采样网格点的高程值上下界，得到所述原始模拟区域中所有网格点的高程值上下界包括以下步骤：

建立所述原始模拟区域中各等高线的等高线树；

根据建立的等高线树以及所述等高线高程值和已知离散点高程值计算所述原始模拟区域中未采样网格点的高程值上下界。

5.根据权利要求1所述的方法，其特征在于，在步骤(1)之后还包括：向外围扩展所述原始模拟区域，获得扩展区域，所述扩展区域包含内部区域和新增区域两部分；其中，所述内部区域范围与所述原始模拟区域范围相同；

利用外推的方法获得所述新增区域中各网格点的高程值上下界；

根据所述原始模拟区域和所述扩展区域的扩展关系，获得所述内部区域中各网格点的高程值上下界；

根据获得的新增区域中各网格点的高程值上下界和内部区域中各网格点的高程值上下界，在所述扩展区域上建立第二不等式约束方程组。

6.根据权利要求5所述的方法，其特征在于，所述步骤(5)中对所述地形曲面的差分微分方程进行迭代时，还包括下述步骤：

首先进行前迭代，所述前迭代为根据所述第二不等式约束方程组在所述扩展区域上的迭代；当得到的相临两次曲面数据的差值S(k)-S(k-1)(k＞＝1)小于预设值时，转入后迭代；所述后迭代为根据所述第一不等式约束方程组在所述原始模拟区域上的迭代。

7.根据权利要求1所述的方法，其特征在于，所述地形曲面为蒙格面片的形式：z＝f(x，y)，其中，x，y分别表示x坐标和y坐标。

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

F＝f_xf_y

其中，所述f_x、f_y分别表示函数f对地形曲面点的横、纵坐标求一阶偏导数。

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

$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}}$

其中，f_x、f_y分别表示函数f对地形曲面点的横、纵坐标求一阶偏导数；f_xx、f_xy、f_yy分别表示函数f对地形曲面点的横、纵坐标求二阶偏导数。

10.一种空间图像数据处理方法，其特征在于，应用权利要求1至9任一项所述的曲面建模方法。

Images