李 豪，劉 濤，徐精文

?

基于混合地理加權(quán)回歸與克里格的區(qū)域降水量空間插值方法*

李 豪，劉 濤，徐精文

（四川農(nóng)業(yè)大學(xué)資源學(xué)院，成都 611130）

基于四川省區(qū)域范圍內(nèi)144個氣象站點(diǎn)的實(shí)測降水?dāng)?shù)據(jù)，在綜合考慮空間位置、地形等影響因素的基礎(chǔ)上，采用改進(jìn)的回歸克里格模型，即混合地理加權(quán)回歸克里格模型（MGWRK）對四川省年降水量的空間分布進(jìn)行空間插值，并與普通克里格（OK）、全局回歸克里格（GRK）和地理加權(quán)回歸克里格（GWRK）等模型的插值效果進(jìn)行對比分析。結(jié)果表明：（1）應(yīng)用逐步回歸法篩選確定的用于回歸分析的影響因子組合為經(jīng)度、緯度和坡度，可有效消除解釋變量間的多重共線性，為后續(xù)的空間插值奠定基礎(chǔ)；（2）同一回歸變量在地理加權(quán)回歸（GWR）與全局回歸（GR）兩種回歸模型中的AICc（修正的赤池信息量準(zhǔn)則，Corrected Akaike Information Criterion）值之差（ΔAICc）可用于定量判定各回歸變量的空間非平穩(wěn)性類型，據(jù)此將變量坡度設(shè)為全局變量，經(jīng)度和緯度設(shè)為局部變量進(jìn)行處理。在此基礎(chǔ)上，通過MGWRK模型對四川省年降水量進(jìn)行空間插值；（3）MGWRK插值模型綜合考慮了空間位置、地形等多個影響因素及其與降水相互關(guān)系的空間非平穩(wěn)性特征，相對于傳統(tǒng)的OK和GRK法具有更高的插值精度。

降水量；混合地理加權(quán)回歸；克里格；空間插值

降水是氣象、農(nóng)業(yè)和環(huán)境等領(lǐng)域研究的重要數(shù)據(jù)源，對農(nóng)業(yè)生產(chǎn)、區(qū)域水資源管理、防洪減災(zāi)、生態(tài)環(huán)境治理等工作均具有重要的應(yīng)用價值。近年來，隨著眾多領(lǐng)域研究工作的深入，對區(qū)域降水?dāng)?shù)據(jù)在空間尺度上的精細(xì)化需求也逐步提升。

在空間尺度上，目前獲取區(qū)域降水?dāng)?shù)據(jù)的方法主要有兩種：一是氣象站的觀測[1]，這是現(xiàn)階段最主要、技術(shù)最成熟的手段，具有測量精度高、穩(wěn)定可靠等優(yōu)點(diǎn)。但該方法獲得的數(shù)據(jù)為離散的點(diǎn)數(shù)據(jù)，無法提供區(qū)域內(nèi)任意位置的降水資料，且由于受到地形條件及資金投入等的限制，現(xiàn)階段仍存在站點(diǎn)間距離較大、密度偏低等不足，由此獲取的降水?dāng)?shù)據(jù)往往難以滿足小尺度研究對數(shù)據(jù)精細(xì)度的需求。另一種方法是衛(wèi)星遙感降水?dāng)?shù)據(jù)。該方法不直接測量降水，而是通過衛(wèi)星上搭載的傳感器測定云頂紅外溫度、雷達(dá)反射率等數(shù)據(jù)，再反演獲取降水?dāng)?shù)據(jù)[2]，TRMM、GSMaP和GPCP等降水?dāng)?shù)據(jù)集均屬此類。該方法獲得的降水?dāng)?shù)據(jù)為連續(xù)的面數(shù)據(jù)，但大多著眼于全球尺度，空間分辨率相對較低，如TRMM數(shù)據(jù)的分辨率僅為0.25°×0.25°，一些局部、細(xì)節(jié)性的降水空間分布特征未能表達(dá)出來，且存在不同區(qū)域反演精度差異大、覆蓋范圍有限等不足[2]，亦無法完全達(dá)到區(qū)域、流域或更小尺度研究對降水?dāng)?shù)據(jù)的精度要求。

以實(shí)測數(shù)據(jù)為基礎(chǔ)、基于GIS的空間插值技術(shù)已廣泛用于獲取氣象、土壤和環(huán)境等要素的高分辨率空間分布數(shù)據(jù)。該方法是氣象站觀測法的延伸，通過離散的氣象站實(shí)測數(shù)據(jù)預(yù)測未知區(qū)域的降水?dāng)?shù)據(jù)，可以快速獲取全面、高空間分辨率的區(qū)域降水資料，有效解決上述兩種方法的不足，是目前獲取區(qū)域降水空間數(shù)據(jù)的重要途徑，特別在具有充足地面觀測資料的地區(qū)具有較好的適用性和較高的估算精度[3?5]。

目前，用于降水?dāng)?shù)據(jù)的空間插值方法很多，不同方法的預(yù)測精度也不盡相同。根據(jù)已有數(shù)據(jù)與研究區(qū)域的空間分布特征，選擇一種最優(yōu)的插值方法是采用該方法獲取區(qū)域降水?dāng)?shù)據(jù)的關(guān)鍵。克里格（Kriging）法是基于地統(tǒng)計(jì)學(xué)的一種空間插值方法，它充分考慮了樣點(diǎn)的空間變異性特征，具有適用性強(qiáng)、預(yù)測精度高等優(yōu)點(diǎn)，目前在氣象、生態(tài)和土壤等領(lǐng)域的應(yīng)用最為廣泛[6?8]。已有研究表明，通過構(gòu)建回歸克里格（Regression Kriging，RK）模型，將克里格插值與回歸分析結(jié)合起來，引入與插值變量具有較好相關(guān)性的影響因子作為輔助變量，建立插值變量與其影響因子之間的回歸模型，可以有效提高克里格法的插值精度[9?12]。

