王雪珊, 沈慶松,, 高鳳杰, 張興義, 張少良, 王 力

(1.東北農(nóng)業(yè)大學(xué) 資源與環(huán)境學(xué)院, 哈爾濱 150030; 2.中國科學(xué)院 東北地理與農(nóng)業(yè)生態(tài)研究所,哈爾濱 150081; 3.西北農(nóng)林科技大學(xué) 水土保持研究所, 陜西 楊凌 712100)

土壤速效磷(Available phosphorus,AP)是衡量土壤供磷能力和磷素流失風(fēng)險(xiǎn)的重要指標(biāo)[1-2]，該指標(biāo)能夠反映農(nóng)田土壤肥力特征，研究其空間分布規(guī)律對(duì)于區(qū)域施肥管理和農(nóng)業(yè)面源污染防控等具有重要意義[3]。受生物地球化學(xué)循環(huán)過程、自然環(huán)境和人類活動(dòng)等因子影響[4-5]，農(nóng)田土壤AP在空間分布上存在高度空間異質(zhì)性。研究證實(shí)AP預(yù)測(cè)的精度直接影響農(nóng)田土壤磷素的管理水平[6]，因此利用有限信息，改進(jìn)模擬方法提高AP空間分布預(yù)測(cè)模型的精度十分迫切。

空間插值方法是利用采樣點(diǎn)數(shù)據(jù)通過數(shù)學(xué)模型對(duì)研究區(qū)內(nèi)未知點(diǎn)進(jìn)行模擬預(yù)測(cè)，并能夠獲得連續(xù)的屬性值空間分布特征[7]?；貧w克里格法[8]將多元線性回歸模型與克里格方法相結(jié)合，同時(shí)引入多個(gè)環(huán)境因子參與建模，并考慮土壤屬性空間分布過程中的結(jié)構(gòu)性因子與隨機(jī)性因子，其空間插值結(jié)果更加精確[9-10]。Li[10]和Zhang等[11]提取地形和土地利用等變量利用回歸克里格法預(yù)測(cè)土壤有機(jī)質(zhì)的空間分布特征，研究結(jié)果表明回歸克里格方法要優(yōu)于單變量空間插值方法。但是空間非平穩(wěn)性的存在導(dǎo)致回歸克里格方法不能有效地捕捉土壤屬性空間變異的局部特征[12-13]。地理加權(quán)回歸模型利用局部加權(quán)最小二乘法進(jìn)行逐點(diǎn)參數(shù)估計(jì)和空間預(yù)測(cè)，能夠在引入環(huán)境變量的同時(shí)基于空間非平穩(wěn)性假設(shè)揭示更多的局部信息，該模型已逐漸成為土壤屬性空間預(yù)測(cè)的有效手段[14-15]。地理加權(quán)回歸克里格將地理加權(quán)回歸模型與克里格方法相結(jié)合，既能夠考慮空間信息的局部特征，還能利用克里格法對(duì)代表隨機(jī)性的殘差進(jìn)行插值，插值結(jié)果能夠揭示可能被空間非平穩(wěn)性所掩蓋的一些局部變化，反映出更加真實(shí)的土壤屬性空間變異情況[16-18]。楊順華等[16]引入相對(duì)高程和徑流強(qiáng)度指數(shù)預(yù)測(cè)鄉(xiāng)鎮(zhèn)尺度(235 km2)土壤有機(jī)質(zhì)空間分布規(guī)律；馬泉來等[17]研究表明海拔和水系距離可以作為輔助變量預(yù)測(cè)黑土區(qū)小流域(119.16 km2)土壤有機(jī)質(zhì)的空間分布規(guī)律，研究證實(shí)相對(duì)于普通克里格和回歸克里格法，地理加權(quán)回歸克里格能夠顯著提高土壤有機(jī)質(zhì)的模擬精度；Kumar等[12]引入海拔、坡度、植被指數(shù)、土壤類型、土地利用類型等變量研究賓夕法尼亞州(117 599 km2)土壤有機(jī)碳的空間分布規(guī)律；Ye等[19]研究表明海拔、坡度和地形濕度指數(shù)能夠模擬區(qū)域尺度下(16 400 km2)不同采樣點(diǎn)密度土壤有機(jī)碳的空間分布格局，以上研究證明考慮了空間非平穩(wěn)性和殘差空間自相關(guān)的地理加權(quán)回歸克里格方法能夠加強(qiáng)土壤有機(jī)碳的模擬精度。然而，應(yīng)用地理加權(quán)回歸克里格時(shí)，不同屬性值在空間插值過程選用的輔助變量不同，即使同一屬性值不同土壤類型或區(qū)域所選擇的輔助變量也不完全相同，因此有必要在不同區(qū)域和不同土壤類型開展空間插值方法的研究。

