一种生成栅格插值的方法和装置 【技术领域】
本申请涉及地理信息领域， 特别是涉及一种生成栅格插值的方法和装置。背景技术 GIS(Geographic Information System， 地理信息系统 ) 经过了 40 年的发展， 到今 天已经逐渐成为一门相当成熟的技术， 并且得到了极广泛的应用。尤其是近些年， GIS 更以 其强大的地理信息空间分析功能， 在 GPS 及路径优化中发挥着越来越重要的作用。 GIS 地理 信息系统是以地理空间数据库为基础， 在计算机软硬件的支持下， 运用系统工程和信息科 学的理论， 科学管理和综合分析具有空间内涵的地理数据， 以提供管理、 决策等所需信息的 技术系统。在 GIS 应用中， 需要对栅格数据进行插值计算， 即需要在较大地理空间范围下快 速进行栅格插值计算。
现有技术中， 地理空间范围较大条件下进行插值计算， 目前有以下做法 ：
方 法 一、 采 用 大 样 本 量 计 算。 假 设 样 本 点 为 P1， P2， ...， Pn. 对 于 每 一 待 计 算的栅格 V， 常用的算法都是根据样本点与栅格 V 的距离来确定一组系数 a1， a2， ...， an， 这系数是与样本点对应的， 表示每个样本点在插值结果中的贡献度， 即 Value(V) ＝ a1*P1+a2*P2+...+an*Pn ； 此种方法精度高， 但当样本量 n 非常大、 待插值的栅格非常多的 时候， 计算速度非常慢。
方法二、 采用小样本量计算。计算原理同上。由于样本量小， 每个栅格的计算明显 加快， 但是， 对于地理空间范围非常大的条件下， 使用较少的采样点， 由于插值的结果没有 实际采样调查准确可靠， 插值的结果的精确度往往非常低。
方法三、 改进的大样本量计算法。采集的样本点还是 P1， P2， ...， Pn， (n 比较大， 譬如 n ＞ 1000)， 但是， 在对每个栅格计算插值时， 采用大样本量中较少的点参与实际计算， 譬如只选用 C1， C2， ...， Cm(m 值较小， 譬如 m ＝ 10)。这样， 可保证精度的同时提高速度， 但 是， 由于只针对每个栅格选用小部分样本点参与计算， 就会出现相邻的 2 个栅格点 A 和 B， 选 择的参与计算的样本不同， 这样 A 和 B 插值的结果就有可能出现差异较大， 在生成的栅格图 存在明显的阶梯斑快。
发明内容
本申请所要解决的技术问题是提供一种生成栅格插值的方法和装置， 消除生成栅 格图的存在明显的阶梯斑块问题。
为了解决上述问题， 本申请公开了一种生成栅格插值的方法， 其特征在于， 包括 ：
从整个栅格阵列一端的 n×n 子栅格阵列开始， 依次对所有 n×n 子栅格阵列进行 处理， 处理完毕后， 将整个栅格阵列进行输出 ； 其中， 相邻近的两个栅格阵列相互交叠至少 一排栅格 ；
其中， 对于每个 n×n 子栅格阵列进行如下处理 ：
所交叠的各栅格的高程值以相互交叠的已进行平滑处理的子栅格阵列的结果为准； 针对该 n×n 子栅格阵列中的未获得高程值的每个栅格， 基于对应各栅格所选择的多个 采样点的各高程值进行插值计算获得每个栅格的高程值 ；
根据所在子栅格阵列的各栅格的高程值， 按顺序对每个栅格进行平滑处理。
优选的， 所述的平滑处理包括如下步骤 ：
步骤 111， 从所述子栅格阵列其中一个栅格开始， 对于该栅格， 取所在子栅格阵列 的所有栅格的高程值的均值作为该栅格当前的高程值 ；
步骤 112， 针对处理顺序中的下一个栅格， 取所在子栅格阵列的所有栅格当前的高 程值的均值作为该栅格的高程值。
优选的， 根据如下步骤选择对应所述栅格的多个采样点 ：
针对子栅格阵列中的一子栅格 Ti， 计算该栅格的中心点到所有采样点的距离 Li， 取述各距离 Li 中最小的 m 个距离及其对应的采样点 Sj。
优选的， 所述的通过插值计算获得每一子栅格的高程值方法包括 ：
根据计算所述栅格的高程值， 其中 Vi 表示该栅格 Ti 的高程值。优选的， 所述的 n 至少为 2。
相应的， 本发明还公开了一种快速生成栅格插值的装置， 其特征在于， 包括 ：
主栅格处理模块和子栅格处理模块，
所述的主栅格处理模块用于从整个栅格阵列一端的 n×n 子栅格阵列开始， 依次 对所有 n×n 子栅格阵列进行处理， 处理完毕后， 将整个栅格阵列进行输出 ； 其中， 相邻近的 两个栅格阵列相互交叠至少一排栅格 ；
所述的子栅格处理模块包括 ：
高程值获取模块， 用于所交叠的各栅格的高程值以相互交叠的已进行平滑处理的 栅格阵列的结果为准 ； 针对该 n×n 子栅格阵列中的未获得高程值的每个栅格， 基于对应各 栅格所选择的多个采样点的各高程值进行插值计算获得每个栅格的高程值 ；
平滑处理模块， 用于根据所在子栅格阵列的各栅格的高程值， 按顺序对每个栅格 进行平滑处理。
优选的， 所述的平滑处理模块包括 ：
子模块一， 用于从所述子栅格阵列其中一个栅格开始， 对于该栅格， 取所在子栅格 阵列的所有栅格的高程值的均值作为该栅格当前的高程值 ； 子模块二， 用于针对处理顺序中的下一个栅格， 取所在子栅格阵列的所有栅格当 前的高程值的均值作为该栅格的高程值。
优选的， 所述的高程值获取模块包括 ：
采样点获取模块， 用于针对子栅格阵列中的一子栅格， 计算该栅格的中心点到所 有采样点的距离 Li， 取述各距离 Li 中最小的 m 个距离及其对应的采样点 Sj。
与现有技术相比， 本申请包括以下优点 ：
本申请通过以整个栅格阵列中 n×n 的子栅格阵列为单元根据由采样点直接获得 的高程值进行平滑处理， 其中相邻近的后一子栅格阵列以相互交叠的已进行平滑处理的栅 格阵列的结果为准， 循环进行平滑处理， 这样使在生成栅格图时， 速度快， 效率高， 并且去除 了阶梯斑块问题， 得到平滑的栅格图。
附图说明 图 1 是本申请一种本申请一种生成栅格插值的方法的对子栅格阵列的处理流程 示意图 ；
图 2 是本申请一种本申请一种快速生成栅格插值的装置的结构示意图 ；
图 3 是采用方法三大样本量计算得到的栅格图 ；
图 4 是采用方法一大样本量计算得到的栅格图 ；
图 5 是采用本申请大样本量计算得到的栅格图。
具体实施方式
为使本申请的上述目的、 特征和优点能够更加明显易懂， 下面结合附图和具体实 施方式对本申请作进一步详细的说明。
本申请通过以整个栅格阵列中 n×n 的子栅格阵列为单元根据由采样点直接获得 的高程值进行平滑处理， 其中相邻近的后一子栅格阵列以相互交叠的已进行平滑处理的栅 格阵列的结果为准， 循环进行平滑处理， 这样使在生成栅格图时， 速度快， 效率高， 并且去除 了阶梯斑块问题， 得到平滑的栅格图。 本申请的一种生成栅格插值的方法包括 ：
从整个栅格阵列一端的 n×n 子栅格阵列开始， 依次对所有 n×n 子栅格阵列进行 处理， 处理完毕后， 将整个栅格阵列进行输出 ； 其中， 相邻近的两个栅格阵列相互交叠至少 一排栅格。
一个待处理的 A*B(A， B 值比较大 ) 栅格阵列包含大量的栅格， 比如 10000×10000 的待处理栅格阵列。一般可以从左下角的一个 n×n 子栅格阵列开始， 从左向右以 n×n 子 栅格阵列为单元进行处理， 其中相邻近的两个栅格阵列相互交叠至少一排栅格， 比如对上 述 A*B 栅格阵列， 从左下角的第一个 2×2 子栅格阵列开始进行处理， 当对该 2×2 子栅格阵 列处理完毕后， 向右选取第二个 2×2 子栅格阵列， 该子栅格阵列的左边一排与前面的子栅 格阵列的右边一排交叠， 即第一个子栅格阵列与第二个子栅格阵列公用一排栅格， 如此对 第二个子栅格阵列进行处理， 如此循环 ；
当最下面一排作为第一排的所有 2×2 子栅格阵列处理完成后， 再往上进行第二 排的 2×2 子栅格阵列的处理， 当从下向上对子栅格阵列进行处理时， 比如说从左下角往上 的第二排第一个 2×2 子栅格阵列， 则该 2×2 子栅格阵列的下边一排与左下角的 2×2 子栅 格阵列的上边一排交叠， 如此对第二排的第一个 2×2 子栅格阵列进行处理 ； 从左至右， 对 于第二排的第二个 2×2 子栅格阵列， 其左边的一排与第二排的第一个 2×2 子栅格阵列的 右边一排交叠， 其下边的一排与第一排的第二个 2×2 子栅格阵列的上边交叠， 即整个第二 排 2×2 子栅格阵列的下边排栅格排的 2×2 子栅格阵列上边一排栅格交叠， 如此循环， 直至 将所有 2×2 子栅格阵列处理完毕。其中， 子栅格的处理顺序可以有其它方式， 比如， 从右上 角开始， 从左至右， 从上之下对子栅格阵列进行处理等方式。
所述的 n 至少为 2， 即子栅格阵列至少为 2×2 阵列， 比如也可以为 3×3、 4×4 或者 其他。一般情况下， 所述的 n×n 的大小选择为整个栅格阵列大小的二千分之一左右， 同时 还可以根据实际精度的需求确定 n×n 的大小。
另外， 相邻近的两个栅格阵列相互交叠的排数优选为一排栅格， 如果 n×n 相对较 大， 也可根据实际情况选择交叠多排栅格阵列。
在以 n×n 子栅格阵列对整个栅格阵列进行处理时， 当处理到边界而原栅格阵列 缺少多列或者多行形成 n×n 子栅格阵列， 那么可将边界的一行或者一列复制多行或者多 列补齐， 在此基础上完成 n×n 子栅格阵列处理后， 将复制的多行或多列的栅格删除。比如 当从左至右， 从下至上依次以 3×3 子栅格阵列进行处理时， 相邻近两个 3×3 子栅格阵列相 互交叠一排栅格， 如果从左至右处理到右边界时， 原栅格阵列只剩余一列栅格， 而要使最右 边形成一个 3×3 子栅格阵列， 原栅格阵列需要剩余两列栅格， 此时可将原栅格阵列最右边 一列复制添补齐到该列栅格右边， 这样即可继续以 3×3 子栅格阵列进行处理， 处理完毕后 将所述补齐的一列栅格删除即可 ； 如果从下至上， 也缺少一行， 可通过以上相同原理进行处 理。其它情况可类似进行处理。
其中， 对于每个 n×n 子栅格阵列进行如下处理 ：
参照图 1， 示出了本申请一种本申请一种生成栅格插值的方法的对子栅格阵列的 处理流程示意图， 包括 ：
步骤 101， 选定与已进行平滑处理的 n×n 子栅格阵列相邻近的相互交叠至少一排 栅格的 n×n 子栅格阵列。 通过前述次序和规则选择了 n×n 子栅格阵列转入步骤 102。
步骤 102， 所交叠的各栅格的高程值以相互交叠的已进行平滑处理的子栅格阵列 的结果为准 ； 针对该 n×n 子栅格阵列中的未获得高程值的每个栅格， 基于对应各栅格所选 择的多个采样点的各高程值进行插值计算获得每个栅格的高程值。
其中， 对于一个 n×n 子栅格阵列中的未获得高程值的栅格， 可以通过如下方式计 算所述栅格的高程值 ： 比如栅格阵列大小为 A*B， 如果样本量为 S1、 S2、 ...， Sn， 样本量比较 大， 确定本方法计算使用的小样本量 m 取值 ( 例如 m ＝ 10)。对于子栅格阵列中每个要计算 高程值的栅格 Ti(1 ＜＝ i ＜＝ A*B)， 计算 Ti 到每个样本 Si 的距离 Li， 取 Li 中最小的 m 个， 并保存这最小的 n 个距离所对应的样本 Sj(1 ＜ j ＜＝ m)。那么， 待插值点网格 Ti 的值就 是由 Sj(1 ＜ j ＜＝ m) 确定。可根据每个样本点 Sj 到 Ti 的距离大小， 按线性进行插值计 算， 计算得到的结果为 Vi， 这样就得到了该子栅格阵列的各栅格的高程值。其中， 还可以通 过其他方法获取采样点。
其 中， 所述的通过插值计算获得每一子栅格的高程值方法包括： 根据
计算所述栅格的高程值， 其中 Vi 表示该栅格 i 的高程值。 还可以通过其 他的插值方法获得高程值， 比如趋势面法、 样条函数法、 如克立格 (Kriging) 插值法、 谢别 德法 (Shepard′ s Method) 等。
步骤 103， 根据所在子栅格阵列的各栅格的高程值， 按顺序对每个栅格进行平滑处 理。
对于已获得高程值的 n×n 子栅格阵列， 对其进行平滑处理， 优选的， 平滑处理步 骤包括 ：
步骤 111， 从所述子栅格阵列其中一个栅格开始， 对于该栅格， 取所在子栅格阵列 的所有栅格的高程值的均值作为该栅格当前的高程值 ；步骤 112， 针对处理顺序中的下一个栅格， 取所在子栅格阵列的所有栅格当前的高 程值的均值作为该栅格的高程值。 ‘
比如对于 4 个栅格 Pa， Pb， Pc， Pd 通过前述步骤获得的插值结果 Va， Vb， Vc， Vd， 然 后对这 4 个栅格分别进行平滑处理， 处理过程可采用均值法 ： 先对栅格 Pa 进行处理， 计算 Va， Vb， Vc， Vd 的均值 Vh， Pa 的高程值取 Vh ； 然后对栅格 Pb 进行处理， 计算 Vh， Vb， Vc， Vd 的均值 Vt， Pb 的高程值取 Vt ； 再对栅格 Pc 进行处理， 计算 Vh， Vt， Vc， Vd 的均值 Vs， Pc 的 高程值取 Vs ； 再对栅格 Pd 进行处理， 计算 Vh， Vt， Vs， Vd 的均值 Vg， Pd 的高程值取 Vg。这 样， 得到比较平滑的 4 栅格插值结果。
然后， 转入步骤 101， 进入下一个子栅格阵列的处理。
下面以前述对于依次以 2×2 子栅格阵列对整个 A*B 栅格进行处理加以说明对子 栅格阵列的循环处理过程 ：
从左下端开始， 从左至右， 从下至上的顺序依次选择 2×2 子栅格阵列 ； 步骤 101 先 选择左下角第一个 ( 也即最下面一排第一个 )2×2 子栅阵列， 其中 4 个栅格均未获得高程 值， 则通过上述步骤 102 计算 4 个栅格的高程值， 然后转入步骤 103， 对该 4 个栅格进行平滑 处理， 得到各栅格处理后的高程值 ； 然后进入步骤 101， 选择最下面一排第二个 2×2 子栅格阵列， 该子栅格阵列的左 边一排与第一个子栅格阵列右边一排交叠， 所以该子栅格阵列的左边一排的栅格已经获得 高程值， 其余的栅格未获得高程值， 则转入步骤 102， 使所述未获得高程值的栅格获得高程 值， 然后转入步骤 103， 对该 2×2 子栅格阵列进行同样的平滑处理 ；
然后转入下一个 2×2 子栅格阵列的处理， 其过程如上所述 ； 如此循环所有 2×2 子 栅格阵列处理完毕后， 整个栅格阵列的处理即完毕， 则将整个栅格阵列进行输出。
为了便于说明本申请的效果， 以实际大小为 4000*4000 象素的栅格进行地面高度 插值计算， 采样点为 130 个点的高程值 ：
其中， 空间范围为一个长方形， 左下角坐标 ： (x， y) ＝ (420000， 252010) ； 右上角坐 标： (x， y) ＝ (600000， 436010)。
若采用方法一， 大样本量计算。130 个采样点全部参与计算， 约需要 1 天得时间完 成得到的结果如图 4 所示， 其精确度比较高， 并且平滑， 但时间非常慢。
若采用方法二， 小样本量计算。显然精度不够， 实际中不采用。
若采用方法三， 改进的大样本量计算法。使用 10 个采样点参与每象素的计算， 通 过距离选择采样点， 计算时间约为 8 分钟， 得到结果为图 3 所示， 虽然时间快， 但存在明显的 阶梯斑块。
若采样本申请的方法， 使用所有样本中 10 个采用点参与每象素的计算， 通过距离 选择采样点， 计算时间约为 9 分钟， 得到的结果为图 5 所示， 该结果与方法一十得到的结果 十分近似， 说明精度高， 并且时间快， 而且生产的栅格图平滑， 在实际应用中效率非常高。
相应的， 本发明还公开了一种生成栅格插值的装置， 参照图 2 所示， 包括 ：
主栅格处理模块 210 和子栅格处理模块 220，
所述的主栅格处理模块 210 用于从整个栅格阵列一端的 n×n 子栅格阵列开始， 依 次对所有 n×n 子栅格阵列进行处理， 处理完毕后， 将整个栅格阵列进行输出 ； 其中， 相邻近 的两个栅格阵列相互交叠至少一排栅格 ；
所述的子栅格处理模块 220 包括 ：
高程值获取模块， 用于所交叠的各栅格的高程值以相互交叠的已进行平滑处理的 栅格阵列的结果为准 ； 针对该 n×n 子栅格阵列中的未获得高程值的每个栅格， 基于对应各 栅格所选择的多个采样点的各高程值进行插值计算获得每个栅格的高程值 ；
平滑处理模块， 用于根据所在子栅格阵列的各栅格的高程值， 按顺序对每个栅格 进行平滑处理。
优选的， 所述的平滑处理模块包括 ：
子模块一， 用于从所述子栅格阵列其中一个栅格开始， 对于该栅格， 取所在子栅格
阵列的所有栅格的高程值的均值作为该栅格当前的高程值 ；
子模块二， 用于针对处理顺序中的下一个栅格， 取所在子栅格阵列的所有栅格当 前的高程值的均值作为该栅格的高程值。
优选的， 所述的高程值获取模块包括 ：
采样点获取模块， 用于针对子栅格阵列中的一子栅格， 计算该栅格的中心点到所 有采样点的距离 Li， 取述各距离 Li 中最小的 m 个距离及其对应的采样点 Sj。其中所述的采 样点获取模块还可以通过其他方法获取采样点。
对于系统实施例而言， 由于其与方法实施例基本相似， 所以描述的比较简单， 相关 之处参见方法实施例的部分说明即可。
本说明书中的各个实施例均采用递进的方式描述， 每个实施例重点说明的都是与 其他实施例的不同之处， 各个实施例之间相同相似的部分互相参见即可。
以上对本申请所提供的一种生成栅格插值的方法和装置， 进行了详细介绍， 本文 中应用了具体个例对本申请的原理及实施方式进行了阐述， 以上实施例的说明只是用于帮 助理解本申请的方法及其核心思想 ； 同时， 对于本领域的一般技术人员， 依据本申请的思 想， 在具体实施方式及应用范围上均会有改变之处， 综上所述， 本说明书内容不应理解为对 本申请的限制。