本研究以地處長江上游、降水空間分異顯著的四川省為研究對象，在充分考慮降水量與各影響因子間相互關(guān)系及其空間非平穩(wěn)性特征的基礎(chǔ)上，嘗試采用一種改進(jìn)的回歸克里格插值模型，即混合地理加權(quán)回歸克里格模型（Mixed Geographically Weighted Regression Kriging，MGWRK）對四川省年降水量的空間分布進(jìn)行預(yù)測，并與全局回歸克里格、普通地理加權(quán)回歸克里格等方法的插值結(jié)果進(jìn)行對比分析，探尋適合本區(qū)域降水量的空間插值方法，從而準(zhǔn)確獲得四川省年降水量空間分布的精細(xì)化數(shù)據(jù)，以期為在該區(qū)域開展農(nóng)業(yè)產(chǎn)業(yè)結(jié)構(gòu)優(yōu)化、水資源配置、生態(tài)環(huán)境治理等工作提供有效數(shù)據(jù)支持。

1 資料與方法

1.1 研究區(qū)概況

四川省位于中國西南，地處長江上游，介于92°21′?108°12′E、26°03′?34°19′N，區(qū)域面積48.5萬km2[13]。本區(qū)位于第一級青藏高原和第二級長江中下游平原的過渡帶，海拔高差懸殊，地貌復(fù)雜多樣，以山地地貌為主，西高東低，全省可分為東部四川盆地、川西北高原和川西南山地等三大部分。

區(qū)內(nèi)季風(fēng)氣候明顯，氣候類型多樣，差異顯著。東部四川盆地為亞熱帶濕潤、半濕潤氣候，雨熱同季，水熱條件好，年均溫16～18℃，年降水量達(dá)1000～1200mm，50%以上集中在夏季；西北部為高山高原高寒氣候，氣候垂直變化明顯，總體上以寒溫帶氣候?yàn)橹?，冬寒夏涼，水熱不足，年均氣?～12℃，年降水量500～900mm；

1.2 數(shù)據(jù)來源

四川省144個氣象站點(diǎn)的年降水量（觀測精度為0.1mm）、經(jīng)緯度等數(shù)據(jù)來自中國氣象局國家氣象信息中心提供的“中國基本、基準(zhǔn)和一般地面氣象觀測站1981?2010年累年值年值數(shù)據(jù)集”。各氣象站點(diǎn)的空間分布如圖1所示。該數(shù)據(jù)集均為地面實(shí)測數(shù)據(jù)，可信度高，且站點(diǎn)數(shù)量多，能夠保證空間插值的精度，滿足研究要求。

研究區(qū)域的地形數(shù)據(jù)（Digital Elevation Model，DEM）來源于目前應(yīng)用最為廣泛的美國太空總署（NASA）空間科學(xué)數(shù)據(jù)中心提供的SRTM數(shù)據(jù)，空間分辨率為90m，研究前對DEM數(shù)據(jù)進(jìn)行了填洼處理，移除其中的小缺陷，并采用Bilinear法重采樣至1km。其它衍生地形數(shù)據(jù)（坡度、坡向等）通過1km分辨率的DEM數(shù)據(jù)在ArcGIS軟件中利用空間分析工具箱的Slope和Aspect工具提取得到。其中坡向指地表一點(diǎn)的切平面的法線矢量n在水平面的投影nxoy與過該點(diǎn)的正北方向的夾角。

圖1 四川省144個氣象站點(diǎn)分布圖

1.3 研究方法

區(qū)域降水的空間分布受諸多因素的影響，而本研究區(qū)地處內(nèi)陸山區(qū)，地形特征是影響其降水量空間分布的主要因素[14?15]。因此，選擇表達(dá)空間位置的經(jīng)度、緯度，以及表達(dá)地形的海拔、坡度和坡向共5個影響因子作為輔助變量，用于回歸分析和空間插值。

1.3.1 普通克里格插值法（OK）

克里格插值法是一種基于地統(tǒng)計(jì)學(xué)的經(jīng)典空間插值方法，目前已在多個領(lǐng)域得到廣泛應(yīng)用。但對于部分受多種因素影響、具有高度空間異質(zhì)性的目標(biāo)變量而言，傳統(tǒng)的普通克里格法（Ordinary Kriging，OK）的插值精度往往受到較大的限制[16]。

1.3.2 回歸克里格插值法（RK）

回歸克里格（RK）法[17]是將克里格插值法與回歸分析結(jié)合起來的一種空間插值方法，已有研究表明，與OK法相比，RK法通常能夠取得更高的預(yù)測精度[18?20]。RK法的基本思想是采用不同方法對目標(biāo)變量的主要影響因素和隨機(jī)因素（即回歸殘差）分別進(jìn)行分析，首先建立目標(biāo)變量和影響因素之間的回歸方程，分離出趨勢項(xiàng)，然后對回歸殘差采用OK法進(jìn)行插值，最后對回歸分析得到的趨勢項(xiàng)和OK法得到的殘差的估計(jì)值進(jìn)行求和，最終得到目標(biāo)變量的預(yù)測值。RK法可表達(dá)為[12]

