地理坐標下網(wǎng)格化數(shù)據(jù)等值線處理具體步驟及編程方法

劉鸝

（愛丁堡大學，英國愛丁堡 EH8 9YL）

在研究地理坐標時需要對大量的數(shù)據(jù)進行處理，用計算機進行網(wǎng)格化數(shù)據(jù)處理，需要運用編程的方法對數(shù)據(jù)進行分析。網(wǎng)格點的數(shù)據(jù)處理需要一些專業(yè)的功能，執(zhí)行功能的專業(yè)軟件有很多，如沖浪軟件、GOLDEN軟件等。但是，在運用專業(yè)的軟件進行編程時需要處理大量的數(shù)據(jù)節(jié)點，并且對其進行進一步處理，目前的專業(yè)軟件在執(zhí)行這些功能價格昂貴，非常不方便，使用不到位還會引起侵權(quán)等糾紛。本文分析了網(wǎng)格數(shù)據(jù)輪廓處理的原理和方法，并給出了用Microsoft C++5.0編寫的源程序的相應功能。

1.等值線型數(shù)據(jù)網(wǎng)格化方法的研究

地學和地球化學勘探中的網(wǎng)格劃分通常包括3類：研究數(shù)據(jù)網(wǎng)格劃分、輪廓數(shù)據(jù)網(wǎng)格劃分和離散點數(shù)據(jù)網(wǎng)格劃分[1]。不同類型的數(shù)據(jù)必須采用不同的網(wǎng)格劃分方法才能取得良好的效果。對于研究數(shù)據(jù)，三次樣條插值網(wǎng)格化方法有效；對于離散點類型數(shù)據(jù)，可以采用平方反比加權(quán)平均grid ding方法或二次樣條插值grid ding方法[2]。當前，等高線數(shù)據(jù)的網(wǎng)格化實用方法很少。根據(jù)輪廓數(shù)據(jù)的特點，模擬了輪廓圖中手動檢索網(wǎng)格的方法和過程。提出了一種分塊存儲等高線數(shù)據(jù)的網(wǎng)格方法，并在8個方向上查找和插值網(wǎng)格點。主要的想法和過程如下。

1.1 網(wǎng)格圖幅條塊劃分

根據(jù)網(wǎng)格地圖范圍的δx，δy網(wǎng)格線距離，計算出N×M網(wǎng)格線的網(wǎng)格數(shù)；然后，根據(jù)正確選擇的r插值搜索半徑，將網(wǎng)格地圖表拆分到的較小數(shù)據(jù)存儲區(qū)域中的網(wǎng)格虛線數(shù)計算為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

則數(shù)據(jù)存儲小方塊的寬度為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 數(shù)據(jù)內(nèi)插加密及坐標轉(zhuǎn)換

為了確保以下柵格點的八方位插值范圍內(nèi)的插值計算點和研究點的準確性，首先，必須在大約1mm間隔（可以選擇特定間隔）的數(shù)據(jù)點處沿每個輪廓插值并加密飛機輪廓的數(shù)字數(shù)據(jù)并且必須將網(wǎng)格地圖范圍內(nèi)的數(shù)據(jù)點坐標轉(zhuǎn)換并旋轉(zhuǎn)到用戶定義的網(wǎng)格地圖坐標系，以便于后續(xù)步驟中的計算工作。

1.3 等值線型數(shù)據(jù)分塊存儲

縮小要搜索的原始數(shù)據(jù)點的范圍，可以加快網(wǎng)格點的搜索和插值速度。I型輪廓數(shù)據(jù)點（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)

將數(shù)據(jù)點 (X, Y, T) 存放在以下幾個內(nèi)存數(shù)組中：

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是一個內(nèi)存數(shù)組，用于計算存儲在小正方形中的數(shù)據(jù)點數(shù)。收集的大量的數(shù)據(jù)點一般情況下無法全部同時間錄入計算機內(nèi)，大量的數(shù)據(jù)點在計算機內(nèi)無法進行分析。從上面列出的公式中可以得到，JY3的搜索范圍是1～3，因此沿柵格y方向有MW+2條帶塊，并且只有3個條帶（每個條帶都是NW+的）的輪廓數(shù)據(jù)點（僅計算整個柵格表的輪廓數(shù)據(jù)的大約3/[MW+2]），輸出到外部存儲設備以按塊順序訪問數(shù)據(jù)文件，并將邊界數(shù)據(jù)點塊存儲在柵格圖中的所有MW+2（NW+2）×（MW+2）小塊中，并成為按塊順序存儲的新邊界數(shù)據(jù)文件。

