顏金彪，吳 波，何清華

1. 江西師范大學(xué)地理與環(huán)境學(xué)院，江西 南昌 330022； 2. 衡陽(yáng)師范學(xué)院傳統(tǒng)村鎮(zhèn)文化數(shù)字化保護(hù)與創(chuàng)意利用技術(shù)國(guó)家地方聯(lián)合工程實(shí)驗(yàn)室，湖南 衡陽(yáng) 421002

在地學(xué)研究中，通常需要具有一定規(guī)則的面狀數(shù)據(jù)。然而，很多可獲得的地學(xué)數(shù)據(jù)都是站點(diǎn)采集的，呈現(xiàn)離散和不規(guī)則的分布[1]。這就經(jīng)常需要利用空間插值技術(shù)將不規(guī)則的離散數(shù)據(jù)轉(zhuǎn)換成規(guī)則的面狀空間數(shù)據(jù)。

空間插值是通過(guò)有限的離散采樣點(diǎn)來(lái)建立某種插值函數(shù)關(guān)系，并將待插值點(diǎn)一定范圍內(nèi)的已知采樣點(diǎn)代入該函數(shù)表達(dá)式來(lái)獲得該插值點(diǎn)的屬性值[2]。常用的插值方法包括克里金插值方法、自然鄰域法、趨勢(shì)面法及反距離加權(quán)插值(inverse distance weighting,IDW)等。由于IDW插值原理簡(jiǎn)單、計(jì)算簡(jiǎn)便且符合地理學(xué)第一定律[3]，在各學(xué)科領(lǐng)域都得到了廣泛應(yīng)用[4-7]，被認(rèn)為是GIS領(lǐng)域中一種標(biāo)準(zhǔn)的空間插值方法。盡管在多數(shù)情況下IDW插值方法都能取得較好的效果，但由于地學(xué)數(shù)據(jù)的復(fù)雜性，其插值效果還有待進(jìn)一步提高。

近年來(lái)，人們從最鄰近點(diǎn)數(shù)、距離衰減系數(shù)以及空間鄰近度[8]等幾方面提出了各種改進(jìn)算法，改善了IDW空間插值精度及其適用性。在如何選擇最鄰近點(diǎn)數(shù)方面，文獻(xiàn)[9]認(rèn)為利用14個(gè)最鄰近的點(diǎn)進(jìn)行IDW插值是合適的，而文獻(xiàn)[10]認(rèn)為用5個(gè)最鄰近點(diǎn)就可以得到滿(mǎn)意的結(jié)果。其次，文獻(xiàn)[11]在評(píng)估城市噪聲地圖時(shí)發(fā)現(xiàn)距離衰減系數(shù)α=4比α=1能獲得更加準(zhǔn)確的噪聲分布地圖；文獻(xiàn)[9]認(rèn)為不同的α值將有不同的插值結(jié)果，并建議通過(guò)多次試驗(yàn)來(lái)選擇參數(shù)α。文獻(xiàn)[10]認(rèn)為在研究區(qū)內(nèi)的不同地點(diǎn)可能存在不同的距離衰減關(guān)系，并開(kāi)發(fā)了一種自適應(yīng)IDW方法(AIDW)，允許參數(shù)α根據(jù)鄰近采樣點(diǎn)的空間模式變化。在算法參數(shù)的優(yōu)化方面，現(xiàn)有方法或者僅考慮單參數(shù)對(duì)IDW插值的影響，或者依次對(duì)各參數(shù)進(jìn)行獨(dú)立優(yōu)化，都難實(shí)現(xiàn)模型插值精度的整體優(yōu)化。文獻(xiàn)[7]固定距離衰減系數(shù)為2，利用交叉檢驗(yàn)(cross validation，CV)[12]去優(yōu)化得到最鄰近點(diǎn)數(shù)；文獻(xiàn)[9]固定最鄰近點(diǎn)數(shù)為14，利用CV去優(yōu)化距離衰減系數(shù)，但候選的距離衰減系數(shù)僅有3個(gè)可選值?？紤]到IDW的插值精度受到多個(gè)參數(shù)的影響，部分參數(shù)獨(dú)立調(diào)優(yōu)難實(shí)現(xiàn)插值模型的整體優(yōu)化，因此，實(shí)踐上需要提供一種智能化程度較高的多參數(shù)協(xié)同優(yōu)化方法。