式中，y(pi)為位置pi處的目標(biāo)變量，βj,(pi)為回歸系數(shù)，xj,(pi)為位置pi處的第j個自變量，ε(pi)為位置pi處的回歸殘差，n為自變量的個數(shù)。

（1）全局回歸（GR）方法。在RK法的應(yīng)用中，回歸分析部分大多采用基于最小二乘法（Ordinary Least Square，OLS）的全局回歸（Global Regression，GR），即全局回歸克里格法（GRK）。由于其僅對解釋變量進(jìn)行“均值”估計(jì)，各回歸系數(shù)為常數(shù)，不能反映變量的空間非平穩(wěn)性，因此，在處理土壤、降水等這類空間非平穩(wěn)性較強(qiáng)的變量時，GRK法插值精度的提高仍受到限制。

（2）地理加權(quán)回歸（GWR）方法。針對上述問題，F(xiàn)ortheringham等基于局部光滑思想，提出了地理加權(quán)回歸（Geographically Weighted Regression，GWR）模型用以處理回歸分析中解釋變量的空間非平穩(wěn)性[21?22]。GWR模型是對普通GR模型的擴(kuò)展，其采用了局部回歸（Local Regression，LR）方法，對空間非平穩(wěn)性進(jìn)行量化。在GWR中，解釋變量的回歸系數(shù)不再是常數(shù)，而是空間位置的函數(shù)（稱之為空間權(quán)重函數(shù)），式（1）中位置(ui,vi)處的回歸系數(shù)β通過下式計(jì)算得到[22]

式中，β為回歸系數(shù)矩陣；X為解釋變量設(shè)計(jì)矩陣，XT為其轉(zhuǎn)置矩陣；Y為因變量矩陣；W(ui,vi)為空間權(quán)重矩陣，由空間權(quán)重函數(shù)W(i)求得，其作用是定量衡量領(lǐng)域內(nèi)不同空間位置j（j=1，2，…，n）樣點(diǎn)的觀測值對于回歸點(diǎn)(ui, vi)回歸系數(shù)估計(jì)的影響程度。GWR模型能有效捕捉空間非平穩(wěn)性特征，更加真實(shí)地反映出目標(biāo)變量的空間變異情況。而基于GWR構(gòu)建的地理加權(quán)回歸克里格（Geographically Weighted Regression Kriging，GWRK）法能使空間插值的預(yù)測精度得到進(jìn)一步提高[11?12]。

（3）混合地理加權(quán)方法（MGWR）。在許多實(shí)際問題中，解釋變量中往往同時包含了全局變量和局部變量：一些解釋變量具有明顯的空間非平穩(wěn)性，適合使用GWR法建模；而另一部分解釋變量可能不具有空間非平穩(wěn)性，或其空間非平穩(wěn)性非常小可忽略不計(jì)，宜用多元GR法建模。為此，Brunsdon等在GWR模型基礎(chǔ)上，提出了混合地理加權(quán)模型（Mixed Geographical Weighted Regression，MGWR）[23]。在MGWR模型中，同時包含了回歸系數(shù)恒定的全局變量和隨地理位置變化而變化的局部變量，MGWR模型的基本形式為

式中，yi、βk,(ui,vi)和εi的含義同式（2），xk,i為位置i處的第k個局部變量，xl,i為位置i處的第l個全局變量，γl為全局回歸的回歸系數(shù)（常數(shù)），m為局部變量的個數(shù)，n為解釋變量的個數(shù)。

運(yùn)用MGWR模型的關(guān)鍵在于通過判定各解釋變量在進(jìn)行回歸分析時的空間非平穩(wěn)性類型，即某個解釋變量屬于全局變量還是局部變量。但現(xiàn)有研究大多通過主觀經(jīng)驗(yàn)、定性判斷來確定，基于定量方法分析該問題的報(bào)道還較少見。

MGWR模型可解決多個解釋變量不同空間平穩(wěn)特性共存的問題，較之單一的GR模型與GWR模型，能更準(zhǔn)確地反映變量間的空間變化關(guān)系，有效提高模型的解釋能力?降水的影響因素中可能同時包含空間平穩(wěn)和空間非平穩(wěn)兩種類型，MGWR法適用研究降水與其影響因素間的相互關(guān)系。

本研究在深入分析各影響因素的空間平穩(wěn)性特征的基礎(chǔ)上，采用混合地理加權(quán)回歸克里格法，即將回歸克里格法中回歸分析部分替換為MGWR模型，其余步驟、方法保持不變，對四川省年降水?dāng)?shù)據(jù)進(jìn)行空間插值，并探討具有高預(yù)測精度、合適本區(qū)域復(fù)雜地形條件的降水?dāng)?shù)據(jù)空間插值方法。

1.4 模型精度的評價

采用五折交叉驗(yàn)證法（Five-fold cross validation）對各模型的插值精度進(jìn)行評價。將144個氣象站點(diǎn)隨機(jī)分為5等份（每份29個站點(diǎn)），每次拿出其中一個（站點(diǎn)總數(shù)的20%）作為驗(yàn)證集，驗(yàn)證數(shù)據(jù)不參與回歸分析和克里格插值過程，只用于模型精度的評價。將其余4份數(shù)據(jù)（站點(diǎn)總數(shù)的80%）作為建模集進(jìn)行空間插值，該過程共進(jìn)行5次，并求取各驗(yàn)證數(shù)據(jù)集站點(diǎn)的年降水量模擬值。通過計(jì)算平均誤差（Mean Error，ME）、平均絕對誤差（Mean Absolute Error，MAE）、均方根誤差（Root Mean Square Error，RMSE）[12]、平均絕對百分誤差（Mean Absolute Relative Error，MARE）和均方根相對誤差（Root Mean Square Relative Error，RMSRE）等指標(biāo)的5次平均值評價不同插值模型的預(yù)測精度。其中，MARE和RMSRE定義為

