萬核洋 齊泓瑋 尚松浩

(清華大學(xué)水沙科學(xué)與水利水電工程國家重點(diǎn)實(shí)驗(yàn)室， 北京 100084)

0 引言

土壤是與人類關(guān)系最為密切的環(huán)境要素之一，其中的元素遷移、分散和累積受到地域特征的制約，既表現(xiàn)出垂向和區(qū)域的結(jié)構(gòu)性，又具有時(shí)空上的漸變性[1]。而土壤質(zhì)地作為土壤的重要物理屬性，在母質(zhì)、氣候、地形、人類活動(dòng)等因素的影響下，也會在空間上呈現(xiàn)出特定的規(guī)律和結(jié)構(gòu)性。研究土壤質(zhì)地的空間分布規(guī)律，對指導(dǎo)土壤改良、灌溉排水管理、生態(tài)農(nóng)業(yè)區(qū)劃以及水土資源合理利用具有重要的現(xiàn)實(shí)意義，進(jìn)而為灌區(qū)農(nóng)業(yè)可持續(xù)發(fā)展奠定基礎(chǔ)[1-2]。

土壤質(zhì)地的空間變異性分析和插值方法有很多，其中基于變異函數(shù)理論的克里金法應(yīng)用最為廣泛，但是克里金法往往會產(chǎn)生平滑效應(yīng)，即高估較小值卻低估較大值。楊雨亭等[3]研究發(fā)現(xiàn)，克里金平滑效應(yīng)會導(dǎo)致結(jié)果無法真實(shí)反映土壤水分空間分布的局部特征，李江等[4]在地下水位的克里金插值研究中發(fā)現(xiàn)平滑效應(yīng)對極高值與極低值點(diǎn)的準(zhǔn)確預(yù)測產(chǎn)生了較大的影響。

決定土壤質(zhì)地類型的顆粒組成屬于典型的成分?jǐn)?shù)據(jù)，在插值中需額外滿足各成分非負(fù)且和為常數(shù) (1或100%)的約束條件。AITCHISON[5]提出了對數(shù)比轉(zhuǎn)換的方法用于解決成分?jǐn)?shù)據(jù)的閉合問題，后來有學(xué)者又進(jìn)行了大量改進(jìn)和研究[6-7]。對數(shù)比轉(zhuǎn)換方法有很多種類，主要包括加和對數(shù)比(ALR)轉(zhuǎn)換[5]、中心對數(shù)比(CLR)轉(zhuǎn)換[5]、等距對數(shù)比(ILR)轉(zhuǎn)換[6]和對稱對數(shù)比(SLR)轉(zhuǎn)換[7]。已有研究表明，與直接對原始數(shù)據(jù)插值相比，對于經(jīng)過對數(shù)比轉(zhuǎn)換的土壤顆粒組成進(jìn)行插值的結(jié)果更加合理，精度也更高[8-11]。WAN等[8]研究表明ALR轉(zhuǎn)換對克里金插值結(jié)果的平滑性產(chǎn)生了一定影響。但目前不同轉(zhuǎn)換方法對克里金插值平滑性的影響分析研究較少。

此外，關(guān)于非成分?jǐn)?shù)據(jù)克里金平滑效應(yīng)的專門處理方法也有很多研究。隨機(jī)模擬是重現(xiàn)區(qū)域化變量空間結(jié)構(gòu)的方法，可以避免平滑效應(yīng)，但是無法保證局部估計(jì)精度且結(jié)果不唯一[12-13]。OLEA等[13]提出的補(bǔ)償克里金法全局最優(yōu)性不如隨機(jī)模擬，局部精度亦低于傳統(tǒng)的克里金估計(jì)。JOURNEL等[14]的后處理方法較好解決了平滑效應(yīng)問題，但也降低了局部估計(jì)精度。YAMAMOTO[15-16]提出的后處理方法會導(dǎo)致插值表面光滑性降低，而且仍局限在克里金法的理論體系內(nèi)。段建軍等[17]提出的極差平滑率和標(biāo)準(zhǔn)差平滑率等定量指標(biāo)較好地評估了克里金的平滑效應(yīng)，并且提出了簡單易行的一元二次回歸模型校正方法，但也存在特征點(diǎn)少、模型穩(wěn)定性不夠的缺點(diǎn)。目前，在土質(zhì)插值中考慮組分?jǐn)?shù)據(jù)克里金插值平滑效應(yīng)及其校正的研究較少。

BP神經(jīng)網(wǎng)絡(luò)是一種應(yīng)用最廣泛的人工神經(jīng)網(wǎng)絡(luò)模型，結(jié)構(gòu)簡單，適用性強(qiáng)，具有較強(qiáng)的非線性函數(shù)逼近能力，在數(shù)據(jù)預(yù)測、樣本分類等研究中均有較好的精度及應(yīng)用價(jià)值[18-21]。在克里金插值中驗(yàn)證點(diǎn)的殘差受到周圍觀測點(diǎn)及其克里金權(quán)重的影響，這也是產(chǎn)生平滑效應(yīng)的原因，而決定克里金權(quán)重的變異函數(shù)一般是非線性模型，因此基于BP模型尋找殘差與鄰近點(diǎn)的關(guān)系就有了可能。

本文以內(nèi)蒙古河套灌區(qū)的土壤質(zhì)地為研究對象，基于普通克里金法進(jìn)行空間插值，分析評估不同的對數(shù)比轉(zhuǎn)換方法對克里金插值平滑效應(yīng)和插值精度的影響，進(jìn)一步建立基于BP神經(jīng)網(wǎng)絡(luò)的后處理模型，以期能夠較好地改善克里金插值結(jié)果的平滑效應(yīng)，提升插值精度，并得到較為合理的研究區(qū)土壤質(zhì)地空間分布。

1 材料與方法

1.1 研究區(qū)概況