1.4 網(wǎng)格點八方位搜索插值

人工輪廓圖中網(wǎng)格數(shù)據(jù)的獲取方法是一種簡單的線性插值方法。對網(wǎng)格正方形中的每個網(wǎng)格數(shù)據(jù)點執(zhí)行插值計算時，只需在內(nèi)存中讀取的3個相鄰條帶正方形之間的網(wǎng)格廣場附近的9個正方形中查找輪廓數(shù)據(jù)點，即可計算整個映射中的9/[NW+2]×（MW+2）輪廓數(shù)據(jù)點因此，isoline數(shù)據(jù)點搜索次數(shù)非常少，數(shù)據(jù)搜索時間也節(jié)省了。傳統(tǒng)搜索方法使用的時間僅為9/[ （NW+2）×（MW+2）]左右，對于包含大量數(shù)據(jù)的等腰數(shù)據(jù)來說，這是相當可觀的。根據(jù)插值加密輪廓的掃描數(shù)據(jù)點密度和柵格圖輪廓密度，選擇合適的插值條帶寬度，將搜索半徑r用作插值條帶長度，并以柵格點為中心形成八向插值范圍。從相鄰的9個小正方形中的輪廓數(shù)據(jù)點開始，查找位于八向補間范圍內(nèi)且距離網(wǎng)格點最近的輪廓數(shù)據(jù)點之一。數(shù)據(jù)點搜索結(jié)果分為以下3種情況：

（1）如果在4個相對于8個方向的插值條集中，一個插值條集中沒有2個等軸測數(shù)據(jù)點，表明柵格點位于柵格圖紙上的空白區(qū)域或等軸測邊緣，則可以為柵格點分配假值。

（2）如果在8個方向上的4個相對補間動畫范圍集中，至少有2個具有不相等值的等軸測數(shù)據(jù)點位于一個補間范圍集中，則可以選擇最接近的補間范圍集中具有不相等值的2個數(shù)據(jù)點，并且可以獲取網(wǎng)格點中的值。

（3）如果在8個方向上有4組相反的補間范圍，則每個補間范圍中2個輪廓數(shù)據(jù)點的值相等，這表示該網(wǎng)格點位于網(wǎng)格圖的輪廓端點值或邊緣的半閉合圓的閉合圓上。在這種情況下，從相鄰的9個小正方形中的輪廓數(shù)據(jù)點開始，我們必須再次搜索位于8個方向上的補間范圍內(nèi)的另一個輪廓數(shù)據(jù)點，該點距離網(wǎng)格點最近，但其值不等于上一個點，因此我們可以選擇最接近方向的點。

采用上述方法，完成NW+1小網(wǎng)目條紋塊網(wǎng)格點8位搜索插值計算，條紋塊網(wǎng)格點計算結(jié)果存儲在外部存儲設備的塊隨機訪問數(shù)據(jù)文件中完成網(wǎng)格圖中所有MW+1方形網(wǎng)格點的八向搜索插值計算，得出網(wǎng)格圖原始區(qū)域的網(wǎng)格數(shù)據(jù)結(jié)果。

1.5 大范圍網(wǎng)格圖幅分幅計算

為解決大數(shù)據(jù)大比例尺等高線圖的網(wǎng)格計算問題，可以根據(jù)正確選擇的地圖尺寸，沿x方向?qū)⒋蟊壤邧鸥駡D分割成多個相對較小的柵格圖，然后重復上述步驟1～4，并在計算每個小網(wǎng)格地圖時獲得的網(wǎng)格數(shù)據(jù)結(jié)果依次存儲在外部存儲設備上的塊隨機訪問數(shù)據(jù)文件中，最終轉(zhuǎn)化為原始大規(guī)模網(wǎng)格地圖區(qū)域的統(tǒng)一網(wǎng)格數(shù)據(jù)結(jié)果[3]。

2.網(wǎng)格化數(shù)據(jù)等值線處理

2.1 均勻網(wǎng)格化數(shù)據(jù)等值線處理

均勻網(wǎng)格化數(shù)據(jù)等值線處理方法是計算等值線f(x，y)=c和分布均勻的矩形網(wǎng)格域的每個矩形網(wǎng)格邊緣的交點坐標，然后將交點進行連接。通常我們先找到島的起點，然后沿著島一點一點地找到交點，同時可以得到連接順序。具體步驟如下。

