兩種空間插值方法的比較研究

1、兩種空間插值方法的比較研究?jī)煞N空間插值方法的比較研究摘要：距離倒數(shù)加權(quán)法算法簡(jiǎn)單，容易實(shí)現(xiàn)，適合分布較均勻的采樣點(diǎn)集，但容易出現(xiàn)“牛眼”現(xiàn)象；克里金法是一種無(wú)偏最優(yōu)估計(jì)法，精度較高，適合空間自相關(guān)程度高的數(shù)據(jù)，但其算法復(fù)雜，實(shí)現(xiàn)較難。這兩種方法各有其適用情形，本文比較了這兩種方法的優(yōu)劣并提出算法優(yōu)化的思路。關(guān)鍵字：距離倒數(shù)加權(quán)，克里金，優(yōu)化1引言空間插值是根據(jù)一組已知的離散數(shù)據(jù)或分區(qū)數(shù)據(jù)，按照某種假設(shè)推求出其他未知點(diǎn)或未知區(qū)域的數(shù)據(jù)的過(guò)程，簡(jiǎn)單的說(shuō)就是由已知空間特性推求未知空間特性。它是地學(xué)研究中的基本問(wèn)題，也是GIS 數(shù)據(jù)處理的重要內(nèi)容。在利用GIS 處理空間數(shù)據(jù)的過(guò)程中，需要進(jìn)行空間插值的

2、場(chǎng)合很多，如采樣密度不夠、采樣分布不合理、采樣存在空白區(qū)、等值線的自動(dòng)繪制、數(shù)字高程模型的建立、區(qū)域邊界分析、曲線光滑處理、空間趨勢(shì)預(yù)測(cè)、采樣結(jié)果的2.5維可視化等1。通過(guò)歸納，空間插值可以簡(jiǎn)化為以下三種情形：（1）現(xiàn)有離散曲面的分辨率、像元大小或方向與所要求的不符，需要重新插值。例如將一個(gè)掃描影像（航空像片、遙感影像）從一種分辨率或方向轉(zhuǎn)換為另一種分辨率或方向的影像。（2）現(xiàn)有連續(xù)曲面的數(shù)據(jù)模型與所需的數(shù)據(jù)模型不符，需要重新插值。如將一個(gè)連續(xù)曲面從一種空間切分方式變?yōu)榱硪环N空間切分方式，從TIN 到柵格、柵格到TIN 或矢量多邊形到柵格。（3）現(xiàn)有數(shù)據(jù)不能完全覆蓋所要求的區(qū)域范圍，需要插值。