河套灌區(qū)位于內(nèi)蒙古自治區(qū)西部的巴彥淖爾市境內(nèi)(40.1°～41.4° N，106.1°～109.4° E)，是我國干旱區(qū)的最大灌區(qū)，也是亞洲最大的一首制灌區(qū)，東西長約500 km，南北寬20～90 km，灌區(qū)占地 1.21×106hm2，灌溉面積約7.3×105hm2(圖1)。

河套灌區(qū)地處典型的中溫帶大陸性季風(fēng)氣候區(qū)，蒸發(fā)強(qiáng)烈，降水稀少。年平均氣溫6～8℃，自東向西升高；平均相對濕度40%～50%；年日照時(shí)數(shù)為3 100～3 300 h；年平均降水量為150～200 mm，自東向西遞減；20 cm蒸發(fā)皿年蒸發(fā)量 1 900～2 300 mm。河套灌區(qū)種植結(jié)構(gòu)十分復(fù)雜，作物呈插花式、細(xì)碎化分布。其中，春玉米、向日葵和春小麥?zhǔn)呛犹坠鄥^(qū)的主要作物，占總種植面積的85%以上[22]。

河套灌區(qū)是由于黃河攜帶大量泥沙在狼山以南的陷落地帶沖刷堆積而形成的沖積平原，土地肥沃，整體上地勢西南高、東北低，高程在1 000～1 090 m之間，該地形特征為自流式引水灌溉提供了天然條件。灌區(qū)土壤類型主要以鹽漬化淺色草甸土、鹽土和灌淤土為主，質(zhì)地偏砂，有機(jī)質(zhì)含量約1%，含鹽量較高。

1.2 土樣采集與顆粒組成測定

本研究采用網(wǎng)格法和隨機(jī)法相結(jié)合的方式進(jìn)行點(diǎn)位設(shè)計(jì)，使得采樣點(diǎn)不僅能均勻覆蓋全灌區(qū)，而且采樣點(diǎn)的間距也有大有小，共計(jì)采樣298個(gè)(圖1)。土樣采集時(shí)間為2018年5月、2019年9月和2021年5月，采樣點(diǎn)主要分布于農(nóng)田內(nèi)，利用手持GPS進(jìn)行定位確定采樣點(diǎn)的經(jīng)緯度，采集表層土 0～20 cm 內(nèi)的混合土樣作為待測樣品進(jìn)行室內(nèi)分析。

采用篩分法和比重計(jì)法測定土壤顆粒組成，根據(jù)美國農(nóng)業(yè)部土壤質(zhì)地分類法，利用土壤黏粒(Clay)、粉粒(Silt)和砂粒(Sand)的含量(質(zhì)量分?jǐn)?shù))，將土壤質(zhì)地類型分為黏土、砂質(zhì)黏土、粉質(zhì)黏土、粉質(zhì)黏壤土、黏壤土、砂質(zhì)黏壤土、壤土、粉壤土、粉土、砂質(zhì)壤土、壤質(zhì)砂土、砂土。由于灌區(qū)內(nèi)土質(zhì)不會在短期內(nèi)發(fā)生變化，因此3年的采樣數(shù)據(jù)一起處理和統(tǒng)計(jì)分析，通過R語言編程以及GS+ 9.0軟件進(jìn)行。

1.3 研究方法

1.3.1成分?jǐn)?shù)據(jù)的對數(shù)比轉(zhuǎn)換方法

AITCHISON[5]最早提出成分?jǐn)?shù)據(jù)的對數(shù)比轉(zhuǎn)換方法，用以解決成分?jǐn)?shù)據(jù)在插值過程中的閉合效應(yīng)和統(tǒng)計(jì)分析問題，后續(xù)又經(jīng)過不同學(xué)者改進(jìn)與發(fā)展。目前，對數(shù)比轉(zhuǎn)換方法主要包括加和對數(shù)比(ALR)、中心對數(shù)比(CLR)、等距對數(shù)比(ILR)和對稱對數(shù)比(SLR)4種轉(zhuǎn)換方法，使用時(shí)均是先將各顆粒組分通過轉(zhuǎn)換公式進(jìn)行變換，然后對變換后的數(shù)據(jù)進(jìn)行插值等操作，最后采用轉(zhuǎn)回公式將結(jié)果轉(zhuǎn)換為原始的成分?jǐn)?shù)據(jù)形式，這樣的變換過程能夠保證插值后的土壤顆粒各組分總和為100%，滿足定和、非負(fù)等約束。

加和對數(shù)比(ALR)轉(zhuǎn)換公式和轉(zhuǎn)回公式為[5]

yi=ln(zi/zD) (i=1,2,…,D-1)

(1)

其中

(2)

式中zi——第i種顆粒的質(zhì)量分?jǐn)?shù)(依次是黏粒、粉粒和砂粒)

yi——轉(zhuǎn)換后的顆粒i質(zhì)量分?jǐn)?shù)

yj——轉(zhuǎn)換后的顆粒j質(zhì)量分?jǐn)?shù)

D——顆粒組分種類數(shù)量，取3

中心對數(shù)比(CLR)轉(zhuǎn)換公式和轉(zhuǎn)回公式為[5]

(3)

(4)

等距對數(shù)比(ILR)轉(zhuǎn)換公式和轉(zhuǎn)回公式為[6]

(5)

其中

(6)

對稱對數(shù)比(SLR)轉(zhuǎn)換公式和轉(zhuǎn)回公式為[7]

(7)

其中

(i=1, 2, …,D)

(8)

式中δi——樣點(diǎn)中第i種顆粒最小質(zhì)量分?jǐn)?shù)(除0外)的1/2

δj——樣點(diǎn)中第j種顆粒最小質(zhì)量分?jǐn)?shù)(除0外)的1/2

1.3.2普通克里金法