式中，ya,i為位置i處的年降水量實(shí)測值，ye,i為位置i處的年降水量模擬值，n為驗(yàn)證數(shù)據(jù)集的站點(diǎn)個數(shù)。

1.5 分析步驟

（1）對各氣象站點(diǎn)的降水?dāng)?shù)據(jù)進(jìn)行預(yù)處理，包括描述性統(tǒng)計(jì)分析、正態(tài)性檢驗(yàn)等。

（2）分析研究區(qū)域站點(diǎn)的年降水量（目標(biāo)變量）與各影響因子變量間的相關(guān)性，采用逐步回歸分析方法，篩選出用于后續(xù)回歸分析的影響因子變量。

（3）采用GR、GWR和MGWR三種方法，以站點(diǎn)年降水量為因變量，上一步篩選的影響因子為解釋變量進(jìn)行回歸分析。

（4）用OK法對上述各回歸分析的殘差進(jìn)行空間插值，并將殘差插值結(jié)果與回歸分析得到的趨勢項(xiàng)相加，得到四川省年降水量的空間分布圖。

（5）通過驗(yàn)證數(shù)據(jù)集評價不同插值模型的預(yù)測精度。

1.6 軟件平臺

描述性統(tǒng)計(jì)分析采用SPSS 23完成，基于OLS的GR、GWR和MGWR等回歸分析通過GWR4完成，半方差函數(shù)計(jì)算通過地統(tǒng)計(jì)學(xué)軟件GS+9進(jìn)行，OK法插值和降水量空間分布圖的繪制由ArcGIS 10.4完成。各變量空間分布圖的分辨率均為1km。

2 結(jié)果與分析

2.1 建模前準(zhǔn)備

2.1.1 建模用降水?dāng)?shù)據(jù)序列的正態(tài)性檢驗(yàn)

對建模用氣象站點(diǎn)的年降水量數(shù)據(jù)進(jìn)行描述性統(tǒng)計(jì)分析和正態(tài)性檢驗(yàn)，表1表明，研究區(qū)域各站點(diǎn)的年降水量介于346.9～1670.0mm，平均968.1mm；變異系數(shù)達(dá)到23.06%，變異程度中等。站點(diǎn)的年降水量數(shù)據(jù)序列通過了J?B（Jarque-Bera）檢驗(yàn)（P＜0.05），符合正態(tài)分布。

表1 研究區(qū)站點(diǎn)的年降水量描述性統(tǒng)計(jì)分析和正態(tài)性檢驗(yàn)（1981?2010年）

2.1.2 模型影響因子變量篩選

年降水量與各影響因子的Pearson相關(guān)系數(shù)見表2。由表可知，除坡向外，年降水量與其余各影響因子間均存在較高的相關(guān)性，具體表現(xiàn)為與經(jīng)度和海拔呈顯著正相關(guān)，與坡度呈極顯著負(fù)相關(guān)，與緯度呈顯著正相關(guān)。另一方面，各影響因子間也存在一定程度的相互作用，如海拔與經(jīng)度呈極顯著負(fù)相關(guān)，相關(guān)系數(shù)達(dá)到?0.788，可見，部分影響因子間存在較強(qiáng)多重共線性（Multicollinearity），將對插值模型的構(gòu)建產(chǎn)生不利影響。

表2 年降水量與各影響因子的相關(guān)系數(shù)（115個建模站點(diǎn)）

注：*、**分別表示相關(guān)系數(shù)通過0.05、0.01水平的顯著性檢驗(yàn)。下同。

Note:*is P<0.05,**is P<0.01. The same as below.

表3 年降水量與各影響因子的逐步回歸過程（n=115）

表4 表3中模型5的回歸計(jì)算結(jié)果（n=115）

2.2 地理加權(quán)回歸模型的構(gòu)建

以上一步篩選的經(jīng)度、緯度和坡度3個影響因子為局部變量進(jìn)行地理加權(quán)回歸，結(jié)果見表5。由表中可見，各解釋變量的回歸系數(shù)處于一定的變化范圍，不再是常數(shù)，如經(jīng)度的回歸系數(shù)的取值范圍介于?118.89～265.63，變異系數(shù)達(dá)到283.21%，說明降水量與各解釋變量間的相互關(guān)系具有顯著的空間非平穩(wěn)性。

通過ΔAICc判定各解釋變量的空間非平穩(wěn)性類型。AICc是用于檢驗(yàn)回歸模型及模型中變量擬合優(yōu)度的常用指標(biāo)之一，其值越小，擬合效果越佳[25]。本研究中ΔAICc是指同一變量在GWR和GR兩種回歸模型的AICc值之差，其值為正，即GWR模型的AICc值大于GR模型的AICc值，說明將該變量視為局部變量、采用GWR模型進(jìn)行回歸的效果不如GR模型。Nakaya等認(rèn)為，若ΔAICc值大于2，可以認(rèn)為該變量不存在明顯的空間非平穩(wěn)性，進(jìn)行回歸分析時應(yīng)將其從局部變量調(diào)整為全局變量[26]。從3個解釋變量的ΔAICc值來看，經(jīng)度和緯度的ΔAICc值均為負(fù)數(shù)，說明在進(jìn)行回歸分析時將這些變量作局部變量處理是合理的；而坡度的ΔAICc值達(dá)11.359，遠(yuǎn)大于2，說明其空間非平穩(wěn)性不明顯，應(yīng)將該變量作為全局變量處理。因此，以站點(diǎn)的經(jīng)度和緯度為局部變量，坡度為全局變量，采用MGWR法再次進(jìn)行回歸分析。由表6的MGWR分析結(jié)果可見，將坡度設(shè)為全局變量后，其回歸系數(shù)也隨之變?yōu)槌?shù)（?3.425），說明總體上坡度與年降水量之間為負(fù)相關(guān)關(guān)系。