東北黑土區(qū)是我國重要商品糧生產(chǎn)基地，其糧食產(chǎn)量關(guān)系到國家糧食安全問題，合理和精準(zhǔn)模擬土壤養(yǎng)分空間分布格局是區(qū)域土壤養(yǎng)分調(diào)控的關(guān)鍵。因此，為了提高黑土區(qū)AP空間分布的模擬精度，本研究嘗試將回歸克里格方法和地理加權(quán)回歸克里格方法應(yīng)用于AP空間模擬精度研究中，綜合考慮景觀格局、區(qū)域尺度、采樣方法等信息，選取東北黑土區(qū)兩個(gè)典型黑土小流域作為研究對(duì)象，比較不同空間插值方法(反距離權(quán)重法、徑向基函數(shù)法、普通克里格法、協(xié)克里格法、多元線性回歸模型、回歸克里格法、地理加權(quán)回歸模型和地理加權(quán)回歸克里格法)對(duì)黑土區(qū)AP的空間預(yù)測(cè)精度，確定黑土區(qū)AP的最優(yōu)插值方法，旨在為研究區(qū)合理施肥和農(nóng)業(yè)面源污染防控等提供理論依據(jù)和技術(shù)支撐。

1 研究區(qū)概況

光榮小流域位于黑龍江省海倫縣光榮村(47.21°—47.23°N,126.50°—126.51°E)，流域總面積1.86 km2，地形為典型的漫川漫崗(高程186—242 m)，土壤類型以典型黑土為主，當(dāng)前以大豆玉米輪作為主，一年一熟制，其基本理化形狀見表1。該地區(qū)以溫帶大陸性季風(fēng)氣候?yàn)橹?，冬季寒冷干燥，夏季溫?zé)岫嘤?，年均降?30 mm，主要集中在6—8月(65%)，年均氣溫1.5℃，年均有效積溫2 450℃，年均日照時(shí)數(shù)為2 600～2 800 h，凍融期長(zhǎng)達(dá)130 d[1](圖1，表1)。

海溝河小流域位于黑龍江省哈爾濱市阿城區(qū)東部，地理坐標(biāo)為126°55′45″—127°10′05″E,45°34′8″—45°40′50″N，屬于典型的“林地—旱田—水田”景觀格局，總面積約119.76 km2。流域東部為高山丘陵區(qū)，主要植被類型為原始的針闊混交林，土壤類型為森林暗棕壤；中部為典型的漫川漫崗區(qū)，主要耕作植物為玉米，土壤類型為典型黑土；西部地區(qū)位于流域下游區(qū)域，地勢(shì)相對(duì)平坦，以水稻種植為主，土壤類型為黑土發(fā)育的水稻土。海溝河小流域以溫帶大陸性季風(fēng)氣候?yàn)橹鳎昃邓?00～800 mm，主要集中在4—10月，占總降雨量的90%，年均氣溫為3.9°C(圖1，表1)[3,17]。

2 材料與方法

2.1 樣點(diǎn)采集與測(cè)定

本研究分別于2012年和2014年秋季在光榮和海溝河小流域采用隨機(jī)和網(wǎng)格設(shè)計(jì)法布點(diǎn)，兼顧不同土地利用方式，五點(diǎn)法采集0—20 cm耕層土壤樣品各120個(gè)(圖1)，同時(shí)記錄采樣點(diǎn)基本信息，包括經(jīng)緯度坐標(biāo)、高程、坡位、耕作方式、前茬植被等。土樣經(jīng)室內(nèi)風(fēng)干、研磨、過篩等處理，采用0.5 mol的NaHCO3浸提—鉬藍(lán)比色法測(cè)定AP含量[20]。

圖1 研究區(qū)地理位置及采樣點(diǎn)分布

2.2 環(huán)境變量提取