克里金法是地質(zhì)統(tǒng)計(jì)學(xué)的主要內(nèi)容之一。普通克里金(Ordinary Kriging, OK)是最早被提出和系統(tǒng)研究的克里金法，它實(shí)質(zhì)上是一個(gè)線性估計(jì)系統(tǒng)，基本原理是利用半變異函數(shù)模型來擬合經(jīng)驗(yàn)半變異函數(shù)，通過無偏估計(jì)和最優(yōu)估計(jì)條件建立克里金線性方程式，求解各樣本線性權(quán)重，最后用各樣本及其權(quán)重表示未知點(diǎn)估計(jì)值[23]。

在二階平穩(wěn)假設(shè)條件下，對土壤顆粒組分的普通克里金估計(jì)公式為

(9)

式中Z(x0)——位置x0處的土壤顆粒組分估計(jì)值

Z(xi)——已知位置xi處的土壤顆粒組分觀測值

λi——權(quán)重系數(shù)

n——用于估計(jì)的已知樣點(diǎn)個(gè)數(shù)

在無偏和最優(yōu)的約束條件下，可得克里金方程組

(10)

式中μ——拉格朗日乘數(shù)

γ——變異函數(shù)

通過式(9)、(10)求解，即可得到權(quán)重系數(shù)λ，進(jìn)而求得估計(jì)值。插值結(jié)果的精度采用交叉驗(yàn)證法評價(jià)，即對任意一個(gè)采樣點(diǎn)，利用除該點(diǎn)之外的采樣數(shù)據(jù)進(jìn)行插值得到該點(diǎn)的插值結(jié)果，其與實(shí)測值的差值作為誤差，然后計(jì)算1.3.5節(jié)的精度評價(jià)指標(biāo)。

1.3.3BP神經(jīng)網(wǎng)絡(luò)平滑性校正模型

BP神經(jīng)網(wǎng)絡(luò)是一種信號向前傳遞、誤差反向傳播的有監(jiān)督人工神經(jīng)網(wǎng)絡(luò)，適用于非線性分類和回歸問題，由輸入層、隱含層和輸出層組成。

輸入層的特征變量為插值點(diǎn)周圍的觀測值序列。首先確定待插值點(diǎn)的位置，然后通過鄰近法尋找距離插值點(diǎn)最近的n個(gè)觀測點(diǎn)，這里的n與克里金插值時(shí)所用的插值點(diǎn)數(shù)量n保持一致，并將這n個(gè)點(diǎn)從近到遠(yuǎn)排序；在交叉驗(yàn)證中，每個(gè)采樣點(diǎn)依次進(jìn)行類似操作，通過匯總，可以得n個(gè)向量，每個(gè)向量是由距離插值點(diǎn)遠(yuǎn)近排序相同的觀測點(diǎn)組成，其元素個(gè)數(shù)為采樣點(diǎn)總數(shù)。比如對于不同插值點(diǎn)，對應(yīng)距離最近的點(diǎn)組成第1個(gè)特征，第2近的點(diǎn)組成第2個(gè)特征，依次類推，可以得到n個(gè)特征向量。

其次確定輸出變量。在交叉驗(yàn)證中，獲得每個(gè)點(diǎn)的克里金插值結(jié)果后，計(jì)算插值殘差，以此作為神經(jīng)網(wǎng)絡(luò)模型的輸出變量。

根據(jù)交叉驗(yàn)證的結(jié)果，得到訓(xùn)練數(shù)據(jù)集，經(jīng)過反復(fù)調(diào)整參數(shù)得到效果較好的BP網(wǎng)絡(luò)模型。對于某個(gè)位置的插值結(jié)果，將插值點(diǎn)周圍觀測值特征序列輸入模型，即可估算該插值點(diǎn)的殘差，然后結(jié)合插值結(jié)果，得到校正后的最終插值結(jié)果。

1.3.4平滑性指標(biāo)

克里金插值的平滑效應(yīng)是指插值后原來的高值被低估，而原來的低值被高估，呈現(xiàn)一種趨中性，因此本研究采用極差平滑率、標(biāo)準(zhǔn)差平滑率指標(biāo)來衡量平滑效應(yīng)[17]，計(jì)算式為

(11)

(12)

式中CR——極差平滑率

CS——標(biāo)準(zhǔn)差平滑率

Ro——采樣點(diǎn)觀測值極差

So——采樣點(diǎn)觀測值標(biāo)準(zhǔn)差

RI——采樣點(diǎn)估計(jì)值極差

SI——采樣點(diǎn)估計(jì)值標(biāo)準(zhǔn)差

其中，CR和CS越大，說明平滑性越嚴(yán)重。

1.3.5精度評價(jià)指標(biāo)

本研究在克里金插值中選擇合適的變異函數(shù)模型，需要通過交叉驗(yàn)證法進(jìn)行插值精度評估；其次在插值后構(gòu)建BP神經(jīng)網(wǎng)絡(luò)平滑性校正模型時(shí)，需要評估模型精度。采用平均絕對誤差(Mean absolute error, MAE)和均方根誤差(Root mean squared error, RMSE)作為精度指標(biāo)來評價(jià)模型，通常兩個(gè)指標(biāo)的值越小，模型的效果越好。

2 結(jié)果與分析

2.1 描述性統(tǒng)計(jì)分析

圖2為土壤實(shí)測樣本顆粒組成的小提琴分布圖。砂粒、粉粒和黏粒表現(xiàn)出不同的分布特點(diǎn)。其中黏粒含量平均值為15.3%，與中位數(shù)接近；粉粒含量平均值為40.5%，低于相應(yīng)的中位數(shù)；砂粒含量平均值為44.2%，高于其中位數(shù)。從數(shù)據(jù)的分布形狀來看，黏粒分布的對稱性更好，粉粒向高值一側(cè)偏斜，而砂粒向低值一側(cè)偏斜。黏粒含量相對較低，分布集中，粉粒和砂粒含量整體較高，說明土壤黏粒不是該地區(qū)主要的控制性組分。另外粉粒和砂粒的分布范圍較廣，幾乎在0～100%內(nèi)均存在，說明其變異性也較大。從極端值分布來看，黏粒中相對來說存在較多，在高值和低值區(qū)均存在，粉粒和砂粒的極端值相對較少。極端值的存在會影響插值結(jié)果的平滑性，一般來說，極端值越多平滑性越嚴(yán)重。在實(shí)際插值中，往往會將這些點(diǎn)視作異常值剔除，但是灌區(qū)中的確有沙地存在，為了反映客觀事實(shí)和對插值平滑性的探索，本研究保留這些極端值。