（1）先找到起始點，搜索域的外邊緣，查看這些邊緣是否有起始點。如果沒有起點，可以根據(jù)從下到上的順序進行搜索，然后再根據(jù)從左到右的順序進行搜索；

（2）查找交點以確定孤島是否與網(wǎng)格邊緣相交，并線性查找交點；

（3）邊界線在網(wǎng)格中移動；

（4）曲線平滑擬合終點節(jié)點的確定、三次樣條曲線擬合和拋物線擬合。

2.2 非均勻網(wǎng)格分布數(shù)據(jù)的均勻網(wǎng)格化處理

實際上，在收集的過程中采樣點分布不均勻，導致數(shù)據(jù)不均勻，結(jié)果是不能滿足繪圖要求，對此，需要對收集到的數(shù)據(jù)進行網(wǎng)格劃分和數(shù)據(jù)點加密，因此采用了曲面處理方法。

2.2.1 按距離加權(quán)的最小二乘法（N—P法）

N-P方法只能計算距離網(wǎng)格點最近的n個數(shù)據(jù)點。提供大量{xi，yi，zi}ni=1i=1參數(shù)時，假設要查找點（a，b）的高度值，可以找到多項式。通常采用二次多項式：

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

使用最小二乘法盡可能地擬合數(shù)據(jù)點。與一般最小二乘法的含義不同，鄰近數(shù)據(jù)點（a，b）的權(quán)重必須大于遠處數(shù)據(jù)點（a，b），也就是說，必須選擇ci系數(shù)以最小化距離加權(quán)最小二乘法的作用。

2.2.2 按方位取點加權(quán)法

上述N-P方法是取離柵格點最近的n數(shù)據(jù)值來確定柵格點的值，而不管方向如何，所取的點很可能集中在一側(cè)或兩側(cè)，而不會在其他方向上取點。

本節(jié)中介紹的方法將面積劃分為多個以柵格點為中心的象限，并將每個象限中的一個點作為加權(quán)平均值，從而克服了N-P方法的缺點?；驹瓌t如下。

要查找給定柵格點（I，j）的函數(shù)值，請將平面拆分為4個基本象限（以（I，j）為原點），然后將每個象限拆分為n0部分。這樣，整架飛機分成4n0等份。然后，查找每個拆分角度上最近的數(shù)據(jù)點（I，j），其值為Zil，其到（I，j）的距離為ril。網(wǎng)格（I，j）中的值為：

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

式中參數(shù)

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

2.2.3 加權(quán)最小二乘法擬合法（M—S法）

M-S方法不僅可以全面分析所有的數(shù)據(jù)點，對數(shù)據(jù)整體趨勢進行分析，而且具有N-P方法反映局部特征的優(yōu)點。

假設必須計算網(wǎng)格點（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是數(shù)據(jù)點的數(shù)目，（xi，yi）是數(shù)據(jù)點的坐標，zi是（xi，yi）中觀察到的值，而裝配值p（xi，yi）和觀察值zi之間的誤差平方和是q。

M-S方法應考慮距離加權(quán)，并對上述公式做一些修改。在這里，僅僅以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]

就是權(quán)，它是距離的函數(shù)。

為了求方程中的Crs系數(shù)，按照最小二乘法的原則，應該是

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

所提供的程序一共給出3種加權(quán)形式：

Ⅰ型

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

式中：d為網(wǎng)點，（a,b）到離散數(shù)據(jù)點（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.編程實現(xiàn)

本文主要運用了gridding數(shù)據(jù)的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.結(jié)語

通過對上述方法的描述和計算實例，可以總結(jié)出適用于8個網(wǎng)格點位置的異構(gòu)數(shù)據(jù)塊存儲和搜索插值grid ding方法的以下結(jié)論。

（1）由于采用了大規(guī)模計算技術(shù)，將等軸測數(shù)據(jù)存儲在條帶中，并具備條帶網(wǎng)格點計算能力，因此可以解決大規(guī)模計算、海量數(shù)據(jù)和大型網(wǎng)絡等問題，使這種方法非常實用。

（2）由于采用了八格線點插值方法，保證了格線點的高插值精度，很好地解決了地圖邊緣輪廓線的判斷和插值，從而確保了的效果。

（3）運用兩種編程方法對網(wǎng)絡化數(shù)據(jù)等值線進行處理。