環(huán)境變量能夠直接或間接地反映地球生物化學(xué)循環(huán)過程，如土壤發(fā)生過程、地表徑流、淋溶、植被分布等進(jìn)而影響土壤養(yǎng)分空間分布特征[11-12]。因此，環(huán)境變量(地形因子和遙感數(shù)據(jù))常被用做輔助變量參與土壤養(yǎng)分空間異質(zhì)性分析以提高模擬精度。本研究共提取30個(gè)環(huán)境變量通過相關(guān)性分析和回歸分析揭示其與AP的關(guān)系，環(huán)境變量提取方法見表2[21]。

表1 研究區(qū)土壤基本屬性及樣點(diǎn)參數(shù)

表2 輔助變量及其提取方法

2.3 驗(yàn)證方法

為比較不同空間插值方法模擬精度，通過ArcGIS 10軟件子集要素模塊隨機(jī)自動(dòng)抽取訓(xùn)練樣本集和驗(yàn)證樣本集，樣本數(shù)量比為4∶1。運(yùn)用驗(yàn)證集樣點(diǎn)的實(shí)測(cè)值與預(yù)測(cè)值進(jìn)行精度評(píng)價(jià)，常用評(píng)價(jià)指標(biāo)包括平均誤差(ME)、平均相對(duì)誤差(MAE)、均方根誤差(RMSE)和相對(duì)提高度(RI)，相應(yīng)的計(jì)算公式如下：

(1)

(2)

(3)

RI=(RMSER-RMSEE)/RMSER×100%

(4)

式中:n為驗(yàn)證點(diǎn)數(shù)量;z(xi,yi)和z*(xi,yi)分別為第i個(gè)驗(yàn)證點(diǎn)的實(shí)測(cè)值與預(yù)測(cè)值。ME是插值無偏性的量度,其結(jié)果越接近0說明結(jié)果越無偏;MAE和RMSE是插值精度的量度,其值越小說明插值方法越精確;RI用來兩兩對(duì)比不同插值方法精度的提高程度。

3 結(jié)果與分析

3.1 土壤速效磷描述性分析

對(duì)隨機(jī)抽取的訓(xùn)練集樣本AP數(shù)據(jù)進(jìn)行經(jīng)典統(tǒng)計(jì)學(xué)分析，結(jié)果顯示光榮小流域AP含量的范圍為2.56～50.01 mg/kg，平均值為14.51 mg/kg，偏度和峰度分別為1.86，3.85，結(jié)合Kolmogorov-Smirnov test驗(yàn)證法進(jìn)行非參數(shù)檢驗(yàn)可知研究區(qū)AP數(shù)據(jù)不符合正態(tài)分布，需進(jìn)行數(shù)據(jù)轉(zhuǎn)化。變異系數(shù)為51%，表明研究區(qū)AP的空間變異程度為中等(表3)。

海溝河小流域AP含量的范圍為6.44～162.56 mg/kg，平均值為49.48 mg/kg，較光榮小流域AP含量高71%，偏度和峰度1.49，3.12，未通過K-S檢驗(yàn)，需進(jìn)行數(shù)據(jù)轉(zhuǎn)化。AP變異系數(shù)為57%，屬中等程度空間變異(表3)。

表3 研究區(qū)AP描述性統(tǒng)計(jì)特征值

3.2 土壤速效磷與輔助變量的相關(guān)性分析

光榮小流域AP與高程呈顯著正相關(guān)關(guān)系(r=0.247)；與坡度、地形起伏度、地表粗糙度、地表切割深度呈極顯著負(fù)相關(guān)關(guān)系；與交通距離和居民點(diǎn)距離呈極顯著負(fù)相關(guān)關(guān)系。由于相關(guān)地形因子之間存在多重共線性，為減少數(shù)據(jù)冗余，提高精度，本研究對(duì)與AP顯著相關(guān)的環(huán)境變量進(jìn)行主成分分析，得到第一主成分(地形因子)和第二主成分(人類活動(dòng)因子)，均與AP呈極顯著負(fù)相關(guān)，可以作為輔助變量參與建模(表4)。海溝河流域AP與高程呈顯著負(fù)相關(guān)關(guān)系(r=-0.257)；與含鐵礦物指數(shù)呈顯著負(fù)相關(guān)關(guān)系(r=0.221)；與其他環(huán)境變量并未達(dá)到顯著相關(guān)(表4)。