經(jīng)典的IDW假設(shè)地理事物的分布呈現(xiàn)各向同性[13]的特點(diǎn)，因而通常采用歐氏距離來(lái)表達(dá)空間鄰近度。由于地理事物分布的復(fù)雜性，其分布往往呈現(xiàn)各向異性[14]的特點(diǎn)，如氣溫呈現(xiàn)出明顯的經(jīng)緯度方向性，使得各向同性的距離度量難以反映氣溫在經(jīng)緯度方向的不同變化。文獻(xiàn)[15]將歐氏距離轉(zhuǎn)變?yōu)榇逯迭c(diǎn)與參考樣本點(diǎn)的最短路徑距離來(lái)重新定義空間鄰近度；文獻(xiàn)[16]提出采用一種“橢圓”距離來(lái)處理各向異性的插值問(wèn)題，文獻(xiàn)[17]提出了一種基于最短風(fēng)場(chǎng)路徑的復(fù)雜距離來(lái)表征城市風(fēng)場(chǎng)對(duì)顆粒物濃度分布的影響；針對(duì)河道各向異性的復(fù)雜地形，文獻(xiàn)[18]采用最短時(shí)間距離來(lái)替代歐氏距離；顧及電離層分布的經(jīng)緯度方向異性，文獻(xiàn)[19]重新定義空間鄰近性后獲得了更高精度的電離層插值效果。目前已有的研究大多情況針對(duì)特定應(yīng)用場(chǎng)景來(lái)改進(jìn)算法，并通過(guò)重新定義空間鄰近度來(lái)提高IDW的插值精度，需要額外的輔助數(shù)據(jù)(如風(fēng)速、流速等)來(lái)定義各向異性的空間鄰近度，制約了這類(lèi)方法的普適性應(yīng)用。

本文提出一種顧及空間各向異性的多參數(shù)協(xié)同優(yōu)化IDW插值方法(PIDW)。首先，引入具有物理意義的距離及方向調(diào)節(jié)參數(shù)來(lái)反映地理現(xiàn)象的各向異性，將典型的歐氏距離擴(kuò)展為各向異性的“橢圓”距離；然后，基于CV準(zhǔn)則構(gòu)造了適合于PIDW插值算法的適應(yīng)度函數(shù)，并采用粒子群算法對(duì)PIDW插值模型的多個(gè)參數(shù)進(jìn)行協(xié)同優(yōu)化，實(shí)現(xiàn)PIDW插值效果最佳的目的。

1 PIDW算法

PIDW是一種顧及各向異性的IDW插值的改進(jìn)模型，通過(guò)引入距離及方向調(diào)節(jié)參數(shù)來(lái)描述空間點(diǎn)間各向異性的鄰近性關(guān)系，并利用粒子群優(yōu)化算法實(shí)現(xiàn)模型的多參數(shù)協(xié)同優(yōu)化。圖1是PIDW方法的主要流程，其算法主要包括四參數(shù)協(xié)同優(yōu)化以及插值兩大過(guò)程。

圖1 PIDW方法的流程Fig.1 The flowchart of the proposed PIDW method

1.1 各向異性表達(dá)模型

IDW插值模型如式(1)所示

(1)

式中，i代表任意待插值點(diǎn)；zj為點(diǎn)j處的實(shí)測(cè)值；Sij是未知點(diǎn)i與已知點(diǎn)j之間的距離；N為最鄰近點(diǎn)數(shù)；α為距離衰減系數(shù)。

(2)

式中，(Xi,Yi)與(Xj,Yj)表示點(diǎn)i與j的笛卡兒坐標(biāo)。

由式(2)可知，當(dāng)λ等于1，那么Sij將等價(jià)于歐氏空間距離；當(dāng)λ>1，如果兩點(diǎn)之間ΔY越大，那么兩點(diǎn)Y方向之間的距離將更遠(yuǎn)，即Y方向的空間鄰近性減弱；反之，當(dāng)λ<1，那么兩點(diǎn)之間X方向的距離相對(duì)將增大，即兩點(diǎn)的X方向的空間鄰近性將減弱。如圖2所示，距離調(diào)節(jié)參數(shù)λ能夠改變最鄰近點(diǎn)的選擇，其中圖2(a)為均質(zhì)空間下所選擇的空間鄰近點(diǎn)，圖2(b)為顧及各向異性時(shí)選擇的空間鄰近點(diǎn)。

