地理坐标下网格化数据等值线处理具体步骤及编程方法

刘鹂

（爱丁堡大学，英国爱丁堡 EH8 9YL）

在研究地理坐标时需要对大量的数据进行处理，用计算机进行网格化数据处理，需要运用编程的方法对数据进行分析。网格点的数据处理需要一些专业的功能，执行功能的专业软件有很多，如冲浪软件、GOLDEN软件等。但是，在运用专业的软件进行编程时需要处理大量的数据节点，并且对其进行进一步处理，目前的专业软件在执行这些功能价格昂贵，非常不方便，使用不到位还会引起侵权等纠纷。本文分析了网格数据轮廓处理的原理和方法，并给出了用Microsoft C++5.0编写的源程序的相应功能。

1.等值线型数据网格化方法的研究

地学和地球化学勘探中的网格划分通常包括3类：研究数据网格划分、轮廓数据网格划分和离散点数据网格划分[1]。不同类型的数据必须采用不同的网格划分方法才能取得良好的效果。对于研究数据，三次样条插值网格化方法有效；对于离散点类型数据，可以采用平方反比加权平均grid ding方法或二次样条插值grid ding方法[2]。当前，等高线数据的网格化实用方法很少。根据轮廓数据的特点，模拟了轮廓图中手动检索网格的方法和过程。提出了一种分块存储等高线数据的网格方法，并在8个方向上查找和插值网格点。主要的想法和过程如下。

1.1 网格图幅条块划分

根据网格地图范围的δx，δy网格线距离，计算出N×M网格线的网格数；然后，根据正确选择的r插值搜索半径，将网格地图表拆分到的较小数据存储区域中的网格虚线数计算为LN×LM：

LN=R/ΔX

IF (LN＊XK.LT.RS) LN=R/ΔX+1

LM=R/ΔY

IF (LM＊YK.LT.RS) LM=R/ΔY+1

则数据存储小方块的宽度为XLN×YLM:

XLN=LN×ΔX

YLM=LM×ΔY

MW=M/LM

IF (MW＊LM.LT.M) MW=M/LM+1

NW=N/LN

IF (NW＊LN.LT.N) NW=N/LN+1

注：＊代表指针变量

/代表注释标志

1.2 数据内插加密及坐标转换

为了确保以下栅格点的八方位插值范围内的插值计算点和研究点的准确性，首先，必须在大约1mm间隔（可以选择特定间隔）的数据点处沿每个轮廓插值并加密飞机轮廓的数字数据并且必须将网格地图范围内的数据点坐标转换并旋转到用户定义的网格地图坐标系，以便于后续步骤中的计算工作。

1.3 等值线型数据分块存储

缩小要搜索的原始数据点的范围，可以加快网格点的搜索和插值速度。I型轮廓数据点（x，y，t）应存储在序列号较小的正方形中（IX+1，JY+1）：

IX=X/XLN+1 (1 ≤ IX+1 ≤ NW+2)

JY=Y/YLM+1 (1≤JY+1≤MW+2)

JY3=MOD (JY, 3) +1 (1 ≤ JY3 ≤ 3)

将数据点 (X, Y, T) 存放在以下几个内存数组中：

MNG (IX+1, JY3) =MNG (IX+1, JY3) +1

XN (MNG (IX+1, JY3) , IX+1, JY3) =X

YM (MNG (IX+1, JY3) , IX+1, JY3) =Y

TMN (MNG (IX+1, JY3) , IX+1, JY3) =T

注： /代表注释标志

MNG是一个内存数组，用于计算存储在小正方形中的数据点数。收集的大量的数据点一般情况下无法全部同时间录入计算机内，大量的数据点在计算机内无法进行分析。从上面列出的公式中可以得到，JY3的搜索范围是1～3，因此沿栅格y方向有MW+2条带块，并且只有3个条带（每个条带都是NW+的）的轮廓数据点（仅计算整个栅格表的轮廓数据的大约3/[MW+2]），输出到外部存储设备以按块顺序访问数据文件，并将边界数据点块存储在栅格图中的所有MW+2（NW+2）×（MW+2）小块中，并成为按块顺序存储的新边界数据文件。

1.4 网格点八方位搜索插值