表4 研究區(qū)AP與輔助變量相關(guān)性分析

3.3 土壤速效磷與輔助變量的回歸分析

多元線性逐步回歸分析能夠保證與AP顯著相關(guān)的環(huán)境變量進(jìn)入回歸模型的同時(shí)去除各環(huán)境變量之間的共線性。結(jié)果表明，坡度和交通距離是光榮小流域AP空間分布規(guī)律的最佳解釋變量，引入主成分分析得到的地形因子和人類活動(dòng)因子參與回歸分析可提高回歸方程模擬精度。高程和含鐵礦物指數(shù)是海溝小流域AP空間分布規(guī)律的最佳解釋變量，各變量之間不存在多重共線性問題(表5)。

表5 研究區(qū)AP多元線性回歸方程及其相關(guān)參數(shù)

為了便于對(duì)不同插值方法模擬精度進(jìn)行比較，在進(jìn)行地理加權(quán)回歸模型建模過程中同樣選取上述因子作為環(huán)境變量(表6)。

表6 研究區(qū)AP地理加權(quán)回歸模型參數(shù)

3.4 土壤速效磷的空間變異特征分析

通過GS+9.0軟件對(duì)AP進(jìn)行地統(tǒng)計(jì)分析揭示其空間變異特征，通常當(dāng)C0/(C0+C)<25%時(shí)，變量具有較強(qiáng)的空間自相關(guān)性，且主要受地形等結(jié)構(gòu)性因子影響；當(dāng)25%75%時(shí)，變量空間自相關(guān)性較弱，主要受施肥等隨機(jī)因子的影響。本結(jié)果表明光榮小流域AP的變程為274 m，可作為該區(qū)域AP空間分布特征研究的參考采樣距離；C0/(C0+C)為32.5%，表明AP主要受結(jié)構(gòu)因子影響，同時(shí)在一定程度上受隨機(jī)因子影響，處于中等程度空間自相關(guān)水平。而海溝河流域C0/(C0+C)為7.2%，表明海溝河流域土壤AP空間自相關(guān)性較強(qiáng)，主要受結(jié)構(gòu)性因子影響，變程為1 650 m(表7)。

表7 研究區(qū)AP地統(tǒng)計(jì)學(xué)模型和參數(shù)

3.5 土壤速效磷的空間分布特征

不同空間插值方法所預(yù)測(cè)的光榮小流域AP空間分布格局整體趨勢(shì)基本一致(圖2)。AP含量在流域上游居民點(diǎn)和道路附近含量較高，并沿水系方向從上游向流域出口處遞減，與研究區(qū)高程變化情況較吻合。從制圖效果來看，反距離權(quán)重法和徑向基函數(shù)法存在明顯的斑塊效應(yīng)，但其預(yù)測(cè)值范圍更接近實(shí)測(cè)值變化范圍。多元線性回歸模型雖然能夠反映相似的AP空間分布格局，但是其預(yù)測(cè)值范圍存在明顯的低估現(xiàn)象，插值結(jié)果并不理想，相比之下，考慮了空間非平穩(wěn)性的地理加權(quán)回歸模型制圖效果與模擬精度均高于多元線性回歸模型。普通克里格、回歸克里格和地理加權(quán)回歸克里格插值范圍及空間預(yù)測(cè)格局相近，其中地理加權(quán)回歸克里格得到的預(yù)測(cè)值范圍更接近于實(shí)測(cè)值，同時(shí)考慮了環(huán)境變量的回歸克里格和地理加權(quán)回歸克里格的高低值界線趨于模糊化，地理加權(quán)回歸克里格還能夠揭示流域內(nèi)AP的局部分布特征。

海溝河流域AP低值區(qū)出現(xiàn)在流域東部林地種植區(qū)，在中部道路兩側(cè)出現(xiàn)明顯的條帶狀富集區(qū)，并沿水系方向在流域出口處達(dá)到最大值，其空間分布特征與高程變化及土地利用方式密切相關(guān)。從制圖效果看(圖3)，反距離權(quán)重法和徑向基函數(shù)法存在明顯的斑塊效應(yīng)，但其預(yù)測(cè)值范圍更接近實(shí)測(cè)值變化范圍。多元線性回歸模型和地理加權(quán)回歸模型出現(xiàn)明顯的低估現(xiàn)象。回歸克里格和地理加權(quán)回歸克里格模擬的AP空間分布情況更加精細(xì)化，能夠反映更多的空間信息。