圖2 土壤實(shí)測樣本顆粒組成分布Fig.2 Distribution of soil particle-size fractions for observed samples

圖3 為對數(shù)比轉(zhuǎn)換前后的土壤組分?jǐn)?shù)據(jù)的正態(tài)性分布。在轉(zhuǎn)換之前，粉粒相對來說更接近正態(tài)性，黏粒由于兩邊極端值較多而偏離正態(tài)曲線較多，砂粒中高值部分偏離了正態(tài)曲線。原始土壤組分?jǐn)?shù)據(jù)經(jīng)過不同的對數(shù)比轉(zhuǎn)換后，在正態(tài)性上表現(xiàn)各有差異。ALR轉(zhuǎn)換后的第1成分ALR1除了兩端少數(shù)點(diǎn)外的大部分點(diǎn)都在正態(tài)曲線上，第2成分ALR2除了較大的極端值偏離較多外，也基本都在正態(tài)曲線上。CLR和SLR轉(zhuǎn)換后的3個(gè)成分均各自對應(yīng)于原始數(shù)據(jù)的黏粒、粉粒和砂粒組分，CLR1和SLR1表現(xiàn)類似，都是低值部分偏離較多，其他大部分點(diǎn)都接近標(biāo)準(zhǔn)線；CLR2和SLR2也是低值部分偏離較多；而CLR3和SLR3則是高值部分偏離較多，主要是因?yàn)樵紨?shù)據(jù)中黏粒和粉粒均存在極端的低值，而砂粒相應(yīng)存在極端的高值。對于ILR轉(zhuǎn)換，其中組分ILR1與ALR1類似，大部分值在標(biāo)準(zhǔn)線附近，而ILR2與ALR2不同，主要是極端的低值部分偏離了標(biāo)準(zhǔn)線。

圖3 原始樣本及對數(shù)比轉(zhuǎn)換后數(shù)據(jù)的正態(tài)分位圖Fig.3 Normal quantile-quantile plots for original and logratio transformed data

為了進(jìn)一步說明正態(tài)性的變化，通過Kolmogorov- Smirnov(KS)檢驗(yàn)得到原始樣本和不同轉(zhuǎn)換后的數(shù)據(jù)的p值結(jié)果，如表1所示。在剔除極端值之前，對數(shù)比轉(zhuǎn)換前后差別不大，遠(yuǎn)不能滿足正態(tài)假設(shè)；但剔除極端值以后，經(jīng)過對數(shù)比轉(zhuǎn)換的數(shù)據(jù)p值更接近甚至大于0.05，正態(tài)性明顯得到了改善。綜上可知，原始數(shù)據(jù)中由于存在離群值的情況，很難完全符合標(biāo)準(zhǔn)正態(tài)分布，即使經(jīng)過對數(shù)比轉(zhuǎn)換也很難使所有數(shù)據(jù)都靠近標(biāo)準(zhǔn)線，但是可以使除極端值之外的大部分?jǐn)?shù)據(jù)更加接近正態(tài)分布。

表1 原始樣本及不同對數(shù)比轉(zhuǎn)換后數(shù)據(jù)的KS檢驗(yàn)的p值結(jié)果Tab.1 p values by KS test for original and logratio transformed data

2.2 基于原始數(shù)據(jù)的OK插值

首先使用普通克里金法直接對原始的土壤顆粒組分?jǐn)?shù)據(jù)進(jìn)行插值計(jì)算，然后對交叉驗(yàn)證結(jié)果進(jìn)行分析和評估。如圖4所示，對原始數(shù)據(jù)直接使用OK插值的結(jié)果具有較明顯的平滑性，平滑率基本都在0.3以上，其中粉粒的平滑性最大，極差平滑率和標(biāo)準(zhǔn)差平滑率均在0.4以上，黏粒和砂粒的平滑性較為接近，極差平滑率在0.3左右，標(biāo)準(zhǔn)差平滑率在0.35左右。

圖4 土壤顆粒組成原始數(shù)據(jù)及對數(shù)比轉(zhuǎn)換數(shù)據(jù)的插值平滑率Fig.4 Smoothing rates for interpolation of original soil particle size fractions and logratio transformed data

圖5為原始數(shù)據(jù)的OK插值精度評價(jià)結(jié)果。從圖5可以看出，不同組分的插值精度具有明顯差異，其中黏粒的插值精度最高，平均絕對誤差和均方根誤差分別為4.4%和5.6%。插值精度最差的是砂粒，兩個(gè)指標(biāo)均大于12%，而粉粒精度介于黏粒和砂粒之間，平均絕對誤差為10.5%，均方根誤差為12.9%。

圖5 不同方法插值的交叉驗(yàn)證精度Fig.5 Comparisons of accuracy indexes of cross-validation for different methods

2.3 基于不同對數(shù)比轉(zhuǎn)換的插值平滑性評估

原始組分?jǐn)?shù)據(jù)經(jīng)過4種不同的對數(shù)比轉(zhuǎn)換后，再用普通克里金插值，然后再使用轉(zhuǎn)回公式得到最終結(jié)果，通過交叉驗(yàn)證結(jié)果評估插值的平滑性(圖4)。與原始數(shù)據(jù)的直接插值結(jié)果不同，經(jīng)過對數(shù)比轉(zhuǎn)換后的插值結(jié)果中，砂粒的平滑性得到了明顯改善，極差平滑率由0.3降低到0.2以下，標(biāo)準(zhǔn)差平滑率也由0.38下降到0.15左右。黏粒的插值平滑性改善效果相對較差，其中極差平滑率甚至大于原始數(shù)據(jù)，在3個(gè)組分中平滑性最嚴(yán)重。粉粒則位于兩者之間，其極差平滑率由0.45下降到0.3左右，標(biāo)準(zhǔn)差平滑率也由0.43下降到0.25以下。