人工轮廓图中网格数据的获取方法是一种简单的线性插值方法。对网格正方形中的每个网格数据点执行插值计算时，只需在内存中读取的3个相邻条带正方形之间的网格广场附近的9个正方形中查找轮廓数据点，即可计算整个映射中的9/[NW+2]×（MW+2）轮廓数据点因此，isoline数据点搜索次数非常少，数据搜索时间也节省了。传统搜索方法使用的时间仅为9/[ （NW+2）×（MW+2）]左右，对于包含大量数据的等腰数据来说，这是相当可观的。根据插值加密轮廓的扫描数据点密度和栅格图轮廓密度，选择合适的插值条带宽度，将搜索半径r用作插值条带长度，并以栅格点为中心形成八向插值范围。从相邻的9个小正方形中的轮廓数据点开始，查找位于八向补间范围内且距离网格点最近的轮廓数据点之一。数据点搜索结果分为以下3种情况：

（1）如果在4个相对于8个方向的插值条集中，一个插值条集中没有2个等轴测数据点，表明栅格点位于栅格图纸上的空白区域或等轴测边缘，则可以为栅格点分配假值。

（2）如果在8个方向上的4个相对补间动画范围集中，至少有2个具有不相等值的等轴测数据点位于一个补间范围集中，则可以选择最接近的补间范围集中具有不相等值的2个数据点，并且可以获取网格点中的值。

（3）如果在8个方向上有4组相反的补间范围，则每个补间范围中2个轮廓数据点的值相等，这表示该网格点位于网格图的轮廓端点值或边缘的半闭合圆的闭合圆上。在这种情况下，从相邻的9个小正方形中的轮廓数据点开始，我们必须再次搜索位于8个方向上的补间范围内的另一个轮廓数据点，该点距离网格点最近，但其值不等于上一个点，因此我们可以选择最接近方向的点。

采用上述方法，完成NW+1小网目条纹块网格点8位搜索插值计算，条纹块网格点计算结果存储在外部存储设备的块随机访问数据文件中完成网格图中所有MW+1方形网格点的八向搜索插值计算，得出网格图原始区域的网格数据结果。

1.5 大范围网格图幅分幅计算

为解决大数据大比例尺等高线图的网格计算问题，可以根据正确选择的地图尺寸，沿x方向将大比例尺栅格图分割成多个相对较小的栅格图，然后重复上述步骤1～4，并在计算每个小网格地图时获得的网格数据结果依次存储在外部存储设备上的块随机访问数据文件中，最终转化为原始大规模网格地图区域的统一网格数据结果[3]。

2.网格化数据等值线处理

2.1 均匀网格化数据等值线处理

均匀网格化数据等值线处理方法是计算等值线f(x，y)=c和分布均匀的矩形网格域的每个矩形网格边缘的交点坐标，然后将交点进行连接。通常我们先找到岛的起点，然后沿着岛一点一点地找到交点，同时可以得到连接顺序。具体步骤如下。

（1）先找到起始点，搜索域的外边缘，查看这些边缘是否有起始点。如果没有起点，可以根据从下到上的顺序进行搜索，然后再根据从左到右的顺序进行搜索；

（2）查找交点以确定孤岛是否与网格边缘相交，并线性查找交点；

（3）边界线在网格中移动；

（4）曲线平滑拟合终点节点的确定、三次样条曲线拟合和抛物线拟合。

2.2 非均匀网格分布数据的均匀网格化处理

实际上，在收集的过程中采样点分布不均匀，导致数据不均匀，结果是不能满足绘图要求，对此，需要对收集到的数据进行网格划分和数据点加密，因此采用了曲面处理方法。

2.2.1 按距离加权的最小二乘法（N—P法）

N-P方法只能计算距离网格点最近的n个数据点。提供大量{xi，yi，zi}ni=1i=1参数时，假设要查找点（a，b）的高度值，可以找到多项式。通常采用二次多项式：

f(x，y)=c1+c2x+c3y+c4xy+c5x2+c6y2。

使用最小二乘法尽可能地拟合数据点。与一般最小二乘法的含义不同，邻近数据点（a，b）的权重必须大于远处数据点（a，b），也就是说，必须选择ci系数以最小化距离加权最小二乘法的作用。