3.6 模擬精度比較

以普通克里格法作為參照標(biāo)準(zhǔn)對(duì)比不同空間插值方法的評(píng)價(jià)參數(shù)，結(jié)果表明光榮小流域反距離權(quán)重法和徑向基函數(shù)法的各項(xiàng)驗(yàn)證指標(biāo)均低于普通克里格。多元線性回歸模型和地理加權(quán)回歸模型得到的均方根誤差較普通克里格方法分別降低了12%，2.1%，而回歸克里格和地理加權(quán)回歸克里格的各項(xiàng)驗(yàn)證指標(biāo)均優(yōu)于普通克里格，模擬精度分別提高了4.4%，4.6%。此外，引入主成分分析結(jié)果得到的插值方法評(píng)價(jià)參數(shù)均有一定程度的提高，其中引入主成分分析的回歸克里格方法模擬精度較普通克里格和回歸克里格分別提高5.7%，4.5%，引入主成分分析的地理加權(quán)回歸克里格方法模擬精度較普通克里格和地理加權(quán)回歸克里格分別提高了6.9%，2.4%(表8)。

在海溝河流域反距離權(quán)重法、徑向基函數(shù)法和協(xié)克里格的各項(xiàng)驗(yàn)證指標(biāo)均低于普通克里格法；而回歸克里格和地理加權(quán)回歸克里格的各項(xiàng)驗(yàn)證指標(biāo)均優(yōu)于普通克里格法，模擬精度分別提高了0.7%，1.7%；地理加權(quán)回歸模型較多元線性回歸模型的模擬精度提高了3.6%(表8)。綜上所述，考慮了環(huán)境變量和殘差空間自相關(guān)的回歸克里格和地理加權(quán)回歸克里格能夠提高AP的空間預(yù)測(cè)精度，地理加權(quán)回歸克里格法最優(yōu)。

注：A—L分別是反距離權(quán)重法、徑向基函數(shù)法、普通克里格、協(xié)同克里格、多元線性回歸模型、回歸克里格、地理加權(quán)回歸模型、地理加權(quán)回歸克里格、多元線性回歸模型(PCA)、回歸克里格(PCA)、地理加權(quán)回歸模型(PCA)、地理加權(quán)回歸克里格(PCA)。

注：A—H分別是反距離權(quán)重法、徑向基函數(shù)法、普通克里格、協(xié)同克里格、多元線性回歸模型、回歸克里格、地理加權(quán)回歸模型、地理加權(quán)回歸克里格。

表8 不同空間插值方法精度評(píng)價(jià)參數(shù)

4 討 論

通過對(duì)比AP空間分布圖和模擬精度評(píng)價(jià)指標(biāo)，我們發(fā)現(xiàn)反距離權(quán)重法和徑向基函數(shù)法制圖效果較差，存在明顯的斑塊效應(yīng)，同時(shí)，兩種方法的各項(xiàng)驗(yàn)證指標(biāo)均處于最低水平，研究表明，反距離權(quán)重法僅可在樣點(diǎn)密度大且地形起伏較小的情況下適用，最好用于較小區(qū)域的柵格化，不宜在較大范圍內(nèi)進(jìn)行插值[22]，徑向基函數(shù)多適用于樣本數(shù)目較多且變化平緩的數(shù)據(jù)[23]，Shen等[21]比較不同空間插值方法對(duì)黑土區(qū)土壤全磷模擬精度的影響，研究表明當(dāng)樣點(diǎn)空間分布不規(guī)則或均勻度較低時(shí)徑向基函數(shù)會(huì)產(chǎn)生精確的插值結(jié)果。因此，反距離權(quán)重法和徑向基函數(shù)法需在特殊情況下才可能應(yīng)用并提高黑土區(qū)小流域土壤AP的空間模擬精度。