3、如將離散的采樣點(diǎn)數(shù)據(jù)內(nèi)插為連續(xù)的數(shù)據(jù)表面2。 現(xiàn)有的空間插值方法多種多樣，但每一種方法都有其適用情形和無(wú)法避免的缺陷，本文分析了距離倒數(shù)加權(quán)法和克里金法的插值結(jié)果，并提出改進(jìn)的思路。 2方法距離倒數(shù)加權(quán)法和克里金法都是建立在地理學(xué)第一定律之上的，即：空間距離越近，地理事物的相似性越大3。它們都是通過(guò)確定待插點(diǎn)周圍采樣點(diǎn)的權(quán)重來(lái)求取待插點(diǎn)的估計(jì)值，可統(tǒng)一表示。設(shè)n x x ,1 為區(qū)域上的一系列觀測(cè)點(diǎn)，)(,),(1n x Z x Z 為相應(yīng)的觀測(cè)值。待插點(diǎn)0x 處的值)(0x Z 可采用一個(gè)線性組合來(lái)估計(jì)：=ni i i x Z x Z 10)()( （1）但距離倒數(shù)加權(quán)法只考慮采樣點(diǎn)與待插點(diǎn)

4、之間的距離，而克里金法不僅考慮距離，還要考慮 采樣點(diǎn)的空間分布及其與待插點(diǎn)的空間方位關(guān)系4。2.1距離倒數(shù)加權(quán)法距離倒數(shù)加權(quán)權(quán)重i 的賦值表達(dá)式為0;,2,1),(),(100=-m i x x d x x d m i ii i （2）式中，冪指數(shù)越小，權(quán)重越趨向取平均值；越大，越近的點(diǎn)權(quán)重越大，越遠(yuǎn)的點(diǎn)權(quán)重越小5。當(dāng)為零時(shí)，就是等權(quán)模型，即m i /1=，等權(quán)雖然簡(jiǎn)單易操作，但忽略了地理學(xué)第一定律。有的文獻(xiàn)6采用下列方案確定權(quán)值i ?=else x x d x x d x x d x x d n i i ,0),(),.,(),(min ),(,121 （3）它相當(dāng)于待插點(diǎn)取最鄰近點(diǎn)的值，即

5、泰森多邊形法（最近鄰點(diǎn)插值法）。2.2克里金法地質(zhì)統(tǒng)計(jì)學(xué)是以區(qū)域化變量為基礎(chǔ)，借助變異函數(shù)，研究既具有隨機(jī)性又具有結(jié)構(gòu)性，或空間相關(guān)性和依賴性的自然現(xiàn)象的一門(mén)科學(xué)4?？死锝鸩逯捣ㄊ堑刭|(zhì)統(tǒng)計(jì)學(xué)的重要組成部分，也是地質(zhì)統(tǒng)計(jì)學(xué)的核心。2.2.1區(qū)域化變量能用空間分布來(lái)表征一個(gè)自然現(xiàn)象的變量稱為區(qū)域化變量，它反映了區(qū)域內(nèi)的某種特征或現(xiàn)象。區(qū)域化變量根據(jù)區(qū)域內(nèi)位置的不同而取不同的值，可以說(shuō)，它是與位置有關(guān)的隨機(jī)變量。在進(jìn)行采樣觀測(cè)以后，可以將其表示為一個(gè)空間點(diǎn)函數(shù)),()(w v u x x x Z x Z = （4）式中， w v u x x x , 為三維直角坐標(biāo)系中的三軸。區(qū)域化變量具有以下幾個(gè)特

6、性：隨機(jī)性：區(qū)域化變量是一個(gè)隨機(jī)變量，它具有局部的、隨機(jī)的、異常的特征；結(jié)構(gòu)性：區(qū)域化變量在一定范圍內(nèi)具有某種程度的相似性，即自相關(guān)性，當(dāng)超出這一范圍時(shí)，自相關(guān)性消失；空間局限性：即這種結(jié)構(gòu)性的表現(xiàn)被限定在一定的空間內(nèi)；連續(xù)性：不同的區(qū)域化變量具有不同程度的連續(xù)性，其連續(xù)性是有變異函數(shù)來(lái)表示的；異向性：區(qū)域化變量可能表現(xiàn)為各向同性，也可能表現(xiàn)為各向異性4,7。2.2.2平穩(wěn)假設(shè)克里金插值法是一種無(wú)偏最優(yōu)估計(jì)法8，其中，無(wú)偏是指偏差為0，即要服從二階平穩(wěn)或本證假設(shè)。1.二階平穩(wěn)當(dāng)區(qū)域化變量)(x Z 滿足下列兩個(gè)條件時(shí)，稱其為二階平穩(wěn)或弱平穩(wěn)，在整個(gè)研究區(qū)內(nèi)有)(x Z 的數(shù)學(xué)期望存在，且等于常

7、數(shù)，即：m h x Z E x Z E =+=)()( （5）在整個(gè)研究區(qū)內(nèi)，)(x Z 的協(xié)方差函數(shù)存在且平穩(wěn)，即只依賴于滯后h ，而與x 無(wú)關(guān)：)()()()()()()()(),(2h C mh x Z x Z E h x Z E x Z E h x Z x Z E h x Z x Z Cov =-+=+-+=+ （6）特殊的，當(dāng)0=h 時(shí)，上式變?yōu)?0()(C x Z Var =，即方差存在且為常數(shù)。2.本征假設(shè)是比二階平穩(wěn)更弱的平穩(wěn)假設(shè)，當(dāng)區(qū)域化變量)(x Z 的增量)()(h x Z x Z +-滿足下列兩條件時(shí)，稱其為滿足本征假設(shè)或內(nèi)蘊(yùn)假設(shè)。在整個(gè)研究區(qū)內(nèi)有0)()(=+-h x

8、 Z x Z E （7）增量)()(h x Z x Z +-的方差函數(shù)存在且平穩(wěn)（即不依賴于x ）：)(2),(2)()()()()()()()(222h h x h x Z x Z E h x Z x Z E h x Z x Z E h x Z x Z Var =+-=+-+-=+- （8）2.2.3變異函數(shù)變異函數(shù)是地統(tǒng)計(jì)學(xué)特有的基本工具，它既能描述區(qū)域化變量的空間結(jié)構(gòu)性變化,又能描述其隨機(jī)性變化。區(qū)域化變量)(x Z 在點(diǎn)x 和h x +處的值)(x Z 與)(h x Z +差的方差的一半稱為區(qū)域化變量)(x Z 的變異函數(shù)，記為),(h x 。在二階平穩(wěn)假設(shè)或本征假設(shè)的條件下，有222

9、)()(21)()()()(21)()(21),(h x Z x Z E h x Z E x Z E h x Z x Z E h x Z x Z Var h x +-=+-+-=+-= （9） 由上式可知，變異函數(shù)依賴于自變量x 和h ，當(dāng)變異函數(shù)),(h x 僅僅依賴于距離h 而與位置x 無(wú)關(guān)時(shí)，),(h x 可改寫(xiě)為)(h ，即2)()(21)(h x Z x Z E h +-= （10） 具體表示為 =+-=)(12)()()(21)(h N i i i h x Z x Z h N h （11）在變異函數(shù)中有四個(gè)參數(shù)，它們的定義如下。變程（range ）：指區(qū)域化變量在空間上具有相關(guān)性的

10、范圍。在變程范圍之內(nèi)，數(shù)據(jù)具有相關(guān)性；而在變程之外，數(shù)據(jù)之間互不相關(guān)，即在變程以外的觀測(cè)值不對(duì)估計(jì)結(jié)果產(chǎn)生影響。塊金值(nugget)：變異函數(shù)如果在原點(diǎn)間斷，在地質(zhì)統(tǒng)計(jì)學(xué)中稱為“塊金效應(yīng)”，表現(xiàn)為在很短的距離內(nèi)有較大的空間變異性，無(wú)論h 多小，兩個(gè)隨機(jī)變量都不相關(guān)。測(cè)量誤差與自然現(xiàn)象的微觀變異性任意一方或兩者共同作用產(chǎn)生了塊金值?；_(tái)值(sill)：代表變量在空間上的總變異性大小。即為變差函數(shù)在h 大于變程時(shí)的值，為塊金值和拱高之和。拱高(partial sill)：基臺(tái)值與塊金值之差。當(dāng)塊金值等于0時(shí)，基臺(tái)值即為拱高。當(dāng)變異函數(shù)確定之后，執(zhí)行克里金系統(tǒng)就只是一個(gè)簡(jiǎn)單的計(jì)算過(guò)程。變異函數(shù)的理

11、論模型主要有以下三個(gè)9,10。球狀模型()? ?-?=? ?=ah c a h a h a h c h a h Sph c h ,2123003 （12） 指數(shù)模型 ()? ?-?=? ?=a h c a h Exp c h 3exp 1 （13） 高斯模型 ()()? ?-?=223exp 1a h c h （14） 式中，c 為基臺(tái)值，a 為變程，h 為滯后距變異函數(shù)的形式是內(nèi)插質(zhì)量的關(guān)鍵，需要注意的是，由于不同的區(qū)域有不同的空間模式，因而也就有不同的變異函數(shù)，而空間內(nèi)插都有一個(gè)隱含的假定，即空間是連續(xù)的，因此，在選擇變異函數(shù)模型之前，檢查數(shù)據(jù)以確定空間連續(xù)性是十分必要的11。2.2.4普

12、通克里金依據(jù)本征假設(shè)和無(wú)偏性要求，有=n i i 11（15）在以上的條件下要達(dá)到估計(jì)方差最小10，即()()()()()()()()()min 200*200*00*2=-=-=x Z x Z E x Z x Z E x Z x Z E k （16）再應(yīng)用拉格朗日乘數(shù)法求極值，進(jìn)一步推導(dǎo)可得()()()?=-=-=n i i j ni i j i n j x x C x x C 1011,1 （17）最小的估計(jì)方差,即克里金方差可用以下公式求解 ()()=-+-=-=n i i i k x x C x x C x Z x Z Var 100000*2)()( （18）或()()00102x

13、x x x ni i i k -+-= （19） 3實(shí)驗(yàn)3.1實(shí)驗(yàn)數(shù)據(jù)福建地處我國(guó)東南部，山地、丘陵占陸域的80%，屬溫暖濕潤(rùn)的亞熱帶海洋性季風(fēng)氣候12，其地形復(fù)雜，降雨豐富，本文選取福建省為空間插值的實(shí)驗(yàn)區(qū)，收集該省67個(gè)氣象站點(diǎn)的歷史數(shù)據(jù)，以2005年12月月平均降水為試驗(yàn)指標(biāo)，進(jìn)行內(nèi)插分析比較研究。3.2距離倒數(shù)加權(quán)法應(yīng)用ArcGIS9.2 Geostatistcal Analyst 功能模塊進(jìn)行插值，搜索半徑取1.2，冪指數(shù)取2， 從CrossValidation 的統(tǒng)計(jì)來(lái)看，最大誤差高達(dá)-200.667774109（武夷山），最小誤差只有-1龍海），出現(xiàn)這一

14、結(jié)果的原因是在武夷山站周圍氣象觀測(cè)站點(diǎn)分布較稀疏，可選點(diǎn)較少；而在龍海站周圍的分布較均勻，可選點(diǎn)較多。在結(jié)果圖中“牛眼”現(xiàn)象比較明顯的是九仙山站（269）（見(jiàn)圖1），與它相鄰的站點(diǎn)分別是永春站（120）、大田站（182）、尤溪站（182）、永泰站（124），產(chǎn)生這一現(xiàn)象的原因是距離倒數(shù)加權(quán)法在站點(diǎn)數(shù)據(jù)較少時(shí)容易受采樣極值點(diǎn)的影響，會(huì)在網(wǎng)絡(luò)區(qū)域產(chǎn)生圍繞采樣點(diǎn)位置的“靶心”，即所謂的“牛眼”現(xiàn)象13。圖1 距離倒數(shù)加權(quán)法，九仙山站及其周圍四縣市的插值結(jié)果圖3.3克里金法在克里金精度評(píng)價(jià)系統(tǒng)中，符合下列標(biāo)準(zhǔn)的結(jié)果是最優(yōu)的：標(biāo)準(zhǔn)平均值（Mean Standardized）最接近于0，均方根預(yù)測(cè)誤差（R

15、oot-Mean-Square）最小，平均標(biāo)準(zhǔn)誤差（Average Mean Error）最接近于均方根預(yù)測(cè)誤差（Root-Mean-Square），標(biāo)準(zhǔn)均方根預(yù)測(cè)誤差（Root-Mean-Square Standardized）最接近于14。本次插值結(jié)果的各參數(shù)如下：表1 預(yù)測(cè)誤差（prediction errors）標(biāo)準(zhǔn)平均值（MS）-0.001175均方根（RMS）52.92平均標(biāo)準(zhǔn)誤差（AME）43.94標(biāo)準(zhǔn)均方根（RMSS）0.9058從表中各數(shù)據(jù)來(lái)看，本次插值精度一般，這主要是因?yàn)椴蓸狱c(diǎn)較稀疏且自相關(guān)程度較低造成的，但從插值結(jié)果圖來(lái)看基本消除了“牛眼”現(xiàn)象（見(jiàn)圖2），插值精度較距離