圖2 各向同性與各向異性環(huán)境下最鄰近點(diǎn)的選擇Fig.2 Illustration of the different selection of proximity points in the cases of isotropic and anisotropic assumptions

式(2)蘊(yùn)含了數(shù)據(jù)的各向異性方向正好沿著數(shù)據(jù)的坐標(biāo)軸，這與大多數(shù)的實(shí)際情況不符。因此，需要進(jìn)一步引入方向參數(shù)θ來(lái)描述一般情形下的各向異性，從而可以建立原始坐標(biāo)與反映各向異性坐標(biāo)的轉(zhuǎn)換關(guān)系，如式(3)所示

(3)

結(jié)合式(2)、式(3)，可以得到式(4)

(4)

式中，Eij表示待插值點(diǎn)i與樣本點(diǎn)j之間顧及各向異性的距離；λ用于調(diào)節(jié)坐標(biāo)系統(tǒng)下不同坐標(biāo)軸的距離；θ∈(0,2π)表示原始坐標(biāo)與各向異性坐標(biāo)系的旋轉(zhuǎn)角度；ΔX與ΔY為原始數(shù)據(jù)中待插值點(diǎn)與已知樣本點(diǎn)之間的歐氏距離。

1.2 基于粒子群的多參數(shù)協(xié)同優(yōu)化

考慮到影響IDW插值效果由距離衰減系數(shù)、最鄰近點(diǎn)數(shù)、距離調(diào)節(jié)及方向參數(shù)共同決定，參數(shù)的確定本質(zhì)上是一個(gè)多元變量的尋優(yōu)過(guò)程。本文提出采用粒子群智能優(yōu)化搜索算法[20]，通過(guò)CV方法構(gòu)造粒子群優(yōu)化的適應(yīng)度函數(shù)，將多變量的優(yōu)化問(wèn)題轉(zhuǎn)化為求解函數(shù)的極值問(wèn)題，從而實(shí)現(xiàn)插值精度偏差與方差的總體滿(mǎn)意效果。

1.2.1 粒子群算法

粒子群算法采用群體智能優(yōu)化策略。首先，給空間中所有粒子分配初始隨機(jī)位置和速度；然后，根據(jù)每個(gè)粒子的速度、問(wèn)題空間中已知的全局最優(yōu)位置和粒子已知的最優(yōu)位置依次更新每個(gè)粒子的位置，通過(guò)種群內(nèi)粒子間的合作與競(jìng)爭(zhēng)機(jī)制來(lái)進(jìn)行全局優(yōu)化獲得最優(yōu)解。在D維空間中，設(shè)定粒子群數(shù)為Q，粒子i的位置表示為xi=(xi1,xi2,…,xiD)，粒子i的速度為vi=(vi1,vi2,…,viD)，pbesti=(pi1,pi2,…,piD)表示粒子i經(jīng)歷過(guò)的最佳位置，gbest=(g1,g2,…,gD)表示所有粒子的全局最佳位置，f(xi)為適應(yīng)度函數(shù)[20]。其中粒子i的第d維位置和速度更新方法分別如式(5)、式(6)所示

(5)

(6)

1.2.2 適應(yīng)度函數(shù)

粒子群算法的多參數(shù)優(yōu)化問(wèn)題關(guān)鍵在于確定適應(yīng)度函數(shù)。為了降低可能的過(guò)擬合現(xiàn)象[21]，導(dǎo)致插值精度偏差與方差的背離，本文利用CV檢驗(yàn)方法來(lái)構(gòu)造PIDW算法的適應(yīng)度函數(shù)，如式(7)、式(8)所示

(7)

xi=argminf(xi)

(8)

1.2.3 算法流程

1.2.3.1 確定4個(gè)參數(shù)(α,λ,θ,N)的組合解