普通克里格方法利用采樣點(diǎn)數(shù)據(jù)和半方差函數(shù)的結(jié)構(gòu)性，對(duì)未采樣點(diǎn)的區(qū)域化變量進(jìn)行最優(yōu)無偏估值，已被廣泛應(yīng)用于土壤養(yǎng)分空間異質(zhì)性分析研究[1]。本研究中，所選擇研究區(qū)樣點(diǎn)均勻度或樣點(diǎn)密度均較高，滿足普通克里格插值方法需要，其模擬精度僅次于回歸克里格和地理加權(quán)回歸克里格方法。

地理加權(quán)回歸模型是一種局部回歸模型，可用于處理目標(biāo)變量與解釋變量之間的空間非平穩(wěn)回歸系數(shù)[14-15]。本研究中，地理加權(quán)回歸模型模擬精度低于普通克里格方法，這可能是由于經(jīng)逐步線性回歸剔除后可用于建模的環(huán)境變量較少，回歸模型擬合精度較低。Wang等[13]比較了地理加權(quán)回歸模型對(duì)土壤有機(jī)質(zhì)和全氮模擬精度的影響，結(jié)果表明地理加權(quán)回歸模型能夠有效提高土壤有機(jī)質(zhì)的模擬精度，但是由于環(huán)境變量數(shù)據(jù)不詳盡且未能解釋所有相關(guān)變量之間的自相關(guān)和互相關(guān)而導(dǎo)致該方法對(duì)全氮模擬精度較差[24]。本研究已提取相關(guān)環(huán)境變量達(dá)30個(gè)，但并未進(jìn)一步提高回歸模型精度，這與Shen等[21]在黑土區(qū)土壤全磷空間模擬精度研究結(jié)果一致。此外，比較不同研究區(qū)回歸分析結(jié)果，并未發(fā)現(xiàn)共同的環(huán)境變量用于AP建模。因此，有關(guān)如何進(jìn)一步提高黑土區(qū)土壤磷素與環(huán)境變量之間的回歸模型擬合精度的問題需進(jìn)一步探討。研究表明，樣點(diǎn)密度會(huì)影響空間插值方法對(duì)土壤屬性法的模擬精度[19]，但本研究發(fā)現(xiàn)，當(dāng)樣點(diǎn)均勻度相似且高于45%時(shí)，樣點(diǎn)密度和采樣距離對(duì)AP模擬精度無顯著影響。

回歸克里格和地理加權(quán)回歸克里格已經(jīng)逐漸成為土壤屬性空間預(yù)測(cè)的優(yōu)秀空間統(tǒng)計(jì)方法，大量研究表明輔助變量的引入和殘差空間自相關(guān)的結(jié)合能夠有效提高土壤屬性的空間預(yù)測(cè)精度[9-10,16-19]，本研究結(jié)果表明，在不同尺度和景觀格局下，回歸克里格和地理加權(quán)回歸克里格各項(xiàng)驗(yàn)證指標(biāo)均優(yōu)于普通克里格插值方法，插值結(jié)果也能夠反映更多的空間信息，其中考慮了空間非平穩(wěn)性的地理加權(quán)回歸克里格方法可以作為黑土區(qū)小流域土壤AP空間預(yù)測(cè)的最優(yōu)插值方法。此外，我們發(fā)現(xiàn)引入主成分分析結(jié)果能夠有效提高回歸克里格方法和地理加權(quán)回歸克里格方法的模擬精度，因此，當(dāng)輔助變量充足且變量間存在多重共線性時(shí)，可以引入主成分分析以進(jìn)一步提高空間插值方法的模擬精度。

5 結(jié) 論

通過引入地形等環(huán)境因子建立的回歸克里格、地理加權(quán)回歸克里格方法較傳統(tǒng)空間插值方法能在一定程度上提高黑土小流域土壤AP的空間模擬精度，將主成分分析結(jié)果與上述方法相結(jié)合能夠進(jìn)一步提高模擬效果。其中，地理加權(quán)回歸克里格方法模擬效果最優(yōu)，可作為黑土小流域土壤AP空間分布模擬的主要方法。本研究選取的輔助變量以地形和植被覆蓋等環(huán)境因子為主，未能充分考慮與AP密切相關(guān)的施肥、耕作措施等農(nóng)田管理因子，可作為提高黑土小流域土壤AP空間模擬效果的重要方向。研究結(jié)果可為黑土小流域尺度土壤養(yǎng)分精準(zhǔn)管控提供技術(shù)支撐。