表5 地理加權(quán)回歸模型（GWR）回歸計(jì)算結(jié)果（n=115）

注：ΔAICc是同一變量在GWR和GR兩種回歸模型的AICc值之差。

Note: ΔAICc is the difference between the value of AICc of the same variable calculated by GWR model and that by GR model.

表6 混合地理加權(quán)模型（MGWR）回歸計(jì)算結(jié)果（n=115）

表7 三種回歸模型的診斷指標(biāo)

2.3 回歸模型殘差的克里格插值

正態(tài)性檢驗(yàn)結(jié)果顯示，GR、GWR和MGWR三種模型的殘差不符合正態(tài)分布。因此，對殘差數(shù)據(jù)進(jìn)行3/5次方轉(zhuǎn)換，經(jīng)轉(zhuǎn)換后使其滿足進(jìn)行克里格插值的條件。在GS+9軟件中，對3種經(jīng)3/5次方轉(zhuǎn)換后的回歸殘差進(jìn)行半方差函數(shù)分析。從表8可看出，各半方差函數(shù)的決定系數(shù)R2均達(dá)到0.60以上，取得了較好的擬合效果；各半方差函數(shù)的塊基比［即C0/ （C0+ C）］均小于20%，說明各回歸殘差存在較強(qiáng)的空間自相關(guān)性，適合使用克里格法進(jìn)行空間插值。

表8 各回歸殘差的半方差函數(shù)擬合結(jié)果

根據(jù)半方差函數(shù)的分析結(jié)果，采用OK法分別對各模型的殘差進(jìn)行空間插值，得到回歸殘差的空間分布圖（回歸殘差進(jìn)行3/5次方轉(zhuǎn)換后進(jìn)行插值，得到插值結(jié)果后再進(jìn)行逆轉(zhuǎn)換）。通過各影響因子變量及其回歸系數(shù)的空間分布圖、殘差插值結(jié)果圖等分別獲得基于GRK、GWRK和MGWRK模型的四川省年降水量的空間分布圖。同時采用OK法直接對各站點(diǎn)的年降水?dāng)?shù)據(jù)進(jìn)行空間插值。各模型的插值結(jié)果如圖2所示。

由圖2可見，各插值模型得到的四川省年降水量空間分布的總趨勢基本一致，均呈現(xiàn)東高西低、盆地高高原低、盆周山地高盆中丘陵低的趨勢，該預(yù)測結(jié)果與前人的研究基本一致[27]。此外，一些降水空間分布的細(xì)節(jié)在插值結(jié)果中也有所反映。如全省年降水的峰值區(qū)大致位于四川盆地西部邊緣的山前丘陵、中山區(qū)，包括雅安的名山、寶興和樂山的峨眉山等地。該區(qū)域呈長軸方向，為西北—東南向的橢圓形，區(qū)內(nèi)各地的年降水量達(dá)到1300～1500mm。形成這一現(xiàn)象的主要原因是該區(qū)域西部、西北部和南部均有一系列高山阻隔，形成了特殊的“喇叭”狀的地形，使東來的太平洋東南暖濕氣流與盆周山地下沉的冷濕氣流交匯于此，形成了著名的“華西雨屏”現(xiàn)象，使該區(qū)成為全省乃至中國內(nèi)陸降水量最大的地區(qū)。

圖2 基于不同模型的四川省年平均降水量空間分布（1981?2010年）

2.4 地理回歸克里格插值模型的精度評價

采用五折交叉驗(yàn)證法，通過ME、MAE和MARE等指標(biāo)評價各插值模型的預(yù)測精度。各評價指標(biāo)的值越接近0，說明模型的插值精度越高。表9表明，各模型的插值精度大致表現(xiàn)為MGWRK＞GWRK ≈GRK＞OK。綜合分析表明，包含空間位置、地形等多個影響因素的3種回歸克里格模型對年降水量的插值精度較OK模型均有較大提高，如RMSRE值從OK的1.351顯著減至MGWRK的0.897；GWRK模型雖然考慮了影響因子空間非平穩(wěn)性，在回歸效果上也優(yōu)于GRK模型，但最終的插值精度與后者相比并沒有明顯的提高；而在4種模型中，進(jìn)一步探究了各影響因子的空間非平穩(wěn)性類型的MGWRK模型的插值精度最優(yōu)。

根據(jù)144個站點(diǎn)的空間分布特征，通過Thiessen多邊形法對研究區(qū)進(jìn)行剖分，得到MGWRK模型插值殘差（即實(shí)測年降水?dāng)?shù)據(jù)與MGWRK模型模擬值之差）的空間分布（圖3）。從圖3可看出，MGWRK模型的殘差值介于?240～258，且分布較為集中，殘差值介于?50～50的站點(diǎn)共有112個，占站點(diǎn)總數(shù)的77.78%，占全省總面積的81.91%。由此可見，采用MGWRK模型獲得的年降水量數(shù)據(jù)在四川省大多數(shù)地區(qū)可取得較好的計(jì)算精度，插值殘差總體而言處于一個較低水平。而殘差高值區(qū)主要分布在研究區(qū)中西部的四川盆地西部邊緣的山前丘陵、中山區(qū)，即上述的“華西雨屏區(qū)”，包括天全、瀘定和峨眉山等站點(diǎn)。