不同的轉(zhuǎn)換方法對平滑性的影響程度也有所不同。從極差平滑率來看，SLR轉(zhuǎn)換結(jié)果的平滑性最大，ALR和CLR相差不大，而ILR轉(zhuǎn)換結(jié)果的平滑性最小，其中黏粒和粉粒分別為0.29和0.30，砂粒僅為0.16。對標(biāo)準(zhǔn)差平滑率而言，規(guī)律類似，ILR具有最小的平滑性，黏粒、粉粒和砂粒的標(biāo)準(zhǔn)差平滑率依次為0.22、0.20和0.15?？傮w而言，ILR轉(zhuǎn)換方法在普通克里金插值中對平滑性的改善效果最好。與直接使用普通克里金插值的結(jié)果相比, 黏粒、粉粒和砂粒的極差平滑率分別減小5.8%、33.8%和45.6%，標(biāo)準(zhǔn)差平滑率分別減小38.6%、53.9%和60.2%。另外從圖5中可以看出，基于ILR轉(zhuǎn)換數(shù)據(jù)的插值(ILR-OK)對各顆粒組分的插值精度影響不大，甚至比直接使用OK還略差，其主要原因是對數(shù)比轉(zhuǎn)換中考慮了組分?jǐn)?shù)據(jù)的特性，相當(dāng)于增加了OK插值的約束條件。

2.4 基于ILR轉(zhuǎn)換和BP神經(jīng)網(wǎng)絡(luò)的平滑效應(yīng)校正模型(ILR-BP)

與基于原始數(shù)據(jù)的普通克里金插值相比，對ILR轉(zhuǎn)換后的數(shù)據(jù)插值可以改善平滑性，但還依然存在較大的平滑率，尤其是黏粒和粉粒仍具有接近0.3的極差平滑率，因此需要進(jìn)一步修正以降低平滑性。在ILR-OK的交叉驗(yàn)證結(jié)果中，將目標(biāo)點(diǎn)的插值殘差作為訓(xùn)練目標(biāo)，同時(shí)將目標(biāo)點(diǎn)插值與周圍鄰近點(diǎn)觀測值的殘差構(gòu)成的特征向量作為輸入，通過構(gòu)建BP神經(jīng)網(wǎng)絡(luò)模型擬合兩者之間的關(guān)系，不斷調(diào)試參數(shù)，將具有較好訓(xùn)練效果的模型作為最終的殘差預(yù)測模型。

圖6為3種顆粒組分的殘差預(yù)測BP模型的擬合效果。由圖6可見，不同組分的殘差分布不同，黏粒的殘差最小，粉粒次之，而砂粒的殘差較大。另外，散點(diǎn)基本均勻分布在1∶1線兩邊，黏粒的擬合效果相對較好，平均絕對誤差和均方根誤差分別為3.24%和4.61%；粉粒和砂粒由于原始的殘差比較大，變異性強(qiáng)，因此整體擬合效果不如黏粒好。大部分散點(diǎn)主要集中在0值附近，說明大部分插值點(diǎn)的平滑性不明顯，而少數(shù)分布在兩端的點(diǎn)是影響平滑性的主要因素，也是決定平滑性修正效果的關(guān)鍵點(diǎn)?？傮w可見，殘差預(yù)測模型的擬合效果在可以接受的范圍內(nèi)，尤其是對于較大的殘差，即平滑性嚴(yán)重的樣點(diǎn)，有較強(qiáng)的擬合能力。

圖6 基于ILR-BP的殘差預(yù)測模型評估結(jié)果Fig.6 Evaluation results of residual prediction model based on ILR-BP

基于殘差預(yù)測模型的預(yù)測結(jié)果，對ILR轉(zhuǎn)換數(shù)據(jù)的交叉驗(yàn)證結(jié)果進(jìn)行校正。根據(jù)平滑性校正前后觀測值與交叉驗(yàn)證結(jié)果的對比(圖7)可以發(fā)現(xiàn)，校正后的極差平滑率均為0，并且標(biāo)準(zhǔn)差平滑率也接近0(0.03 ～ 0.07)，說明平滑性得到了很大改善，相比只使用ILR轉(zhuǎn)換，平滑性得到了進(jìn)一步明顯的削弱。其中粉粒的平滑性效果改善最大，從原來的0.4以上直接減小為0，而黏粒和砂粒的標(biāo)準(zhǔn)差平滑率比粉粒稍大，均為0.07，但總體上平滑性也已很小。從散點(diǎn)分布來看，校正后的散點(diǎn)比校正前分布更加離散，但是更加均勻地分布在1∶1線的兩邊，靠近均值的點(diǎn)變化不大，高值區(qū)和低值區(qū)的點(diǎn)則均有向1∶1線靠攏的趨勢，更說明了校正模型對極端值平滑性的改善效果。從平滑性校正前后的插值精度對比(圖5)可以發(fā)現(xiàn)，經(jīng)過ILR-BP模型校正后各組分的插值精度均有所提高。黏粒的平均絕對誤差和均方根誤差分別減小29.4%和19.8%；粉粒的精度提高也較為明顯，兩個(gè)指標(biāo)分別減小27.6%和21.0%；砂粒的精度提升效果不如黏粒和粉粒，兩個(gè)指標(biāo)分別減小17.1%和14.6%。因此，校正后的結(jié)果在平滑性得到極大改善的同時(shí)，還能保證插值精度有所提升。

圖7 平滑性校正前后觀測值與交叉驗(yàn)證結(jié)果對比Fig.7 Comparisons between observations and cross-validation results before and after smoothing correction