2.2.2 按方位取点加权法

上述N-P方法是取离栅格点最近的n数据值来确定栅格点的值，而不管方向如何，所取的点很可能集中在一侧或两侧，而不会在其他方向上取点。

本节中介绍的方法将面积划分为多个以栅格点为中心的象限，并将每个象限中的一个点作为加权平均值，从而克服了N-P方法的缺点。基本原则如下。

要查找给定栅格点（I，j）的函数值，请将平面拆分为4个基本象限（以（I，j）为原点），然后将每个象限拆分为n0部分。这样，整架飞机分成4n0等份。然后，查找每个拆分角度上最近的数据点（I，j），其值为Zil，其到（I，j）的距离为ril。网格（I，j）中的值为：

Z(i‚j)=∑i1=14n0Ci1 ·Zi1Z(i‚j)=∑i1=14n0Ci1 ·Zi1,

式中参数

Ci1=Πj=14n0j≠i1r2j∑K=14n0ΠL=14n0L≠Kr2LCi 1=Πj=14n0j≠i1rj2∑K=14n0ΠL=14n0L ≠KrL2。

2.2.3 加权最小二乘法拟合法（M—S法）

M-S方法不仅可以全面分析所有的数据点，对数据整体趋势进行分析，而且具有N-P方法反映局部特征的优点。

假设必须计算网格点（a，b）处曲面f（a，b）的高度，则需要得到p（x，y）多项式（通常为二次多项式）。

p（x，y）=c00+C10x+c01y+C20x2+c11xy+C02y2

总的趋势分析应该是

q=I=1n[p（xi-yi）-zi]2q = I = 1n[φ（xi-yi）-zi]2

是最小值，其中：n是数据点的数目，（xi，yi）是数据点的坐标，zi是（xi，yi）中观察到的值，而装配值p（xi，yi）和观察值zi之间的误差平方和是q。

M-S方法应考虑距离加权，并对上述公式做一些修改。在这里，仅仅以p（x，y）二次多项式为例。

Q=∑i=1n[P（xi‚yi）-zi]2·W[（xi-a）2+（yi-b）2]Q=∑i=1n[p（xi‚yi）-zi]2 ·W[（xi-a）2+（yi-b）2]

式中：

W[(xi-a)2+(yi-b)2]W[(xi-a)2+(yi-b)2]

就是权，它是距离的函数。

为了求方程中的Crs系数，按照最小二乘法的原则，应该是

∂Q∂Crs=0,（r,s=0,1,2）∂Q∂Crs=0,（r,s=0，1,2）

所提供的程序一共给出3种加权形式：

Ⅰ型

W（d2）=1（d2+ε）W（d2）=1（d2+ε）

式中：d为网点，（a,b）到离散数据点（xi,yi）的距离，即d=(xi-a)2+(yi-b)2。

Ⅱ型

W（d2）=1（d2+ε）2，W(d2)=1(d2+ε)2，d 与 ε的意义同Ⅰ型。

Ⅲ型

W(d2)=exp(-ad2)(d2+ε)W(d2)=exp（-ad2）（d2+ε）

3.编程实现

本文主要运用了gridding数据的isolina处理原理和方法，下面是具体编程。

//N—P法

float nds (int n, float a, float b, float＊x, float＊y,float＊z, int m)

{int i, j, ic;float sl, s2, d, dis[5000];

//S—M法

float wlsa (float a,float b, int np,float ＊ x,float＊y,float ＊z, int k)

{int i, j;float x1, y1, x2, y2, term, xt, yt, xxt, yyt,xyt, zt, e[6][7], u[6], zzz;

注：＊代表指针变量；/代表注释标志。

4.结语

通过对上述方法的描述和计算实例，可以总结出适用于8个网格点位置的异构数据块存储和搜索插值grid ding方法的以下结论。

（1）由于采用了大规模计算技术，将等轴测数据存储在条带中，并具备条带网格点计算能力，因此可以解决大规模计算、海量数据和大型网络等问题，使这种方法非常实用。

（2）由于采用了八格线点插值方法，保证了格线点的高插值精度，很好地解决了地图边缘轮廓线的判断和插值，从而确保了的效果。

（3）运用两种编程方法对网络化数据等值线进行处理。