表9 不同模型的插值精度分析

圖3 MGWRK模型插值殘差的空間分布

3 結(jié)論與討論

3.1 討論

區(qū)域降水的空間分布受到多方面因素的影響。前人研究表明，空間位置、大氣環(huán)流以及地形因素均會對區(qū)域降水有一定的制約作用[14?15]。本研究在進(jìn)行分析時充分考慮區(qū)域降水形成的物理過程，將空間位置和地形等主要影響因子納入插值過程。研究結(jié)果表明，四川省年降水的整體空間分布主要受空間位置因素（即經(jīng)度和緯度）的影響，GR模型的計(jì)算結(jié)果也證明了這一點(diǎn)：經(jīng)度（55.804）和緯度（?47.083）的回歸系數(shù)均遠(yuǎn)大于代表地形因子的坡度的回歸系數(shù)（?7.915），說明整體上看，空間位置是影響四川省降水空間分布的主控因素，在較大空間尺度下，其對降水的影響占據(jù)主要地位，并導(dǎo)致地形對降水分布的影響相對不明顯。另一方面，地形因素對降水的影響主要表現(xiàn)在局部地區(qū)，并往往會出現(xiàn)異常的地形降水分布。如四川盆地西部邊緣的山前丘陵、中山區(qū)特殊的“喇叭”狀地形，導(dǎo)致東來的太平洋東南暖濕氣流與盆周山地下沉的冷濕氣流交匯于此，形成了著名的“華西雨屏”現(xiàn)象，使得該區(qū)的年降水量要明顯高于周圍地區(qū)。

綜合各種評價誤差的計(jì)算結(jié)果可以看出，由于深入考慮了各影響因子與降水量相互關(guān)系的空間非平穩(wěn)性，MGWRK法對提高降水空間插值的精度作用是明顯的，其插值效果優(yōu)于OK法等傳統(tǒng)插值方法，能較準(zhǔn)確地反映四川省降水分布實(shí)際狀況，尤其是一些降水分布的細(xì)節(jié)性特征能較好地表達(dá)出來；在氣象站點(diǎn)相對較少的四川省西部地區(qū)，其插值精度也能得到較好的保證。另一方面，不同地區(qū)的插值精度也存在一定差異。MGWRK法在降水分布符合主控因素（即空間位置）的地區(qū)具有較好的插值精度；而當(dāng)降水分布出現(xiàn)異常，即地形因素對降水分布的作用增強(qiáng)時，如在上述年降水量遠(yuǎn)高于周圍地區(qū)的華西雨屏帶，插值效果則有所下降，誤差相對較大，這與鄔倫等學(xué)者的研究結(jié)論也是一致的[28?29]。

從時間尺度來看，本研究結(jié)果表明，就四川省多年平均氣候狀態(tài)下的降水而言，其空間分布的變化規(guī)律較明顯，影響因素較少且明確，因此取得了較好的插值結(jié)果。但在實(shí)際應(yīng)用時，對于較小時間尺度如日降水、逐小時降水等，影響降水空間分布的因素往往會增多，且隨機(jī)性因素的作用也可能增強(qiáng)，MGWRK模型是否仍然適用，此外，對于制約降水分布的影響因素，本研究僅考慮了海拔、坡度、坡向等地形因子，尚未考慮大氣環(huán)流等其它因子對降水的影響，加入這些影響因素是否有助于插值精度的進(jìn)一步提高，這些問題都有待進(jìn)一步深入研究，對插值模型進(jìn)行改進(jìn)。

3.2 結(jié)論

（1）用于回歸分析的影響因子為經(jīng)度、緯度和坡度，在MGWR回歸分析中，坡度變量不存在明顯的空間非平穩(wěn)性，作為全局變量處理，而經(jīng)度和緯度作為局部變量處理。

（2）綜合考慮了空間位置、地形等多個影響因素的3種回歸克里格模型（GRK、GWRK和MGWRK）的插值精度均優(yōu)于OK模型，而進(jìn)一步考慮了各影響因子與降水量相互關(guān)系的空間非平穩(wěn)性，并深入探究了空間非平穩(wěn)性類型的MGWRK模型的插值效果最佳。

（3）MGWRK插值模型以氣象站點(diǎn)實(shí)測數(shù)據(jù)為基礎(chǔ)，在考慮解釋變量空間平穩(wěn)性特征的基礎(chǔ)上，進(jìn)一步定量區(qū)分了解釋變量的空間平穩(wěn)性類型，由此獲得的降水空間數(shù)據(jù)具有較高的插值精度，并能較準(zhǔn)確地表達(dá)局部區(qū)域降水發(fā)生突變的細(xì)節(jié)信息。在空間插值研究中，輔助信息的加入有助于一定程度上提高預(yù)測精度。在此基礎(chǔ)上，進(jìn)一步厘清不同輔助信息在插值過程中所扮演的“角色”，如不同輔助變量的空間平穩(wěn)性類型，可使插值精度得到進(jìn)一步提高。

[1] 葛朝霞,曹麗青.氣象學(xué)與氣候?qū)W教程[M].北京:中國水利水電出版社,2009:213-215.Ge Z X,Cao L Q.A course on meteorology and climatology[M].Beijing:China Water and Power Press,2009:213-215.(in Chinese)