(1)給定粒子群算法的初始化條件，搜索維度D=4(α,λ,θ,N)，α、λ、θ、N的搜索區(qū)間分別為[1，5]、[0.01，100]、[0，2π]及[4，30]，其中參數(shù)N和α的搜索范圍是在參考文獻(xiàn)[9—11]的基礎(chǔ)上設(shè)定的，距離調(diào)節(jié)參數(shù)λ的范圍根據(jù)經(jīng)驗(yàn)確定，而角度參數(shù)θ的范圍為所有可能的取值。采用隨機(jī)初始化種群的方法產(chǎn)生一組包含(α,λ,θ,N)的解集合。

(2)計(jì)算粒子的適應(yīng)度值。將每個(gè)粒子的xi代入式(7)的適應(yīng)度函數(shù)，得到每個(gè)粒子的f(xi)。

(3)獲取粒子i的個(gè)體最優(yōu)值pbesti。對(duì)每個(gè)粒子，用其f(xi)值和個(gè)體最優(yōu)位置Pbest(i)比較，如果f(xi)優(yōu)于Pbest(i)，則用當(dāng)前待優(yōu)化參數(shù)的位置xi更新個(gè)體歷史最佳位置Pbest(i)。

(4)獲取群體的全局最優(yōu)值。對(duì)每個(gè)粒子，將其當(dāng)前f(xi)與全局最佳位置Gbest做比較，如果當(dāng)前的適應(yīng)值更佳，則用當(dāng)前粒子的位置更新全局最佳位置Gbest。

(5)根據(jù)式(5)和式(6)對(duì)粒子的速度和位置(包括距離衰減系數(shù)、最鄰近點(diǎn)數(shù)、距離調(diào)節(jié)參數(shù)以及各向異性方向)進(jìn)行更新。

(6)如未滿(mǎn)足粒子群算法的終止條件，則返回步驟(2)。本文設(shè)定迭代次數(shù)達(dá)到200或者最佳適應(yīng)度值的增量比值少于1e-6時(shí)算法停止。

(7)輸出最佳位置xi，即使得適應(yīng)度函數(shù)值最佳時(shí)的距離衰減系數(shù)、最鄰近點(diǎn)數(shù)、距離調(diào)節(jié)及各向異性方向參數(shù)。

1.2.3.2 插值

根據(jù)粒子群算法優(yōu)化得到的參數(shù)組合解(α,λ,θ,N)，代入IDW插值模型(式(1))中獲取待插值點(diǎn)屬性值。式(1)中Si的計(jì)算方式采用本文顧及各向異性方向的“橢圓”距離計(jì)算方式Eij(式(4))，其中α、N、λ、θ分別為優(yōu)化得到的距離衰減系數(shù)、最鄰近點(diǎn)數(shù)、距離調(diào)節(jié)及各向異性方向參數(shù)。

2 試 驗(yàn)

本文使用不同分辨率的全國(guó)氣溫站點(diǎn)與黃河堤壩數(shù)據(jù)來(lái)驗(yàn)證PIDW算法的可行性，并將PIDW插值算法與經(jīng)典的反距離平方加權(quán)(classical IDW,CIDW)、四象限IDW(four quadrant IDW,FIDW)[22]、自適應(yīng)的反距離平方權(quán)重(adaptive IDW,AIDW)[23]、K近鄰反距離平方加權(quán)(K-nearest neighbor adaptive IDW,KAIDW)[2]、普通克里金(ordinary Kriging，OK)以及顧及各向異性的克里金(anisotropic Kriging,AnisOK)進(jìn)行對(duì)比試驗(yàn)。

(9)

(10)

(11)

2.1 年均氣溫插值試驗(yàn)

數(shù)據(jù)來(lái)源于中國(guó)氣象中心2005年的682個(gè)氣象站點(diǎn)的觀測(cè)數(shù)據(jù)。通過(guò)前期的數(shù)據(jù)預(yù)處理，剔除少量觀測(cè)臺(tái)站儀器出現(xiàn)故障的數(shù)據(jù)，對(duì)每個(gè)臺(tái)站的日平均氣溫求和，最終獲取各個(gè)臺(tái)站的年平均氣溫。由于氣溫與海拔存在約0.6℃/100 m的比例關(guān)系，本文根據(jù)氣象站臺(tái)的海拔數(shù)據(jù)，將每個(gè)氣象站臺(tái)的年平均氣溫歸算到海平面。這樣處理的理由是剔除不同高程對(duì)年均氣溫的影響，著重討論經(jīng)緯度變化對(duì)氣溫差異的影響[26](圖3)。