16、倒數(shù)加權(quán)法有所提高，總的來(lái)看，基本符合降雨量西北高，東南低的特征。從CrossV alidation的統(tǒng)計(jì)來(lái)看，最大誤差仍在武夷山站（-137.071350687），最小誤差為1.0930*（漳州）。從點(diǎn)的權(quán)值來(lái)看，對(duì)待插點(diǎn)影響較大的點(diǎn)并不多，完全可以通過(guò)對(duì)選點(diǎn)方式的改進(jìn)來(lái)減少計(jì)算量。此外，選點(diǎn)過(guò)多時(shí)會(huì)出現(xiàn)負(fù)權(quán)值的情況，Clayton V Deutsch提出了一種校正負(fù)權(quán)值的方法8。圖2 克里金法，九仙山站及其周圍四縣市的插值結(jié)果圖4分析從以上的結(jié)果分析來(lái)看，每一種空間插值方法都有其適用的條件和各自的優(yōu)缺點(diǎn)，沒(méi)有一種方法具有放之四海而皆準(zhǔn)的普適性。距離倒數(shù)加權(quán)法算法簡(jiǎn)單，容易實(shí)現(xiàn)，當(dāng)采樣點(diǎn)較密

17、集且均勻分布時(shí)精度較高，但是在采樣點(diǎn)稀疏時(shí)則容易產(chǎn)生“牛眼”現(xiàn)象，而且冪指數(shù)的選取對(duì)采樣點(diǎn)的權(quán)重賦值有很大的影響，在全局僅采用一個(gè)冪指數(shù)是不符合實(shí)際情形的，George Y Lu，David W Wong對(duì)此進(jìn)行了改進(jìn)，提出一種用點(diǎn)模式分布函數(shù)來(lái)確定值的方法5?？死锝鸱ㄊ且环N無(wú)偏估計(jì)方法，在采樣點(diǎn)密集且具有高度空間自相關(guān)的情況下插值效果非常好，但它算法復(fù)雜，計(jì)算量大，對(duì)數(shù)據(jù)的要求很高，而且在變異函數(shù)的各參數(shù)設(shè)置和模型選取方面，對(duì)使用者的經(jīng)驗(yàn)要求比較高。基于對(duì)各插值方法的原理及其插值結(jié)果的思考，本文就插值算法的優(yōu)化提出以下兩點(diǎn)建議。(1)算法自身的優(yōu)化不管是距離倒數(shù)加權(quán)法，還是克里金法，它們選取采樣點(diǎn)的方式都