劉鸝

（愛丁堡大學(xué)，英國(guó)愛丁堡 EH8 9YL）

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

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

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

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

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

注：＊代表指針變量

/代表注釋標(biāo)志

1.2 數(shù)據(jù)內(nèi)插加密及坐標(biāo)轉(zhuǎn)換

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

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

縮小要搜索的原始數(shù)據(jù)點(diǎn)的范圍，可以加快網(wǎng)格點(diǎn)的搜索和插值速度。I型輪廓數(shù)據(jù)點(diǎn)（x，y，t）應(yīng)存儲(chǔ)在序列號(hào)較小的正方形中（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ù)點(diǎn) (X, Y, T) 存放在以下幾個(gè)內(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

注： /代表注釋標(biāo)志

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

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

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

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

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

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

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

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

為解決大數(shù)據(jù)大比例尺等高線圖的網(wǎng)格計(jì)算問題，可以根據(jù)正確選擇的地圖尺寸，沿x方向?qū)⒋蟊壤邧鸥駡D分割成多個(gè)相對(duì)較小的柵格圖，然后重復(fù)上述步驟1～4，并在計(jì)算每個(gè)小網(wǎng)格地圖時(shí)獲得的網(wǎng)格數(shù)據(jù)結(jié)果依次存儲(chǔ)在外部存儲(chǔ)設(shè)備上的塊隨機(jī)訪問數(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ù)等值線處理方法是計(jì)算等值線f(x，y)=c和分布均勻的矩形網(wǎng)格域的每個(gè)矩形網(wǎng)格邊緣的交點(diǎn)坐標(biāo)，然后將交點(diǎn)進(jìn)行連接。通常我們先找到島的起點(diǎn)，然后沿著島一點(diǎn)一點(diǎn)地找到交點(diǎn)，同時(shí)可以得到連接順序。具體步驟如下。

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

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

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

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

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

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

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

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

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

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

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

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

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

要查找給定柵格點(diǎn)（I，j）的函數(shù)值，請(qǐng)將平面拆分為4個(gè)基本象限（以（I，j）為原點(diǎn)），然后將每個(gè)象限拆分為n0部分。這樣，整架飛機(jī)分成4n0等份。然后，查找每個(gè)拆分角度上最近的數(shù)據(jù)點(diǎn)（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ù)點(diǎn)，對(duì)數(shù)據(jù)整體趨勢(shì)進(jìn)行分析，而且具有N-P方法反映局部特征的優(yōu)點(diǎn)。

假設(shè)必須計(jì)算網(wǎng)格點(diǎn)（a，b）處曲面f（a，b）的高度，則需要得到p（x，y）多項(xiàng)式（通常為二次多項(xiàng)式）。

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

總的趨勢(shì)分析應(yīng)該是

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

是最小值，其中：n是數(shù)據(jù)點(diǎn)的數(shù)目，（xi，yi）是數(shù)據(jù)點(diǎn)的坐標(biāo)，zi是（xi，yi）中觀察到的值，而裝配值p（xi，yi）和觀察值z(mì)i之間的誤差平方和是q。

M-S方法應(yīng)考慮距離加權(quán)，并對(duì)上述公式做一些修改。在這里，僅僅以p（x，y）二次多項(xiàng)式為例。

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ù)，按照最小二乘法的原則，應(yīng)該是

?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)點(diǎn)，（a,b）到離散數(shù)據(jù)點(diǎn)（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.編程實(shí)現(xiàn)

本文主要運(yù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;

注：＊代表指針變量；/代表注釋標(biāo)志。

4.結(jié)語(yǔ)

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

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

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

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