圖3 氣象站點(diǎn)的分布Fig.3 Spatial distribution of national weather stations

為減少樣本選擇的影響，本文重復(fù)進(jìn)行了3次隨機(jī)抽樣，其中每次從總體中抽取30個(gè)作為測(cè)試集，其余作為訓(xùn)練集。圖3中“*”“·”及“+”形代表3次隨機(jī)抽樣的樣點(diǎn)分布。經(jīng)過(guò)PSO優(yōu)化，PIDW的距離衰減系數(shù)α為1，最鄰近點(diǎn)數(shù)N為7，距離調(diào)節(jié)參數(shù)λ為0.26，各向異性方向θ為0.02°。將計(jì)算結(jié)果的平均值作為插值精度的最終評(píng)價(jià)值。7種算法的對(duì)比結(jié)果如圖4所示。

圖4 利用MAE、RMSE及MAX指標(biāo)評(píng)價(jià)7種算法的插值精度Fig.4 Comparisons the accuracies of seven different algorithms using MAE,RMSE and MAX measures

從圖4(a)可以發(fā)現(xiàn)：CIDW與FIDW插值算法的MAE隨著插值點(diǎn)數(shù)的增多呈現(xiàn)震蕩式變化，并且FIDW的MAE顯著低于CIDW；KAIDW、AIDW與FIDW的MAE大致保持一致，但OK和PIDW的MAE顯著低于KAIDW、AIDW、FIDW和CIDW；OK的MAE低于PIDW，但最大偏差僅為0.007℃；AnisOK的MAE誤差為0.52，低于PIDW的0.54。由圖4(b)中可以看出：CIDW、FIDW和OK的RMSE隨著點(diǎn)數(shù)的增多而下降，在參考點(diǎn)數(shù)量為13以上時(shí)趨于穩(wěn)定；OK的RMSE稍好于AIDW與KAIDW方法，但劣于CIDW和FIDW算法；AnisOK的RMSE誤差顯著低于CIDW、FIDW、OK、AIDW、KAIDW；PIDW的RMSE為0.78℃，比次優(yōu)的AnisOK方法低0.02℃。圖4(c)中的MAX反映了與RMSE類(lèi)似的規(guī)律：CIDW與FIDW插值算法隨著點(diǎn)數(shù)的增多穩(wěn)步下降，其中FIDW的MAX略低于CIDW，而AIDW與KAIDW的MAX表現(xiàn)最差；OK的MAX總體上優(yōu)于AIDW和KAIDW，但明顯差于CIDW和FIDW算法；PIDW的MAX顯著低于CIDW、FIDW、AIDW、KAIDW以及OK；AnisOK的MAX具備一定的優(yōu)勢(shì)，相較于次優(yōu)的PIDW方法，其最大值誤差還要低0.4℃。

圖4的結(jié)果表明：①顧及各向異的PIDW和AnisOK都顯著優(yōu)于CIDW、FIDW、OK、AIDW和KAIDW算法；②不存在哪種插值算法，其MAE、RMSE與MAX的3個(gè)指標(biāo)上都占據(jù)一致性的優(yōu)勢(shì)；③顧及樣點(diǎn)分布均勻性的FIDW插值算法精度普遍優(yōu)于CIDW；④采用一階鄰近點(diǎn)作為參考樣點(diǎn)的AIDW與KAIDW插值算法并未體現(xiàn)出較大的優(yōu)勢(shì)，可能是由測(cè)試數(shù)據(jù)集中一些位于邊界區(qū)域的測(cè)試點(diǎn)所導(dǎo)致，顧及空間異質(zhì)性的KAIDW未呈現(xiàn)比AIDW插值算法的優(yōu)勢(shì)，可能是由于KAIDW對(duì)于局部空間異質(zhì)性的分類(lèi)存在較多的誤判；⑤OK相比于CIDW、FIDW、AIDW、KAIDW總體上并未呈現(xiàn)出明顯優(yōu)勢(shì)，這可能是由于研究區(qū)域不能較好地滿(mǎn)足OK插值算法所需的平穩(wěn)性假設(shè)條件，導(dǎo)致變異函數(shù)未能較好地反映氣溫變化的變異規(guī)律；⑥PIDW插值算法除在MAE指標(biāo)稍劣于OK插值算法，在其余指標(biāo)上相較于未顧及各向異性影響插值算法具備明顯的優(yōu)勢(shì)；⑦除RMSE指標(biāo)外，AnisOK在MAE和MAX指標(biāo)上低于PIDW插值算法。