由于土壤顆粒各組分是單獨(dú)進(jìn)行插值，平滑性會直接體現(xiàn)在各組分自身上，但土壤質(zhì)地類型是黏粒、粉粒和砂粒的含量共同決定的，因而各組分的插值平滑性會進(jìn)一步影響到土質(zhì)類型的預(yù)測結(jié)果。為了進(jìn)一步比較ILR-BP模型校正后與校正前的土質(zhì)類型變化情況，根據(jù)美國制標(biāo)準(zhǔn)，通過土壤質(zhì)地三角來展示原始樣本數(shù)據(jù)以及插值平滑性校正前后的交叉驗(yàn)證結(jié)果，如圖8所示。由圖8a可知，原始樣本數(shù)據(jù)中，共有7種土質(zhì)，主要為壤質(zhì)砂土、砂質(zhì)壤土、壤土和粉壤土，另外還有少量的砂土、黏壤土和粉質(zhì)黏壤土。直接經(jīng)過OK插值的交叉驗(yàn)證結(jié)果卻只有4種主要的土質(zhì)(圖8b)，正是因?yàn)椴逯档钠交?yīng)，各組分的結(jié)果都向中部集中，使得原始樣本中存在的3種土質(zhì)消失了。為了避免這種現(xiàn)象的發(fā)生，需要進(jìn)行平滑性校正，校正后的土質(zhì)分布與原始樣本高度一致(圖8c)，首先是分布的離散程度，相比校正前的趨中性，校正后明顯更加分散；其次土質(zhì)類型也包括了原始樣本中的所有7種土質(zhì)，但是也出現(xiàn)了個(gè)別原始樣本中沒有的砂質(zhì)黏壤土，主要是因?yàn)樵紭颖局写嬖诳拷百|(zhì)黏壤土與砂質(zhì)壤土界限的一部分樣點(diǎn)，分類的不確定性大。因此，通過土質(zhì)類型的對比，可以看出ILR-BP模型校正后的結(jié)果能更準(zhǔn)確地描述原始樣本中的土壤質(zhì)地類型。

圖8 原始樣本點(diǎn)及平滑性校正前后交叉驗(yàn)證結(jié)果的土質(zhì)類型對比Fig.8 Comparisons of soil texture types for original samples and cross-validation results before and after smoothing correction

2.5土壤顆粒組成及質(zhì)地空間分布基于ILR-BP模型可以對研究區(qū)內(nèi)表層土壤顆粒組分的空間分布插值結(jié)果進(jìn)行校正，如圖9所示，圖中各組分的范圍是根據(jù)原始樣本數(shù)據(jù)設(shè)定的。圖9a為直接用OK插值后沒有經(jīng)過校正的結(jié)果，由于平滑效應(yīng)，插值的預(yù)測值都在相應(yīng)閾值內(nèi)，而且基本無法顯示出極端值的分布區(qū)域。圖9b為經(jīng)過平滑性校正后的空間分布，黏粒在灌區(qū)中部和西部含量較低，在這兩個(gè)地方也確實(shí)存在黏粒含量幾乎為0的測量樣本，但是校正前的黏粒分布中幾乎不能體現(xiàn)這一特點(diǎn)。另外灌區(qū)大部分地區(qū)黏粒含量在15%～20%，東部存在零散的地區(qū)偏高。粉粒的分布特點(diǎn)與黏粒類似，校正后的分布更加精細(xì)，局部特征更明顯，尤其是存在極端低值的樣本區(qū)域，大部分地區(qū)的粉粒含量在40%左右。對于砂粒而言，幾乎與黏粒和粉粒的分布趨勢相反，主要是因?yàn)槌煞謹(jǐn)?shù)據(jù)的定和約束，砂粒含量在灌區(qū)西部和中部偏南的部分地區(qū)較高，這也與實(shí)際采樣的結(jié)果一致。校正前的砂粒分布雖然也有此規(guī)律，但是原本的極端高值由于平滑效應(yīng)而消失了，另外灌區(qū)東部存在砂粒含量較低的區(qū)域，大部分區(qū)域的砂粒含量在50%以下?？傮w而言，校正后的顆粒含量空間分布更加明晰準(zhǔn)確，局部變異特征更加明顯，更能反映采樣數(shù)據(jù)的空間特征。基于土壤顆粒組分的空間分布，根據(jù)美國制標(biāo)準(zhǔn)，可以計(jì)算出相應(yīng)的土質(zhì)類型的空間分布，如圖10所示。由圖10a可以看出，主要有4種土質(zhì)類型，即壤質(zhì)砂土、砂質(zhì)壤土、壤土和粉壤土，與原始樣本的主要土質(zhì)一致，但是沒有反映出其余3種分布較少的土質(zhì)類型的采樣點(diǎn)特征，這正是由于平滑效應(yīng)造成的結(jié)果。由圖10b可知，灌區(qū)中部和西部主要是偏砂性的土壤，與校正前相比，還顯示出了砂土的分布特征；壤土和粉壤土的面積最廣，其中壤土主要集中分布在中西部、東北部和東南部，粉壤土主要分布在北部和南部靠近黃河的地區(qū)；除此之外，灌區(qū)內(nèi)還存在零星分布的黏壤土和粉質(zhì)黏壤土，這與采樣點(diǎn)的土質(zhì)類型基本一致。因此，經(jīng)過平滑性校正后，土質(zhì)的空間分布既能保證整體規(guī)律不變，又能反映更加真實(shí)的土質(zhì)類型，呈現(xiàn)出更加明顯的交錯(cuò)分布特征。

圖9 平滑性校正前后土壤顆粒組成空間分布Fig.9 Spatial distributions of soil particle size fractions before and after smoothing correction

圖10 平滑性校正前后土壤質(zhì)地類型空間分布Fig.10 Spatial distributions of soil texture types before and after smoothing correction

