-
随着对地观测技术的发展，激光探测与测量(light detection and ranging, LiDAR)、合成孔径雷达干涉(synthetic aperture radar interferometry, InSAR)等技术在地形测量领域已得到广泛的应用，获取得到的地表高程数据的精度越来越高，数据量也越来越大[1]。如何根据海量地表高程观测数据快速、准确地生成高精度数字高程模型(digital elevation model，DEM)，为数字地形分析及地理信息空间表达提供数据支持，已成为GIS领域重要课题之一[2]。为了快速、准确地根据已知采样数据重建复杂地形，引入并发展了反距离权重方法、多项式方法、克里金方法、径向基函数(radial basis function，RBF)方法等[3]，由美国测量学家Hardy最早提出的RBF方法[4]，具有较优的插值性能[5]，在地形建模、LiDAR点云建模、数值逼近等领域被广泛使用。随着采样数据数量的增加，RBF插值模型求解速度变慢，同时，插值矩阵过于庞大，导致插值模型求解困难甚至求解失败[6-7]。随着地形采样数据的增加，采用RBF插值模型对大量采样数据进行地形重建同样会受到此限制。
针对大量采样数据进行RBF插值模型求解速度慢甚至不可解的问题，Beatson等人和Wendland进行了一系列的研究，提出了快速多极化方法[8]、残差迭代方法[9]、单元分解方法[10]、区域分解原理(domain decomposing method, DDM)[11]等方法，分别从不同的角度为RBF插值模型的快速、准确求解提供了解决方法。文献[12]采用多线程技术对等值线数据局部逐点并行插值，但需对每个局部模型的边界进行无缝连接，权函数的选取对插值结果的精度有较大的影响，并且输入数据为等值线数据，应用范围有限。文献[13]提出限制性加性施瓦兹(restricted additive Schwarz method，RASM)方法，将大规模线性矩阵方程求解问题转换为局部求解，文献[14]基于该理论，将RBF求解过程进行了并行求解和数值试验，为基于大量数据RBF地形快速插值重建提供了参考，但其采用的截断高斯(Gauss)基函数，在理论上使得插值结果在数值空间存在不完备性。另一方面，紧支撑径向基函数(compact support RBF，CSRBF)[15]使RBF基函数的影响范围控制在有限的区域，提高了插值模型的可解性，在一定程度上能够提高RBF方法的求解效率，但未能完全解决大数据量的CSRBF插值模型求解问题，且存在基函数紧支撑半径选取的问题。结合区域分解原理与RASM，将基于大量数据RBF插值转换为分块可解，选取合适的CSRBF基函数紧支撑半径[15]，根据RASM的并行特性，对插值模型进行并行求解，可提高RBF插值模型的求解效率和插值精度，为海量地表高程数据快速、准确地生成DEM提供理论参考。
本文提出一种基于CSRBF基函数插值模型的自适应并行地形插值方法，根据地形采样数据的特点选取合适的紧支撑半径，结合区域分解和RASM并行原理，构建地表高程数据并行生成DEM的方法，并通过DEM数据进行实验，对方法的求解精度、运算效率进行验证和分析。
-
假定实验区域Ω内有N个已知点p1, p2, …, pN，及其对应的属性值f1, f2, …, fN，则CSRBF插值模型F(p)可表示为由N个已知点处基函数的线性组合：
$$ \mathit{F}\left( \mathit{p} \right){\rm{ = }}\sum\nolimits_{\mathit{i} = 1}^\mathit{N} {{\mathit{c}_\mathit{i}}\mathit{\varphi }\left( {\left\| {{\mathit{p}_\mathit{i}}{\rm{ - }}\mathit{p}} \right\|} \right)} $$ (1) 选取常见的Wendland C2光滑度的紧支撑基函数[15] φ(r)=(1－r)+4(4r+1)。其中，(1－r)+为紧支撑域的控制函数，当待插值点与已知基函数节点的距离rji大于基函数紧支撑半径rcs时，φ(r)=0。紧支撑半径rcs的选取直接影响插值模型的求解效率和插值结果的精度。紧支撑半径的选取，需要根据研究对象的具体情况确定，使得选取的紧支撑半径既能提高插值结果的求解效率又可保证模型求解精度。
-
本文采用一种基于局部子区域采样点之间距离分布的方法，自适应选取各子区域插值模型的最佳紧支撑半径。其计算过程如下：
假定实验区域Ω的局部子区域Ωk内nk个已知点，两点间的距离函数为dist(•, •)：
步骤1 计算实验区域Ω内的任意两点之间的平均最临近距离dmm:
$$ {\mathit{d}_{{\rm{mm}}}}{\rm{ = mean(min(dist(}}{\mathit{p}_\mathit{i}}{\rm{, }}{\mathit{p}_\mathit{j}}{\rm{))), }}{\mathit{p}_\mathit{i}}{\rm{, }}{\mathit{p}_\mathit{j}} \in {\rm{\Omega }} $$ (2) 式中，mean为求均值函数；min为求最小值函数。
步骤2 对每个Ωk∈Ω，计算Ωk内任意点pjk距离点pik的最小距离的平均值dmkk和Ωk内所有点之间的平均距离dmeank；
将点Ωk最佳紧支撑半径rcsk设置为：
$$ \mathit{r}_{\mathit{cs}}^\mathit{k}{\rm{ = (}}\frac{{\mathit{d}_{\mathit{mk}}^k}}{{\mathit{d}_{{\rm{mm}}}^\mathit{k}}}{\rm{)}}\mathit{d}_{{\rm{mean}}}^\mathit{k} $$ 通过该计算过程，求解每个子区域的最佳紧支撑半径。
-
为使全局CSRBF插值模型能够并行地求解，根据区域分解原理，将实验区域Ω划分为K个相互重叠的大小相同的子区域Ωk(k=1, 2, …, K)，如图 1(a)。
图 1 区域分解图示
Figure 1. DDM Diagram
以原始区域中的一个局部子区域Ωk图 1(a)、图 1(b)中虚线圆所示，图 1(c)中Ωξk即为对Ωk进行放大，如图 1(b)中，从区域中心向外依次为黑色区域Ωψk、灰色条带区域$ \tilde \Omega _{^\mathit{\xi }}^\mathit{k}$、浅灰色条带区域$\tilde \Omega _{^\mathit{\zeta }}^\mathit{k}$，并行求解的重叠子区域为黑色区域Ωψk与灰色带区域$ \tilde \Omega _{^\mathit{\xi }}^\mathit{k}$构成的区域Ωξk，其中区域Ωψk并行计算的核心区域(长度dψl，宽度dψw)，区域$\tilde \Omega _{^\mathit{\xi }}^\mathit{k}$为并行运算的辅助区域，用于控制相邻子区域无缝连接，其条带宽度dξk影响CSRBF插值模型的求解精度和求解效率，浅灰色区域$\tilde \Omega _{^\mathit{\zeta }}^\mathit{k} $为子区域Ωξk的CSRBF插值模型所能影响到的区域，其条带宽度dζk等于Ωξk中基函数节点的紧支撑半径的最大值argmax{rcsi}, i∈Ωξk。
浅灰色条带区域$ \tilde \Omega _{^\mathit{\zeta }}^\mathit{k}$、灰色条带区域$ \tilde \Omega _{^\mathit{\xi }}^\mathit{k}$、黑色区域Ωψk共同构成区域Ωζk，为局部子区域Ωξk进行插值模型迭代求解的残差影响区域。区域划分的准则如下。
(1) 区域Ωk内点数nk≤N×0.01%；
(2) 浅灰色重叠区域$ \tilde \Omega _{^\mathit{\zeta }}^\mathit{k}$的条带宽度dζk=max(rcsi), rcsi∈Ωξk；
(3) 灰色重叠区域$ \tilde \Omega _{^\mathit{\xi }}^\mathit{k}$的条带宽度dξk介于一定范围3dmmk≤dξk≤max(dψl, dψw)；
(4) 每个子区域大小相同：$ \tilde \Omega _{^\mathit{\xi }}^\mathit{i}\left\{ {\mathit{d}_{^\mathit{\xi }}^\mathit{i}} \right\} = \tilde \Omega _{^\mathit{\xi }}^\mathit{j}\left\{ {\mathit{d}_{^\mathit{\xi }}^\mathit{j}} \right\}$，Ωψi{dψl dψw}=Ωψj{dψl dψw}，i, j=1, 2, …K。
由于并行运算过程是将全局插值矩阵Φ基于区域分解原理映射到每个局部子区域Ωk，通过不同的处理器单元processornp对每个子区域的局部插值矩阵Φk进行局部求解，将局部区域的解ck映射到全局解c，然后通过线性方程迭代方法，求取全局区域插值模型的最终解cfinal。因此，区域划分完成后，需要建立局部子区域的点pik坐标与全局区域的点pi坐标索引映射关系IΩk{Ii:pik∈Ωk→pi∈Ω，并根据映射关系IΩk建立局部子区域Ωk的插值矩阵Φk，其中，重叠子区域插值矩阵Φk与非重叠子区域插值矩阵Φk共同构成了插值矩阵$ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}= \left[\begin{array}{l} {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^\mathit{k}}\\ {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^{\mathit{\hat k}}} \end{array} \right]$，而映射关系表IΩk可表示为IΩk=[1 0]，因此，Φk与Φ的关系可通过映射关系IΩk表示为：
$$ {\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^\mathit{k}}{\rm{ = }}{\mathit{\boldsymbol{I}}_{{{\rm{\Omega }}_\mathit{k}}}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} $$ (3) $$ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\rm{ = }}\mathit{\boldsymbol{I}}_{{{\rm{\Omega }}_\mathit{k}}}^{\rm{T}}{\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}^\mathit{k}} $$ (4) -
随着实验区域内数据量的增加，线性矩阵Φc=F中矩阵Φ变得过于庞大，直接求解线性矩阵的解c比较困难，可利用矩阵预条件方法，选取非奇异矩阵P将线性矩阵Φc=F, 等价转换为：
$$ {\mathit{\boldsymbol{P}}^{{\rm{ - 1}}}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} c}}{\rm{ = }}{\mathit{\boldsymbol{P}}^{{\rm{ - 1}}}}\mathit{\boldsymbol{F}} $$ (5) 式(5) 的等价变换，降低了原始线性矩阵中矩阵Φ的条件数，使得矩阵Φ变得可解，提高了矩阵求解效率。RASM求解矩阵预条件，利用区域分解将全局矩阵分解为多个局部矩阵单元，这些局部矩阵单元在空间上具有一定的重叠，从而保证了区域分解后局部单元边界处的连续性。
将局部非重叠子区域Ωψk与全局区域Ω的映射关系的矩阵表示为$ {{\mathit{\boldsymbol{\hat I}}}_{{{\rm{\Omega }}_\mathit{k}}}}$，包含重叠区域的局部子区域Ωk(即Ωξk)与全局区域Ω的映射关系的矩阵表示为IΩk，则局部非重叠子区域插值矩阵ΦΩψkk与全局区域插值矩阵Φ的映射关系可表示为：
$$ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{{\rm{\Omega }}_\mathit{\psi }^\mathit{k}}^\mathit{k}{\rm{ = }}{{\mathit{\boldsymbol{\hat I}}}_{{{\rm{\Omega }}_\mathit{k}}}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} $$ (6) 包含重叠区域的局部子区域插值矩阵ΦΩξkk与全局区域插值矩阵Φ的映射关系可表示为：
$$ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_{{\rm{\Omega }}_\mathit{\xi }^\mathit{k}}^\mathit{k}{\rm{ = }}{\mathit{\boldsymbol{I}}_{{{\rm{\Omega }}_\mathit{k}}}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }} $$ (7) 因此，根据式(3)~式(7) 局部非重叠子区域Ωψk的矩阵预条件表示为：
$$ {\mathit{\boldsymbol{P}}_\mathit{k}}{\rm{ = }}\mathit{\boldsymbol{\hat I}}_{{{\rm{\Omega }}_\mathit{k}}}^{\rm{T}}\mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_\mathit{k}^{ - 1}{\mathit{\boldsymbol{I}}_{{{\rm{\Omega }}_\mathit{k}}}} $$ (8) 式中，$ {{\mathit{\boldsymbol{\hat I}}}_{{{\rm{\Omega }}_\mathit{k}}}}$将矩阵映射范围限制在非重叠子区域Ωψk，因此式(8) 的计算过程不需要与相邻子区域进行数据通信，并可对各子区域进行并行计算。结合RASM预条件，采用灵活广义最小残差(flexible generalized minimal residual, FGMRES)线性求解方法[16]，根据设定的残差阈值条件对线性矩阵进行迭代求解，获得原始线性矩阵的解c。基于RASM预条件的FGMRES迭代求解过程如下，本文将可调紧支撑半径引入到插值矩阵的生成算法为：
输入 原始区域Ω点集{p1, p2, …, pN}及其对应的高程属性值{f1, f2, …, fN}，迭代残差阈值resthreshold；
初始化 (1) 将Ω区域分解为K个相互重叠的子区域{Ω1, Ω2, …, ΩK}；(2) 初始化残差向量res1{res1q=f1, res2q=f2, …, resnq=fn}；(3) 初始化线性方程解c1=0；(4) 求解区域各基函数节点最佳紧支撑半径rcsi，i=1, 2, …, N；(5) 基于最佳紧支撑半径rcsi生成插值矩阵Φ；其中，插值矩阵Φ是区域内已知点之间通过CSRBF基函数φ(r)构成的关系矩阵，将第i个点与第j个点的关系表示为φi, j=φ(‖pi－pj‖)，则矩阵Φ可表示为：
$$ \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\rm{ = }}\left[{\begin{array}{*{20}{c}} {{\mathit{\varphi }_{{\rm{1, 1}}}}}&{{\mathit{\varphi }_{{\rm{1, 2}}}}}& \cdots &{{\mathit{\varphi }_{{\rm{1, }}\mathit{N}}}}\\ {{\mathit{\varphi }_{{\rm{2, 1}}}}}&{{\mathit{\varphi }_{{\rm{2, 2}}}}}& \cdots &{{\mathit{\varphi }_{{\rm{2, }}\mathit{N}}}}\\ \vdots & \vdots & \ddots & \vdots \\ {{\mathit{\varphi }_{\mathit{N}{\rm{, 1}}}}}&{{\mathit{\varphi }_{\mathit{N}{\rm{, 2}}}}}& \cdots &{{\mathit{\varphi }_{\mathit{N}{\rm{, }}\mathit{N}}}} \end{array}} \right] $$ (9) 对求解过程进行迭代，根据每一步的解$ {\mathit{c}_{\mathit{q}{\rm{ + 1}}}}{\rm{ = }}{\mathit{c}_\mathit{q}}{\rm{ + }}\sum\nolimits_{\mathit{k}{\rm{ = 1}}}^\mathit{K} {{{\mathit{\boldsymbol{\hat I}}}_{{{\rm{\Omega }}_\mathit{k}}}}} \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}_\mathit{k}^{ - 1}{\mathit{\boldsymbol{I}}_{{{\rm{\Omega }}_\mathit{k}}}}\left( {\mathit{\boldsymbol{F}} - \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\mathit{c}_\mathit{q}}} \right)$计算迭代残差resq+1，直至残差小于阈值，然后，输出线性矩阵求解结果cfinal。
基于单元分解方法对矩阵进行子区域划分，根据CSRBF紧支撑特性，结合RASM方法，使得插值模型线性矩阵转化为局部可解，基于多处理器并行架构，对该计算过程进行并行处理，不仅提高了矩阵求解的效率，而且大大降低了矩阵求解时计算机连续内存的使用要求，提高了方法的可行性。
-
如果基函数计算范围固定，则基于RASM的问题求解的计算复杂度只和实验区域的基函数节点的数量N有关，算法的复杂度为O(N)。本文算法选取具有紧支撑半径的基函数作为RBF插值模型的基函数，具有有限的局部影响范围，最佳紧支撑半径的求解过程与求解生成RBF插值模型距离矩阵的过程一致，不会增加算法的复杂度，因此，本文算法的时间复杂度为O(N)。
-
本文使用C++编程语言，利用可移植可扩展科学计算工具包(portable extensible toolkit for scientific computation，PETSc)编程实现本文方法，以4核64位处理器，3.1 G主频，4 G内存的计算机为实验平台，验证该方法的求解精度和运行效率。
-
采用美国先进的太空散射和反射辐射计项目(Advanced Spaceborne Thermal Emission and Reflection Radiometer, ASTER)网站(http://asterweb.jpl.nasa.gov/gdem.asp)提供的DEM(30 m×30 m)，选取中国西南某地区数字高程模型(DEM)DEM数据从中随机选取范围为100×100 km2的区域(如图 2)进行实验。对本文方法的精度和效率进行实验分析，假定实验所用的DEM数据是精确的。
图 2 实验区域与插值结果
Figure 2. Experiment Area and Interpolation Result
对实验区域随机采样250 000个高程点进行实验，采用Wendland C2光滑度的基函数对实验区域进行插值模型求解，将实验区域内所有DEM网格点作为误差验证点。根据本文方法得到的插值模型对误差验证点进行插值，计算插值结果与采样点处实际高程值之间的均方根误差(root mean square error，RMS error)，计算方式为：
$$ {\varepsilon _{{\rm{RMS}}}}{\rm{ = }}\sqrt {\frac{1}{M}{{\sum\nolimits_{i = 1}^M {\left\| {\left( {{{\tilde f}_i} - {f_i}} \right)} \right\|} }^2}} ,i{\rm{ = 1}},{\rm{2}}, \ldots ,M $$ (10) 式中，M为验证点的数量；$ {{\mathit{\tilde f}}_\mathit{i}}$为第i个验证点的插值结果；fi为第i个验证点的实际高程。为获取较优的求解速度，本文方法选取子区域的最小点数为30~50，子区域大小根据实验区域点的数量进行设定。子区域核心区域Ωψk划分完成后，CSRBF基函数的紧支撑半径、子区域间的重叠率、MPI并行处理单元的数量直接影响本文方法插值精度和计算效率，需要对这些因素分别进行实验和分析。
-
为了对本文紧支撑半径计算方法进行实验和分析，以并行处理单元np=4，子区域间重叠率为0.9，分别采用不同的紧支撑半径和本文方法的紧支撑半径进行实验，计算不同紧支撑半径和本文方法的CSRBF插值模型的求解速度和误差情况，结果如表 1(“-”表示无计算结果)。由表 1可见，采用不同的紧支撑半径，插值模型求解时间和插值误差不同。随着紧支撑半径的增加，子区域中采样点的数量增加，局部子区域的求解效率降低，插值模型整体求解时间增加，插值结果的误差随着紧支撑半径的增加而减小然后增大。采用较小的紧支撑半径，使得插值模型求解效率增加，但较小的紧支撑半径使得已知采样点之间的高程变得不连续，会导致插值结果精度降低，较大的紧支撑半径使得插值模型的求解变慢甚至导致插值模型不可解。选择合适的紧支撑半径，既能得到较优的插值精度，又可提高插值模型的求解效率。本文紧支撑半径求解方法得到的插值结果精度和运算效率相对较优。
表 1 不同紧支撑半径对插值效率和精度的影响
Table 1. The Relationship Between Compact Support Radius and the Computation Efficiency, RMS error
rcsi/m 280 560 840 1 120 1 400 1 680 1 960 本文方法 时间/s 2.363 0 3.492 4 4.573 8 5.211 0 8.032 8 23.886 1 - 4.630 6 εRMS/m 482.640 2 31.357 3 17.115 1 18.953 0 22.984 0 28.755 7 - 17.325 1 -
设定子区域间的重叠率$ \mathit{\alpha }{\rm{ = }}\frac{{\mathit{d}_{\rm{ \mathit{ ξ} }}^\mathit{k}}}{{{\rm{max(}}\mathit{d}_\mathit{\psi }^\mathit{l}{\rm{, d}}_\mathit{\psi }^\mathit{w}{\rm{)}}}}$，dζk可根据子区域内基函数节点紧支撑半径的大小确定，一般情况下dζk=max(rcsi)，rcsi∈Ωξk, rcsi为子区域Ωk中第i个基函数对应的紧支撑半径。区域分解参数包括局部子区域的大小(dψl dψw)，子区域间重叠区域宽度(dζk dξk)以及子区域内的点数nk。为了便于进行并行运算，要保证子区域的大小相同和重叠区域dξk宽度相同。较大的重叠率会导致插值速度的降低，如表 2。
表 2 不同的区域重叠率对插值效率的影响
Table 2. Relationship Between Overlap Ratio and the Computation Efficiency
α= 0.2 0.4 0.6 0.8 1 1.2 时间/s 3.837 9 3.615 5 3.681 0 4.080 5 4.670 0 5.573 2 由表 2可见，当子区域重叠率较小时，子区域间重叠率较小导致子区域间边界重叠区域存在较大的误差，使得插值模型求解时全局迭代次数增加，插值模型求解时间相对较长，随着子区域间重叠率的增加，子区域之间边界误差降低，插值模型求解时间变短，继续增加重叠率，参与局部区域插值模型求解的已知点数量增加，插值求解时间增加，导致并行计算的效率降低。因此，为了获取较高的计算效率，最佳子区域间重叠率设为0.5。
-
根据PETSc并行计算函数的特性，分别选取MPI并行处理单元的数量np为1、2、4、6、8、16、20、24、28进行实验，运行时间及误差情况如表 3。由表 3可见，计算机处理器物理内核数目固定，随着MPI并行计算单元的增加，本文方法的运算速度变快，超过一定的MPI并行计算单元数量后，运算速度保持稳定。并行运算的速度与计算机物理内核的数量有直接关系，但在一定程度上增加MPI并行计算单元的数量，能够进一步加快插值模型的求解速度，因此，本文选取MPI并行计算单元的数量为10。
表 3 不同数量并行计算单元对插值效率的影响
Table 3. Relationship Between the Number of MPI Process and the Computation Efficiency
np 1 2 4 6 8 16 20 24 28 时间/s 26.721 0 11.403 5 7.181 4 5.863 2 5.672 3 5.456 2 5.481 0 5.421 5 5.509 3 -
分别采用一般CSRBF方法、自适应分层CSRBF方法、本文方法对实验区域进行插值模型求解，RBF基函数统一选取Wendland C2光滑度的基函数，对比几种方法的计算效率和插值精度进行实验，以本文实验区域为研究对象，分别选取不同数量的已知点，进行插值模型的求解，其中，求解线性方程的方法均为FGMRES迭代方法。对插值模型的求解时间进行统计，结果如下(“-”表示无计算结果)。
由表 4、表 5可见，随着已知点数量的增加，不能够直接采用CSRBF方法进行插值模型的求解，自适应分层CSRBF方法能够对大量数据实验区域进行插值模型求解，但计算效率较低，而本文方法能够快速对大量已知地形采样数据进行插值模型的求解。另一方面，CSRBF方法与本文方法获得的插值模型精度一致，自适应分层CSRBF方法需要对子区域边界连接权值进行取舍，影响了插值结果的精度。本文方法在计算效率和插值精度方面性能较优。
表 4 不同插值方法运算效率对比/s
Table 4. Computation Efficiency Comparison Between Different Kinds of CSRBF Interpolation Methods/s
N 1 000 10 000 50 000 100 000 200 000 300 000 500 000 一般CSRBF 165.062 7 - - - - - - 自适应分层CSRBF 23.663 6 24.979 5 57.568 5 190.464 4 841.548 7 2 205.531 3 - 本文方法 0.290 9 0.382 1 1.333 5 2.183 7 3.817 4 5.830 6 9.268 7 表 5 不同插值方法计算误差对比/m
Table 5. RMS Error Comparison Between Different Kinds of CSRBF Interpolation Methods/m
N 1 000 10 000 50 000 100 000 200 000 300 000 500 000 一般CSRBF 107.766 0 - - - - - - 自适应分层CSRBF 133.244 0 106.554 0 97.856 1 101.636 7 94.724 0 97.112 4 - 本文方法 108.692 0 107.222 4 25.875 5 21.982 3 18.011 3 24.941 1 23.212 2 -
基于本文方法，选取500 000个采样点数进行CSRBF插值模型的求解，1 000×1 000个规则格网点进行对实验区域DEM进行插值，区域重叠率为0.5，区域内最小点数30~50，并行计算节点np为10。插值结果如图 2(b)。从图中可以看出，采用500 000个采样点作为基函数节点进行插值得到的结果与原始实验区域基本一致。
-
地形数据是地学分析的基础数据，快速、准确地生成地形数据对地学分析具有重要意义。本文根据区域分解原理，结合CSRBF插值模型基函数的紧支撑特性，采用区域分解和RASM解决了CSRBF方法无法直接对大规模数据进行直接求解的问题，根据地形采样数据快速、准确地获取插值模型并进行地形重建，同时根据局部地形数据自适应选取CSRBF基函数的最佳紧支撑半径，对DEM数据进行实验，获得了较优的实验结果，表明本文方法可行，可作为一种大数据量并行地形插值方法。
-
摘要: 快速、准确地对地形进行重建以生成数字高程模型是地理信息表达的重要研究内容，径向基函数（radial basis function，RBF）作为一种插值性能较优的空间插值方法，特别适合于重建复杂的地形模型，但随着已知地形采样点数量的增加，RBF插值模型求解速度变慢，同时插值矩阵过于庞大而导致插值模型求解困难甚至求解失败。针对这个问题，本文基于区域分解和施瓦兹并行原理进行地形插值，以紧支撑径向基函数（compact support RBF，CSRBF）构建基于所有地形采样数据的全局插值矩阵，并自适应求解子区域CSRBF插值节点紧支撑半径，基于限制性加性施瓦兹方法（restricted additive Schwarz method，RASM）采用多核并行架构对各局部子区域的插值矩阵进行求解。以某地区数字高程模型（DEM）数据进行插值实验，结果表明，本文方法能够对大规模地形数据进行准确重建，并且具有较高的求解效率。Abstract: Fast and accurate 3D terrain reconstruction to acquire a high resolution Digital Elevation Model (DEM) is one of the most important research areas in geographic information representation. As a kind of spatial interpolation method that is superior to other accurate interpolation methods, the Radial Basis Function (RBF) is particularly suitable for the reconstruction of complex 3D terrain models. However, the efficiency when calculating this interpolation model becomes lower and lower with an increasing number of sampling points, and the interpolation equation becomes too difficult or even fails as the interpolation matrix become bigger and bigger. To address this issue, a parallel interpolation method based on the principle of domain decomposition and restricted additive schwarz method (referred to as RASM) is proposed. A compact support RBF (CSRBF) based global interpolation matrix was built by taking all of the known sampling points, and the optimal local compact support radius is calculated for each of the local domains. The interpolation procedure operates in parallel through a message passing interface (MPI) based on the RASM. DEM data are used in an interpolation experiment. The results show that the method proposed in this paper could accurately reconstruct the terrain with massive terrain sampling data enabling a high efficiency solution.
-
Key words:
- terrain reconstruction /
- domain decomposing method /
- parallel /
- adaptive /
- RASM /
- CSRBF
-
表 1 不同紧支撑半径对插值效率和精度的影响
Table 1. The Relationship Between Compact Support Radius and the Computation Efficiency, RMS error
rcsi/m 280 560 840 1 120 1 400 1 680 1 960 本文方法 时间/s 2.363 0 3.492 4 4.573 8 5.211 0 8.032 8 23.886 1 - 4.630 6 εRMS/m 482.640 2 31.357 3 17.115 1 18.953 0 22.984 0 28.755 7 - 17.325 1 
下载: 导出CSV
表 2 不同的区域重叠率对插值效率的影响
Table 2. Relationship Between Overlap Ratio and the Computation Efficiency
α= 0.2 0.4 0.6 0.8 1 1.2 时间/s 3.837 9 3.615 5 3.681 0 4.080 5 4.670 0 5.573 2
下载: 导出CSV
表 3 不同数量并行计算单元对插值效率的影响
Table 3. Relationship Between the Number of MPI Process and the Computation Efficiency
np 1 2 4 6 8 16 20 24 28 时间/s 26.721 0 11.403 5 7.181 4 5.863 2 5.672 3 5.456 2 5.481 0 5.421 5 5.509 3
下载: 导出CSV
表 4 不同插值方法运算效率对比/s
Table 4. Computation Efficiency Comparison Between Different Kinds of CSRBF Interpolation Methods/s
N 1 000 10 000 50 000 100 000 200 000 300 000 500 000 一般CSRBF 165.062 7 - - - - - - 自适应分层CSRBF 23.663 6 24.979 5 57.568 5 190.464 4 841.548 7 2 205.531 3 - 本文方法 0.290 9 0.382 1 1.333 5 2.183 7 3.817 4 5.830 6 9.268 7
下载: 导出CSV
表 5 不同插值方法计算误差对比/m
Table 5. RMS Error Comparison Between Different Kinds of CSRBF Interpolation Methods/m
N 1 000 10 000 50 000 100 000 200 000 300 000 500 000 一般CSRBF 107.766 0 - - - - - - 自适应分层CSRBF 133.244 0 106.554 0 97.856 1 101.636 7 94.724 0 97.112 4 - 本文方法 108.692 0 107.222 4 25.875 5 21.982 3 18.011 3 24.941 1 23.212 2
下载: 导出CSV
-
[1] 刘坡, 龚建华.大规模遥感影像全球金字塔并行构建方法[J].武汉大学学报·信息科学版, 2016, 41(1):117-122 http://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201601020.htm Liu Po, Gong Jianhua. Parallel Construction of Global Pyramid for Large Remote Sensing Images[J].Geomatics and Information Science of Wuhan University, 2016, 41(1):117-122 http://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201601020.htm [2] 李志林, 朱庆.数字高程模型[M].武汉:武汉大学出版社, 2001:125-139 Li Zhilin, Zhu Qing. Digital Elevation Model[M]. Wuhan, China:Wuhan University Press, 2001:125-139 [3] 汤国安, 杨昕. ArcGIS空间分析实验教程[M].北京:科学出版社, 2006, 363-422 Tang Guoan, Yang Xin. ArcGIS Spatial Analysis Experiment Tutorial[M]. Beijing, China:Science Press, 2006:363-422 [4] Hardy R L.Theory and Applications of the Multiquadric-biharmonic Method 20 Years of Discovery 1968-1988[J]. Computers & Mathematics with Applications, 1990, 19(8):163-208 http://www.sciencedirect.com/science/article/pii/089812219090272L [5] Franke R. Scattered Data Interpolation:Tests of Some Methods[J]. Mathematics of Computation, 1982, 38(157):181-200 http://www.ams.org/mathscinet-getitem?mr=637296 [6] 段平, 盛业华, 张思阳, 等.顾及异向性的局部径向基函数三维空间插值[J].武汉大学学报·信息科学版, 2015, 40(5):632-637 http://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201505013.htm Duan Ping, Sheng Yehua, Zhang Siyang, et al. A 3D Local RBF Spatial Interpolation Considering Anisotropy[J]. Geomatics and Information Science of Wuhan University, 2015, 40(5):632-637 http://www.cnki.com.cn/Article/CJFDTOTAL-WHCH201505013.htm [7] 吕海洋, 盛业华, 段平, 等.局部最优形态参数的RBF分块地形插值方法与实验[J].地球信息科学学报, 2015, 17(3):260-267 http://www.cnki.com.cn/Article/CJFDTOTAL-DQXX201503003.htm Lv Haiyang, Sheng Yehua, Duan Ping, et al.A Hierarchical RBF Interpolation Method Based on Local Optimal Shape Parameters[J]. Journal of Geo-information Science, 2015, 17(3):260-267 http://www.cnki.com.cn/Article/CJFDTOTAL-DQXX201503003.htm [8] Carr J C, Beatson R K, Cherrie J B, et al. Reconstruction and Representation of 3D Objects with Radial Basis Functions[C]. Proceedings of the 28th Annual Conference on Computer Graphics and Interactive Techniques. New York, USA, 2001 [9] Beatson R K, Cherrie J B, Mouat C T. FastFitting of Radial Basis Functions:Methods Based on Preconditioned GMRES Iteration[J]. Advances in Computational Mathematics, 1999, 11(2-3):253-270 http://www.ams.org/mathscinet-getitem?mr=1731700 [10] WendlandH. Fast Evaluation of Radial Basis Functions:Methods Based on Partition of Unity[M]//Approximation Theory X:Wavelets, Splines, and Applications. Nashville, Tennessee, USA:Vanderbilt University Press, 2002 [11] Beatson RK, Light W A, Billings S. Fast Solution of the Radial Basis Function Interpolation Equations:Domain Decomposition Methods[J]. SIAM Journal on Scientific Computing, 2001, 22(5):1717-1740 doi: 10.1137/S1064827599361771 [12] Wang Q, Pan Z, Bu J, et al. Parallel RBF-based Reconstruction from Contour Dataset[C]. Computer Aided Design and Computer Graphics, 200710th IEEE International Conference, Piscataway, USA, 2007 [13] Cai X C, Sarkis M. A Restricted Additive Schwarz Preconditioner for General Sparse Linear Systems[J]. SIAM Journal on Scientific Computing, 1999, 21(2):792-797 doi: 10.1137/S106482759732678X [14] Yokota R, Barba L A, Knepley M G. PetRBF-A Parallel O(N) Algorithm for Radial Basis Function Interpolation with Gaussians[J]. Computer Methods in Applied Mechanics and Engineering, 2010, 199(25):1793-1804 http://www.sciencedirect.com/science/article/pii/S0045782510000666 [15] Fasshauer G E. Meshfree Approximation Methods with MATLAB[M]. Singapore:World Scientific, 2007 [16] Saad Y. A Flexible Inner-outer Preconditioned GMRES Algorithm[J]. SIAM Journal on Scientific Computing, 1993, 14(2):461-469 doi: 10.1137/0914028 -
点击查看大图
图(2) / 表(5)
计量
- 文章访问数: 1909
- HTML全文浏览量: 178
- PDF下载量: 316
- 被引次数: 0