2.2 DEM插值

DEM插值數(shù)據(jù)來(lái)源于黃河山東濟(jì)南段大堤的1∶500大比例尺地形圖，共采集了639個(gè)高程點(diǎn)。黃河大堤為“梯形”構(gòu)造，其中此段大堤南北方向起伏大，而東西向高程則變化平緩，其中平均高程為43.03 m，最大高程為49.30 m，位于大堤的頂部，最低高程36.31 m。

試驗(yàn)中各算法的初始值與氣溫?cái)?shù)據(jù)保持一致，PIDW算法得到最終的優(yōu)化參數(shù)結(jié)果為[1.28，6，99，0]，各種算法的結(jié)果對(duì)比如圖5所示。

圖5 對(duì)比7種算法在DEM數(shù)據(jù)插值的MAE、RMSE和MAX結(jié)果Fig.5 Comparisons of seven different algorithms in DEM interpolation with MAE, RMSE and MAX measures

從圖5(a)中可以發(fā)現(xiàn)：隨著插值點(diǎn)數(shù)的增多，CIDW、FIDW與OK的MAE逐漸下降，在采樣點(diǎn)為12時(shí)達(dá)到最低，分別為1.43、1.37和1.42 m；AIDW與KAIDW的MAE低于CIDW算法，但大于其余3種算法；AnisOK的MAE誤差明顯低于CIDW、FIDW、OK、AIDW及KAIDW；PIDW表現(xiàn)出明顯的優(yōu)勢(shì)，其MAE在7種插值算法中的誤差最小，比次優(yōu)AnisOK的最小MAE值再少27.3 cm。圖5(b)中的RMSE表現(xiàn)出類(lèi)似的結(jié)果：CIDW、FIDW和OK的RMSE隨著點(diǎn)數(shù)的增多呈現(xiàn)下降，在樣點(diǎn)為14以上時(shí)基本達(dá)到穩(wěn)定狀態(tài)，并且OK的RMSE值最??；AIDW與KAIDW的RMSE的結(jié)果最差；AnisOK的RMSE低于除PIDW外的5種插值方法；PIDW的RMSE指標(biāo)在所有算法中同樣呈現(xiàn)出明顯的優(yōu)勢(shì)，比次優(yōu)的AnisOK低27 cm左右。圖5(c)表明CIDW、FIDW與OK的MAX隨著點(diǎn)數(shù)增多而呈現(xiàn)下降趨勢(shì)，在參考點(diǎn)達(dá)到17以上基本保持穩(wěn)定；AIDW和KAIDW比CIDW與FIDW能夠取得更優(yōu)的效果，但PIDW的MAX顯著優(yōu)于其他插值方法?？偨Y(jié)以上兩個(gè)試驗(yàn)結(jié)果，可得出以下結(jié)論：①顧及各向異性的AnisOK與PIDW幾乎在所有指標(biāo)中都取得最小的誤差，表明在計(jì)算過(guò)程中顧及各向異性的空間屬性，將有助于提高模型的插值精度；②總體上PIDW的插值效果稍好于AnisOK方法，而且PIDW方法對(duì)數(shù)據(jù)無(wú)需假設(shè)前提，但地統(tǒng)計(jì)的AnisOK方法要求數(shù)據(jù)具有空間的二階平穩(wěn)性。

2.3 參數(shù)分析

為了解釋PIDW插值算法有效性的原因，下面對(duì)模型的參數(shù)優(yōu)化作進(jìn)一步試驗(yàn)分析。

2.3.1 各向異性對(duì)空間鄰近性的影響