[2] 劉元波,傅巧妮,宋平,等.衛(wèi)星遙感反演降水研究綜述[J].地球科學(xué)進(jìn)展,2011,26(11):1162-1172.Liu Y B,Fu Q N,Song P,et al.Satellite retrieval of precipita- tion:an overview[J].Advances in Earth Sciences,2011, 26(11):1162-1172.(in Chinese)

[3] 王智,吳友均,梁鳳超,等.新疆地區(qū)年降水量的空間插值方法研究[J].中國農(nóng)業(yè)氣象,2011,32(3):331-337.Wang Z,Wu Y J,Liang F C,et al.Study on spatial interpolation method of annual precipitation in Xinjiang[J].Chinese Journal of Agrometeorology,2011,32(3):331-337.(in Chinese)

[4] 趙煜飛,朱江.近50年中國降水格點(diǎn)日值數(shù)據(jù)集精度及評估[J].高原氣象,2015,34(1):50-58.Zhao Y F,Zhu J.Assessing quality of grid daily precipitation datasets in China in recent 50 years[J].Plateau Meteoro- logy,2015,34(1):50-58.(in Chinese)

[5] 李月,齊實(shí),程伯涵,等.哀牢山山區(qū)降水時空分布的影響因素及插值方法比較[J].地球與環(huán)境,2017,45(6):600-611.Li Y,Qi S,Cheng B H,et al.A study on factors of space-time distributions of precipitation in Ailao Mountain area and comparison of interpolation methods[J].Earth and Environ- ment,2017,45(6):600-611.(in Chinese)

[6] 王紅霞,柳小妮,李純斌,等.甘肅省近42年降水量變化時空分布格局分析[J].中國農(nóng)業(yè)氣象,2013,34(4):384-389.Wang H X,Liu X N,Li C B,et al.Spatial temporal distribution of precipitation in Gansu Province last 42 years[J].Chinese Journal of Agrometeorology,2013,34(4):384-389.(in Chinese)

[7] 龍軍,張黎明,沈金泉,等.復(fù)雜地貌類型區(qū)耕地土壤有機(jī)質(zhì)空間插值方法研究[J].土壤學(xué)報(bào),2014,51(6):81-93.Long J,Zhang L M,Shen J Q,et al.Spatial interpolation of soil organic matter in farmlands in areas complex in landform[J].Acta Pedologica Sinica,2014,51(6):81-93.(in Chinese)

[8] 趙晨曦,王云琦,王玉杰,等.北京地區(qū)冬春PM2.5和PM10污染水平時空分布及其與氣象條件的關(guān)系[J].環(huán)境科學(xué),2014,35(2):418-427.Zhao C X,Wang Y Q,Wang Y J,et al.Temporal and spatial distribution of PM2.5 and PM10 pollution status and the correlation of particulate matters and meteorological factors during winter and spring in Beijing[J].Environmental Science,2014,35(2):418-427.(in Chinese)

[9] 謝云峰,張樹文.基于數(shù)字高程模型的復(fù)雜地形下的黑龍江平均氣溫空間插值[J].中國農(nóng)業(yè)氣象,2007,28(2):205-211.Xie Y F,Zhang S W.Spatial interpolation of mean temperature of Heilongjiang Province based on digital elevation model[J].Chinese Journal of Agrometeorology,2007,28(2):205-211.(in Chinese)

[10] 耿廣坡,高鵬,呂圣橋,等.魯中南山區(qū)馬蹄峪小流域土壤有機(jī)質(zhì)和全氮空間分布特征[J].中國水土保持科學(xué),2011,9(6):99-105.Geng G P,Gao P,Lv S Q,et al.Spatial distribution of soil organic matter and total nitrogen in Matiyu small watershed in hilly area of middle southern Shandong Province[J].Science of Soil and Water Conservation,2011,9(6):99-105.(in Chinese)

[11] 張國峰,楊立榮,瞿明凱.基于地理加權(quán)回歸克里格的日平均氣溫插值[J].應(yīng)用生態(tài)學(xué)報(bào),2015,26(5):1531-1536.Zhang G F,Yang L R,Qu M K.Interpolation of daily mean temperature by using geographically weighted regression-Kriging[J].Chinese Journal of Applied Ecology,2015,26(5):1531-1536.(in Chinese)

[12] 楊順華,張海濤,郭龍,等.基于回歸和地理加權(quán)回歸Kriging的土壤有機(jī)質(zhì)空間插值[J].應(yīng)用生態(tài)學(xué)報(bào),2015,26(6):1649-1656.Yang S H,Zhang H T,Guo L,et al.Spatial interpolation of soil organic matter using regression Kriging and geographically weighted regression Kriging[J].Chinese Journal of Applied Ecology,2015,26(6):1649-1656.(in Chinese)

[13] 四川省地方志工作辦公室.四川年鑒2016[M].成都:四川年鑒社,2016:21-22.Office of Local Chronicles Compilation of Sichuan Province.Sichuan yearbook 2016[M].Chengdu:Sichuan Yearbook Press,2016:21-22.(in Chinese)

[14] Spreen W C.A determination of the effect of topography upon precipitation[J].Eos Transactions American Geophy- sical Union,1947,28:285-290.

[15] Smith R B.The influence of mountains on the atmosphere[J].Advances in Geophysics,1979,21(4):87-230.

[16] Odeh I O A,Mcbratney A B,Chittleborough D J.Spatial prediction of soil properties from landform attributes derived from a digital elevation model[J].Geoderma,1994,63:197-214.

[17] Juang K W,Lee D Y.A comparison of three Kriging methods using auxiliary variables in heavy-metal contaminated soils[J].Journal of Environmental Quality,1998,27(2):355-363.