在河套灌區(qū)，平均2～3年農(nóng)田就會深耕一次，土壤經(jīng)過上下的頻繁翻動(dòng)，可以認(rèn)為作物根系區(qū)的土壤質(zhì)地較為均勻，研究中采用表層0～20 cm土樣具有代表性。根據(jù)ILR-BP模型校正之后的河套灌區(qū)各土質(zhì)類型面積及占比情況，如表2所示，壤土和粉壤土的面積占比分別達(dá)到41.7%和27.3%，二者覆蓋了2/3以上的灌區(qū)面積；砂質(zhì)壤土次之，占比15.8%；砂土和壤質(zhì)砂土不到10%；其余3種占比很小，共占1.4%左右。

表2 河套灌區(qū)各土壤質(zhì)地類型面積及相應(yīng)占比Tab.2 Areas and proportions of different soil texture types in Hetao Irrigation District

3 討論

3.1 插值平滑效應(yīng)與插值精度的關(guān)系

克里金插值中的平滑效應(yīng)主要指高值低估和低值高估現(xiàn)象，為了定量表示，采用段建軍等[17]提出的極差平滑率和標(biāo)準(zhǔn)差平滑率，從數(shù)據(jù)離散程度的變化反映平滑效應(yīng)。插值精度是用來衡量插值結(jié)果整體準(zhǔn)確度的指標(biāo)，本研究中采用使用廣泛的平均絕對誤差和均方根誤差來表示。平滑效應(yīng)和插值精度之間存在一定聯(lián)系，但不完全同步變化。

以形式較為簡單的極差平滑率和平均絕對誤差為例。其中極差平滑率主要是通過最大值和最小值的插值前后變化來計(jì)算，本質(zhì)上反映了極端值的插值精度。平均絕對誤差的計(jì)算是把所有樣本點(diǎn)的插值誤差都考慮在內(nèi)，當(dāng)然也包括極端值的誤差。從這個(gè)角度看，極差平滑率是構(gòu)成平均絕對誤差的一部分。如果經(jīng)過插值方法改進(jìn)，極端值的插值誤差雖然減小但幅度不大，而其他樣本點(diǎn)的總誤差增大，那么就有可能導(dǎo)致平滑效應(yīng)減小而插值精度卻沒有提高的現(xiàn)象，如ILR-OK的插值結(jié)果所示。但如果極端值的插值誤差減小的同時(shí)，其他樣本點(diǎn)的誤差也都不同程度地減小，那么平滑效應(yīng)和插值精度就同時(shí)得到改善，如ILR-BP的插值結(jié)果所示。因此平滑率并不能完全決定插值精度，兩者也不存在絕對的同步變化關(guān)系。

3.2 對數(shù)比轉(zhuǎn)換對平滑效應(yīng)的影響

對數(shù)比轉(zhuǎn)換是處理成分?jǐn)?shù)據(jù)常用的方式。由于成分?jǐn)?shù)據(jù)具有定和、非負(fù)的基本特點(diǎn)，在插值中為了保持這種性質(zhì)，需要使用對數(shù)比轉(zhuǎn)換公式。此外，轉(zhuǎn)換前后數(shù)據(jù)的分布和數(shù)據(jù)維度會有所改變。根據(jù)前人的研究[9,24]，對數(shù)比轉(zhuǎn)換一般能夠改善原始數(shù)據(jù)的正態(tài)性，圖3也說明了這一結(jié)論。從圖3可以看到，除個(gè)別極端值之外的數(shù)據(jù)，正態(tài)性確實(shí)得到了改善，但是同時(shí)也說明了對數(shù)比轉(zhuǎn)換不能很好處理存在極端值的情況。

對數(shù)比轉(zhuǎn)換還會影響插值結(jié)果的平滑性。根據(jù)圖4，經(jīng)過對數(shù)比轉(zhuǎn)換后，土壤顆粒含量插值結(jié)果的平滑性大都得到了一定的改善(黏粒的極差平滑率除外)?？赡苁且?yàn)檗D(zhuǎn)換后的數(shù)據(jù)具有更好的正態(tài)性，使得插值結(jié)果的統(tǒng)計(jì)特征與原始數(shù)據(jù)更加一致，極差平滑率和標(biāo)準(zhǔn)差平滑率也更小。同時(shí)不同的對數(shù)比轉(zhuǎn)換對平滑性影響的程度也不同，但差別不大。從本研究可以看出，ALR和CLR表現(xiàn)相似，主要是因?yàn)閮烧叩霓D(zhuǎn)換形式類似，SLR也是基于ALR和CLR的原理進(jìn)行改進(jìn)，表現(xiàn)卻更差一些，但是ILR對平滑性的改善效果非常明顯，可能是因?yàn)镮LR轉(zhuǎn)換方法符合Aitchison空間中同度量的轉(zhuǎn)換，是直接對數(shù)據(jù)處理的方法[6]。

3.3 ILR-BP模型校正的作用及與其他校正方法對比

克里金插值的平滑效應(yīng)主要原因包括：首先是克里金插值的本質(zhì)是局部線性加權(quán)平均，插值結(jié)果會更靠近局部平均值，因此使用的鄰近點(diǎn)個(gè)數(shù)會影響插值結(jié)果。鄰近點(diǎn)個(gè)數(shù)越多，結(jié)果就會更加穩(wěn)定地趨向平均值，平滑效應(yīng)就越明顯，所以不宜使用過多鄰近點(diǎn)，但是鄰近點(diǎn)又不能太少，否則插值結(jié)果不能反映空間變異規(guī)律。其次與地統(tǒng)計(jì)理論有關(guān)，克里金插值的算法基礎(chǔ)是變異函數(shù)，它能反映區(qū)域化變量空間變異的規(guī)律，但在實(shí)際計(jì)算中，是通過擬合經(jīng)驗(yàn)變異函數(shù)散點(diǎn)得到理論的變異函數(shù)模型，而且模型往往是非線性的。在變異函數(shù)擬合過程中，會抵消原始數(shù)據(jù)中的一部分變異性和隨機(jī)性，使得插值系數(shù)會具有一定的平滑性，最終導(dǎo)致插值結(jié)果的平滑性。但總體上還是遵循地理學(xué)第一定律，即越近的點(diǎn)對插值結(jié)果影響越大。事實(shí)上，平滑效應(yīng)是克里金法的一種固有特點(diǎn)，李超[25]研究指出，應(yīng)當(dāng)客觀看待平滑性，要針對具體的研究目的進(jìn)行分析。