由于各向異性能夠引起空間鄰近度的改變，因此設(shè)計(jì)對(duì)照試驗(yàn)來(lái)分析各向異性對(duì)空間鄰近度的影響。圖6(a)和圖6(b)分別表示未顧及各向異性和顧及各向異性時(shí)對(duì)空間鄰近度的影響。由圖6(a)可知，當(dāng)不顧及x或y方向的各向異性時(shí)，樣本點(diǎn)與參考點(diǎn)之間的距離為同心圓模式，而當(dāng)顧及各向異性的影響時(shí)，“同心圓”距離模式被壓縮為“橢圓”，位于東西方向的樣本點(diǎn)具有更強(qiáng)的空間鄰近性。這表明PIDW的距離調(diào)節(jié)參數(shù)符合堤壩高程南北方向起伏大，而東西向變化平緩的事實(shí)。

圖6 比較各向同性與各向異性環(huán)境下的空間鄰近度Fig.6 Comparisons the spatial proximity between isotropic and anisotropic measures

2.3.2 各向異性的方向?qū)臻g鄰近性的影響

考慮到各向異性的方向并不一定正好沿著坐標(biāo)軸方向，而是存在一定的角度θ，本文設(shè)計(jì)一組對(duì)照試驗(yàn)來(lái)分析各向異性的方向?qū)臻g鄰近性的影響。試驗(yàn)2的原始數(shù)據(jù)高程在南北方向變化速率比東西方向大得多。為了體現(xiàn)各向異性方向不沿x或y軸時(shí)對(duì)空間鄰近性造成的影響，將試驗(yàn)2的數(shù)據(jù)逆時(shí)針旋轉(zhuǎn)45°，即變換后的各向異性方向x_new與y_old之間的夾角約為135°(圖7)。

由圖7(a)可知，試驗(yàn)1得到的參考點(diǎn)與其余各樣本點(diǎn)之間的距離基本為同心圓模式，不能真實(shí)地反映高程變化的空間鄰近性，同時(shí)試驗(yàn)1得到的CV殘差和為2 894.4也再一次證明了其擬合精度較差。試驗(yàn)2得到的CV殘差為1 400.4，其計(jì)算得到的各向異性方向與原始y軸間的夾角為132.66°(圖7(b))。可以發(fā)現(xiàn)，試驗(yàn)2的CV殘差和較之于試驗(yàn)1低很多，同時(shí)試驗(yàn)2得到的各向異性方向與真實(shí)的方向相差僅為2.34°。同時(shí)由圖7(b)可知，正確反映各向異性的方向能夠得到準(zhǔn)確的空間鄰近關(guān)系。

圖7 各向異性方向?qū)臻g鄰近度度量的影響Fig.7 Illustrations of the anisotropic effects on spatial proximity measurement

2.3.3 多參數(shù)優(yōu)化對(duì)比分析

影響PIDW插值精度的參數(shù)主要有各向異性(距離調(diào)節(jié)參數(shù)λ與各向異性方向參數(shù)θ)，距離衰減系數(shù)α和最鄰近點(diǎn)數(shù)N。本文基于試驗(yàn)1的數(shù)據(jù)，設(shè)計(jì)如下5個(gè)試驗(yàn)來(lái)表明PIDW的插值精度由各參數(shù)共同決定，而單參數(shù)或少量參數(shù)的優(yōu)化并不能取得插值的整體滿(mǎn)意效果：①固定距離調(diào)節(jié)參數(shù)為1，方向值為0，優(yōu)化最鄰近點(diǎn)數(shù)與距離衰減系數(shù)；②固定距離衰減系數(shù)為2(ArcMap軟件默認(rèn)值)，優(yōu)化最鄰近點(diǎn)數(shù)、距離調(diào)節(jié)及各向異性方向參數(shù)；③固定最鄰近點(diǎn)數(shù)為12(ArcMap軟件默認(rèn)值)，優(yōu)化距離調(diào)節(jié)、各向異性方向參數(shù)及距離衰減系數(shù)；④固定方向值為0，同時(shí)優(yōu)化距離衰減系數(shù)、最鄰近點(diǎn)數(shù)及距離調(diào)節(jié)參數(shù)；⑤距離衰減系數(shù)、最鄰近點(diǎn)數(shù)、距離調(diào)節(jié)與各向異性方向參數(shù)協(xié)同優(yōu)化。試驗(yàn)采用MAE、RMSE和CV殘差作為精度的評(píng)價(jià)指標(biāo)，結(jié)果如圖8所示。