[18] 張慧智,史學(xué)正,于東升,等.中國土壤溫度的空間插值方法比較[J].地理研究,2008,27(6):1299-1307.Zhang H Z,Shi X Z,Yu D S,et al.Spatial prediction of soil temperatures in China using different methods[J].Geograp- hical Research,2008,27(6):1299-1307.(in Chinese)

[19] 王庫.回歸克里格在土壤全氮空間預(yù)測上的應(yīng)用[J].中國農(nóng)學(xué)通報(bào),2013,29(20):142-147.Wang K.Application of regression Kriging on the spatial prediction of total soil nitrogen[J].Chinese Agricultural Science Bulletin,2013,29(20):142-147.(in Chinese)

[20] 曹祥會,龍懷玉,周腳根,等.河北省表層土壤有機(jī)碳和全氮空間變異特征性及影響因子分析[J].植物營養(yǎng)與肥料學(xué)報(bào),2016,22(4):937-948.Cao X H,Long H Y,Zhou J G,et al.Analysis of spatial variability and influencing factors of topsoil organic carbon and total nitrogen in Hebei Province[J].Plant Nutrition and Fertilizer Science,2016,22(4):937-948.(in Chinese)

[21] Fotheringham A S,Charlton M E,Brunsdon C.The geography of parameter space:an investigation into spatial non-station- arity[J].International Journal of Geographical Information Systems,1996,10(5):605-627.

[22] Fotheringham A S,Brunsdon C,Charlton M E.Geographi- cally weighted regression:the analysis of spatially varying relati- onships[M].Chichester:John Wiley & Sons Ltd.,2002:16-24.

[23] Brunsdon C,Fotheringham A S,Charlton M E.Some notes on parametric significance tests for geographically weighted regression[J].Journal of Regional Science,1999,39(3):497-524.

[24] Akaike H.A new look at the statistical model identification[J].IEEE Transactions on Automatic Control,1974,19(6):716-723.

[25] McQuarrie A D R,Tsai C L.Regression and time series model selection[J].World Scientific,1998,42(2):480.

[26] Nakaya T,Fotheringham A S,Brunsdon C,et al.Geographic- ally weighted Poisson regression for disease association mapping[J].Statistics in Medicine,2005,24(17):2695-2717.

[27] 周長艷,岑思弦,李躍清,等.四川省近50年降水的變化特征及影響[J].地理學(xué)報(bào),2011,66(5):619-630.Zhou C Y,Cen S X,Li Y Q,et al.Precipitation variation and its impacts in Sichuan in the last 50 years[J].Acta Geograph- ica Sinica,2011,66(5):619-630.(in Chinese)

[28] 鄔倫,吳小娟,肖晨超,等.五種常用降水量插值方法誤差時空分布特征研究:以深圳市為例[J].地理與地理信息科學(xué),2010,26(3):19-24.Wu L,Wu X J,Xiao C C,et al.On temporal and spatial error distributions of five precipitation interpolation models: a case of Shenzhen[J].Geography and Geo-Information Scien- ce,2010,26(3):19-24.(in Chinese)

[29] 熊秋芬,黃玫,熊敏詮,等.基于國家氣象觀測站逐日降水格點(diǎn)數(shù)據(jù)的交叉檢驗(yàn)誤差分析[J].高原氣象,2011,30(6):1615-1625.Xiong Q F,Huang M,Xiong M Q,et al.Cross-validation error analysis of daily gridded precipitation based on China meteorological observation[J].Plateau Meteorology,2011,30(6):1615-1625.(in Chinese)

Spatial Interpolation of Regional Precipitation Based on Mixed Geographical Weighted Regression Combined with Kriging Interpolation

LI Hao, LIU Tao, XU Jing-wen

(College of Resources Science and Technology, Sichuan Agricultural University, Chengdu 611130, China)

Based on the precipitation data of 1981?2010 from 144 meteorological stations in Sichuan province, using mixed geographical weighted regression Kriging interpolation (MGWRK) model, and considering the impact of topographic factors, the spatial distribution of the average annual precipitation was obtained in this paper. The effect of interpolation value was compared with those values from OK, GRK, and GWRK methods. The result showed that the optimal influencing factors combination was longitude, latitude and slope, determined by using the stepwise regression method, could decrease the multi-collinearity among the explanatory variables significantly. The types of spatial variability of the explanatory variables were analyzed quantitatively based on the index ΔAICc, which was the difference between the value of AICc (Corrected Akaike Information Criterion) of the same variable calculated by GWR model and by GR model. Then set the slope variable as global variable, and the longitude and latitude variables as local variables, the interpolation of the average annual precipitation in Sichuan province was conducted by the MGWRK model. The MGWRK method presented in this paper showed higher accuracy than those of the ordinary Kriging (OK) and global regression Kriging (GRK), because the method has taken into consideration of various influence factors of the spatial position and topography, and the variability of the relationship between these factors and precipitation.

Precipitation; Mix geographically weighted regression; Kriging interpolation; Spatial Interpolation

10.3969/j.issn.1000-6362.2018.10.006

李豪,劉濤,徐精文.基于混合地理加權(quán)回歸與克里格的區(qū)域降水量空間插值方法[J].中國農(nóng)業(yè)氣象,2018,39(10):674?684

2018?02?08

國家自然科學(xué)基金項(xiàng)目（41501291）；四川省教育廳科研基金項(xiàng)目（14ZB0009）

李豪（1980?），博士，講師，主要從事3S技術(shù)在水土資源可持續(xù)利用方面研究。E-mail: lihao@sicau.edu.cn