地理坐標下網格化數據等值線處理具體步驟及編程方法

劉鸝

（愛丁堡大學，英國愛丁堡 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方法的缺點?；驹瓌t如下。

要查找給定柵格點（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）運用兩種編程方法對網絡化數據等值線進行處理。