圖8 對(duì)比不同情形下參數(shù)優(yōu)化的插值精度Fig.8 comparison the interpolating accuracies with optimized parameters for the different cases

從圖8可以看出方案1的效果最差，其MAE、RMSE以及訓(xùn)練樣本的CV殘差值都最大。這是由于氣溫的分布存在較強(qiáng)的各向異性特征，如果插值過(guò)程中將距離調(diào)節(jié)參數(shù)固定為1，導(dǎo)致重要的參數(shù)不能正確地反映真實(shí)情況。方案2和方案3的MAE、RMSE和CV相較于方案1都有一定程度的提升，且方案3略?xún)?yōu)于方案2，表明本例中最鄰近點(diǎn)數(shù)對(duì)影響插值精度的重要性要低于距離衰減系數(shù)。方案4的CV殘差低于方案1、方案2和方案3，表明顧及各向異性影響能夠得到較好的插值效果。從圖8中還可以發(fā)現(xiàn)方案5的預(yù)測(cè)結(jié)果幾乎與方案4完全一致，僅有CV的殘差和相差0.1。這是由于氣溫變化沿著緯度方向變化速率較之于經(jīng)度更快，而方案5優(yōu)化的各向異性方向與方案4所固定的方向相同，因此其插值結(jié)果一致。由于粒子群算法的隨機(jī)性，導(dǎo)致它們的CV殘差和存在微小的差異。以上試驗(yàn)表明如果僅優(yōu)化部分參數(shù)均難以獲得插值效果的整體滿(mǎn)意解，采取多參數(shù)整體優(yōu)化方法將有效增強(qiáng)模型插值精度。由于PIDW算法也可能存在過(guò)擬合的現(xiàn)象，導(dǎo)致預(yù)測(cè)精度的降低，因而在實(shí)際應(yīng)用中需要注意避免PIDW算法的過(guò)擬合問(wèn)題。

3 結(jié) 論

本文提出了顧及各向異性的多參數(shù)協(xié)同優(yōu)化IDW插值算法PIDW，通過(guò)引入描述距離調(diào)節(jié)的參數(shù)λ以及角度的參數(shù)θ，將傳統(tǒng)的歐氏距離擴(kuò)展為任意方向的“橢圓”距離，從而顧及各向異性對(duì)插值結(jié)果的影響；構(gòu)造了基于CV檢驗(yàn)的粒子群優(yōu)化的適應(yīng)度函數(shù)，將影響IDW插值精度的4個(gè)關(guān)鍵參數(shù)進(jìn)行協(xié)同優(yōu)化。利用不同尺度的空間數(shù)據(jù)驗(yàn)證了PIDW算法的插值效果，結(jié)果表明顧及各向異性的AnisOK及PIDW能夠顯著提高各向異性環(huán)境下插值算法的精度，PIDW幾乎在所有指標(biāo)中優(yōu)于已有的IDW及OK插值方法。其次，PIDW總體插值效果稍好于AnisOK，并且不需要數(shù)據(jù)滿(mǎn)足二階平穩(wěn)性假設(shè)條件，因此PIDW比AnisOK更具應(yīng)用的普適性。

PIDW插值算法獲得的協(xié)同優(yōu)化參數(shù)是一種全局意義下的滿(mǎn)意解。然而，當(dāng)局部地理現(xiàn)象的變化特征與整體變化規(guī)律不一致時(shí)，那么PIDW算法得到的局部插值結(jié)果可能不如傳統(tǒng)的IDW插值方法；其次，由于采用粒子群算法對(duì)影響IDW插值精度的多個(gè)參數(shù)同步優(yōu)化，處理效率比傳統(tǒng)IDW方法低。今后將在以下兩個(gè)方面開(kāi)展研究：①利用GPU加速或并行處理技術(shù)提升PIDW插值算法的效率[27]；②研究能識(shí)別局部各向異性特征的LPIDW插值算法。