由此可見，鄰近點(diǎn)值的大小和距離插值點(diǎn)的遠(yuǎn)近是影響插值結(jié)果平滑性的主要因素。所以本研究提出使用按一定空間距離排列的鄰近點(diǎn)序列作為輸入特征，通過BP神經(jīng)網(wǎng)絡(luò)模型訓(xùn)練，建立鄰近點(diǎn)特征與插值結(jié)果殘差的關(guān)系，并結(jié)合對數(shù)比轉(zhuǎn)換處理，建立了ILR-BP校正模型。ILR轉(zhuǎn)換主要是為了預(yù)處理成分?jǐn)?shù)據(jù)，在插值后再使用轉(zhuǎn)回公式，這個(gè)非線性過程改變了數(shù)據(jù)的統(tǒng)計(jì)特征，導(dǎo)致極端值的插值預(yù)測比原始數(shù)據(jù)直接插值好，一定程度上削弱了平滑性。進(jìn)一步削弱平滑性主要依靠BP神經(jīng)網(wǎng)絡(luò)模型的后處理方法，由于其直接基于殘差規(guī)律修正插值結(jié)果，大部分樣本點(diǎn)插值誤差都會不同程度減小，尤其是極端值插值結(jié)果得到明顯改善，因此平滑性得到了明顯削弱，同時(shí)插值精度也有所提升。

ILR-BP模型有以下特點(diǎn)：首先，理論簡單，不需要復(fù)雜的統(tǒng)計(jì)學(xué)理論，只需要找到合適的特征值預(yù)測殘差的分布即可。而YAMAMOTO[15]提出的方法是基于標(biāo)準(zhǔn)差數(shù)來校正，需要通過一系列的統(tǒng)計(jì)計(jì)算，稍顯復(fù)雜。其次，ILR-BP模型校正效果明顯，在保證插值精度的基礎(chǔ)上，能在很大程度上改善插值的平滑性。YAMAMOTO[15]校正方法中還需要使用克里金法去插值標(biāo)準(zhǔn)差數(shù)，從本質(zhì)上來講還是局限在克里金法之中；段建軍等[17]使用一元二次回歸模型來模擬殘差的規(guī)律，雖然較為簡單，但是特征點(diǎn)較少，模型不一定能反映真實(shí)的規(guī)律；而本文的ILR-BP模型方法是基于機(jī)器學(xué)習(xí)，擬合能力更強(qiáng)。最后，ILR-BP模型方法也存在一定的局限性，即對采樣點(diǎn)的個(gè)數(shù)有要求，如果采樣點(diǎn)太少，可能無法訓(xùn)練出較穩(wěn)定的模型。

3.4 河套灌區(qū)土質(zhì)空間分布的影響因素分析

河套灌區(qū)是由于黃河沖積作用形成的沖積平原，因此其中的母質(zhì)大部分來自于上游的泥沙，整體上土質(zhì)偏砂性。灌區(qū)西部位于河流上游，由于河流運(yùn)動(dòng)的分選和沉積作用，上游的顆粒偏粗，而中下游的顆粒更細(xì)，因此砂粒主要分布在灌區(qū)西部，而中東部更多的是粉粒和黏粒。但是由于人類活動(dòng)和其他因素的影響，局部也會有異?，F(xiàn)象，比如中部存在砂粒集中的區(qū)域，這在河套灌區(qū)也較為常見，灌區(qū)中零散分布著大小不一的沙丘，這正是黃河沖積作用的見證，后來由于人類耕作活動(dòng)使得土壤得到了一定的改良。

4 結(jié)論

(1)基于對數(shù)比轉(zhuǎn)換的普通克里金法對平滑效應(yīng)有一定的改善作用，不同轉(zhuǎn)換方法的改善程度也不相同。ALR和CLR的效果一般，作用類似，SLR的改善作用最小，ILR的效果最好，改善后的黏粒、粉粒和砂粒的標(biāo)準(zhǔn)差平滑率依次為0.22、0.20和0.15。

(2)基于ILR轉(zhuǎn)換和BP神經(jīng)網(wǎng)絡(luò)相結(jié)合的校正模型(ILR-BP)能夠較好地削弱克里金插值中的平滑效應(yīng)。校正后的土壤組分交叉驗(yàn)證結(jié)果的極差平滑率為0，標(biāo)準(zhǔn)差平滑率為0.03～0.07，同時(shí)平均絕對誤差和均方根誤差都有一定的下降，這說明平滑性得到校正的同時(shí)，插值精度也有所提高。

(3)校正后的土壤顆粒組成和質(zhì)地類型的空間分布更加精細(xì)，局部特征更加明顯，并且能夠很好地還原樣本的特征，在極端值的表現(xiàn)上更加符合實(shí)際。河套灌區(qū)中部和西部主要是偏砂性土壤，壤土主要集中分布在中西部、東北部和東南部，粉壤土主要分布在北部和南部接近黃河的地區(qū)。

(4)基于ILR-BP的普通克里金法不但能夠解決成分?jǐn)?shù)據(jù)空間插值問題，同時(shí)能夠削弱插值的平滑效應(yīng)，得到較為真實(shí)的區(qū)域土壤質(zhì)地分布，為河套灌區(qū)的水鹽運(yùn)移研究和水土資源規(guī)劃提供了基礎(chǔ)資料。