王銘鑫,范 超,高秉博**,任周鵬,李發(fā)東

(1.中國農(nóng)業(yè)大學土地科學與技術(shù)學院 北京 100193;2.武漢大學遙感信息工程學院 武漢 430079;3.中國科學院地理科學與資源研究所 北京 100011)

土壤是農(nóng)業(yè)生產(chǎn)的基礎,土壤環(huán)境屬性影響耕地的生產(chǎn)能力與糧食安全。土壤受到重金屬污染后不僅會影響作物的產(chǎn)量和農(nóng)產(chǎn)品的安全,危害人類健康,還會進入水體和空氣,造成次生生態(tài)環(huán)境問題。我國由于近些年工業(yè)迅猛發(fā)展、農(nóng)業(yè)化學水平提高、城鎮(zhèn)快速擴張以及對環(huán)境保護的忽視,使得大量的有毒有害物質(zhì)進入土壤,造成了嚴峻的污染形勢。精確的環(huán)境污染空間分布圖是開展各項污染風險管控和防治工作的基礎。由于目前土壤環(huán)境污染以實地采樣調(diào)查為主要手段,成本較高,樣點分布相對稀疏。因此,需要采用空間插值方法基于樣本數(shù)據(jù)完成區(qū)域土壤環(huán)境屬性分布圖。然而,由于土壤污染具有較強的空間異質(zhì)性,為空間插值精度的提升帶來較大的挑戰(zhàn)。

目前,土壤環(huán)境屬性空間插值方法主要分為基于空間位置關(guān)系信息的插值、基于屬性相似信息的插值以及融合空間位置關(guān)系信息和屬性相似信息的插值方法?；诳臻g位置關(guān)系信息的插值以描述空間自相關(guān)性的地理學第一定律和描述空間異質(zhì)性的第二定律為基礎。基于地理學第一定律包括最鄰近點插值、反距離加權(quán)、徑向基函數(shù)和局部多項式等確定性插值方法,以及普通克里金、簡單克里金和指示克里金等地統(tǒng)計方法。與確定性插值方法相比,地統(tǒng)計方法能在給出無偏最優(yōu)插值結(jié)果的同時,得到對插值結(jié)果的不確定性估計,因此應用較為廣泛,但是難以解決土壤環(huán)境屬性的空間異質(zhì)性問題。各向異性半變異函數(shù)方法、三明治方法、空間分層異質(zhì)性方法等空間插值方法在考慮空間自相關(guān)性的基礎上,考慮到對空間異質(zhì)性的處理,能一定程度上解決土壤環(huán)境屬性空間異質(zhì)性插值問題。

基于屬性相似信息的插值方法以地理學第三定律為基礎。地理學第三定律認為:具有相似地理環(huán)境的兩個點(區(qū)域),其目標變量具有相似(近)的取值。包括相似度度量法、Soil-Land Inference Model(SoLIM,土壤-環(huán)境推理模型)、回歸方法和可以處理多變量的嶺回歸、lasso 回歸、隨機森林、提升樹和支持向量機等監(jiān)督學習方法?；趯傩韵嗨菩畔⒌姆椒梢猿浞掷猛寥牢廴緛碓醇捌溆绊懸蛩氐容o助變量數(shù)據(jù),改善土壤環(huán)境屬性的插值精度。但是由于土壤環(huán)境污染源及其影響因素數(shù)據(jù)難以全部精確獲得,同時為了防止模型過擬合,因此模型會產(chǎn)生殘差。由于污染物的擴散遷移特性,殘差往往具有空間自相關(guān)性。因此在考慮屬性相似性的基礎上,綜合應用空間位置關(guān)系信息能夠提高空間插值精度。

融合空間位置關(guān)系信息和屬性相似信息的插值方法包括空間滯后模型、空間誤差模型、空間變系數(shù)回歸和地理加權(quán)回歸等空間回歸方法,以及協(xié)克里金、回歸克里金和因子克里金等地統(tǒng)計方法。其中回歸克里金是數(shù)字土壤制圖領(lǐng)域應用最為廣泛的一種插值方法。但是當輔助變量較多且包含較多類型變量時,回歸克里金處理多變量的能力不足,難以有效利用土壤環(huán)境屬性空間插值所涉及的多維輔助變量信息。

機器學習方法在處理高維輔助變量方面具有獨特的優(yōu)勢,但是主流的機器學習方法僅能利用多維輔助變量信息,并沒有融合空間位置關(guān)系信息進一步改進插值精度。如何將空間位置關(guān)系信息融入機器學習方法改進土壤環(huán)境空間插值精度成為目前數(shù)字土壤制圖領(lǐng)域最前沿的研究方向之一。Sekuli?等提出了隨機森林空間插值法(random forest spatial interpolation,RFSI)。該方法將鄰近若干樣點的空間橫、縱坐標和目標變量作為單獨的屬性加入隨機森林,以改進空間插值結(jié)果。但是該方法分離了樣點位置與目標變量之間的關(guān)系,未能融合空間位置關(guān)系信息和屬性相似信息。Hengl 等提出的隨機森林空間預測框架(random forest for spatial predictions framework,RFsp),將到樣點的距離作為預測屬性加入隨機森林方法,將位置關(guān)系信息有效反映在預測屬性中,完成了空間位置信息和屬性相似信息的融合,可用于改善土壤環(huán)境屬性插值精度。但是土壤環(huán)境變量的空間變異并不一定與距離呈現(xiàn)線性關(guān)系。因此,為了進一步改進土壤環(huán)境屬性插值精度,本文在隨機森林空間預測框架的基礎上,將空間距離替換為空間半變異值,融合隨機森林與地統(tǒng)計中的空間半變異函數(shù),提出了融合半變異函數(shù)的空間隨機森林插值方法(spatial random forest with semi-variogram interpolation,SRFsei),以充分融合空間位置信息和屬性相似信息,提高土壤環(huán)境屬性的空間插值精度。并以湖南省湘潭縣北部部分鄉(xiāng)鎮(zhèn)的土壤重金屬數(shù)據(jù)作為案例,對比分析了該方法的空間插值效果。

1 面向土壤環(huán)境插值的空間隨機森林方法

1.1 隨機森林

隨機森林是一種基于分類回歸樹的集成學習方法,由多個具有高預測精度和低相關(guān)性的分類回歸樹組成。該方法通過隨機選取樣點和特征的重采樣訓練多個分類回歸樹,在預測時綜合多個分類回歸樹的結(jié)果,提高預測精度并解決過擬合問題。其基本流程如下:

1)從原始訓練集中采用bootstrap 抽樣抽取個樣本和個屬性；

2)將個樣本分別建立個分類回歸樹,完成模型訓練；

3)對于待預測點,將其帶入個分類回歸樹,得到每個分類回歸樹的預測結(jié)果；

4)綜合個分類回歸樹的預測結(jié)果,給出對待預測點的最終預測結(jié)果。

隨機森林(RF)是一種集成學習方法,最終預測結(jié)果一般是個分類回歸樹結(jié)果的平均值,如公式(1):

1.2 隨機森林空間預測框架

RF 在建模的過程中忽略待估計點與樣本點之間的空間自相關(guān)性,可能會導致預測結(jié)果不準確,出現(xiàn)過高估計和過低估計,尤其是在目標變量的具有較強的空間自相關(guān)性時。為了彌補這項不足,Hengl 等提出了隨機森林空間預測框架(RFsp):

式中:()為空間點上的目標變量取值;為輔助變量;為遙感圖像的光譜信息;為空間距離屬性,即距離空間樣點的距離。RFsp 將與樣點的距離作為輔助變量加入預測變量,并利用隨機森林方法處理多維預測變量,實現(xiàn)了輔助變量信息與位置關(guān)系信息的融合。Hengl 等已經(jīng)證明,隨機森林與線性地統(tǒng)計建模和克里金等技術(shù)相比,其優(yōu)點在于不需要遵循平穩(wěn)性假設,且當樣點具有代表性,而目標變量、協(xié)變量和空間依賴結(jié)構(gòu)之間的關(guān)系是復雜的、非線性的情況下,可以改進空間插值結(jié)果。

1.3 融合半變異函數(shù)的空間隨機森林方法

隨機森林空間預測框架(RFsp)采用空間距離關(guān)系反映樣點之間的空間位置信息。但是由于模型構(gòu)建的最終目標是實現(xiàn)土壤環(huán)境屬性的精確插值,而空間位置信息并未直接與土壤環(huán)境屬性的空間變異關(guān)聯(lián),同時土壤環(huán)境變量的空間變異并不一定和距離呈線性關(guān)系,因此本文提出結(jié)合地統(tǒng)計空間半變異函數(shù)的改進方法,即融合半變異函數(shù)插值的空間隨機森林方法(spatial random forest with semi-variogram interpolation,SRFsei)。該方法的定義如公式(3):

式中:是基于空間位置關(guān)系和空間半變異函數(shù)所得的變量組,其中每一個變量的取值為待估計點到單個樣本點之間半方差的平方根,定義如公式(4):

式中:h為待估計點到樣點的距離,γ()為空間半變異函數(shù),定義如公式(5):

式中:()為一區(qū)域化隨機變量,為兩空間點之間的距離,(s)和(s+)分別是區(qū)域化變量()在空間位置s和s+處的值,()表示在空間上所有距離為的離散點對的數(shù)量。通過樣本點計算半方差后,可以采用球狀模型、高斯模型和指數(shù)模型等擬合,獲得空間半變異模型的函數(shù)形式。空間半變異函數(shù)由模型類型、基臺值(Sill)、變程(Range)和塊金值(Nugget)唯一確定?！盎_值”是半變異函數(shù)隨著距離增加達到一個相對穩(wěn)定的常數(shù)的值;“塊金值”指的是間隔距離=0 時的半變異函數(shù)值,塊金值是由于混合樣的存在和測量誤差造成的;“變程”指的是變異函數(shù)達到基臺值時的間隔距離,當兩點間的間距大于等于變程(range)時,區(qū)域化變量(s)的空間相關(guān)性消失。

2 案例研究

2.1 研究區(qū)描述

研究區(qū)為湖南省湘潭縣北部5 個鄉(xiāng)鎮(zhèn)(112°57′E,27°46′N),如圖1所示,總面積為522.46 km。研究區(qū)位于湘江東岸,主要地貌類型為低山、丘陵和平原,屬亞熱帶季風濕潤氣候。2011年共采集688 個農(nóng)田土壤樣本點,并使用火焰原子吸收光譜法分析了土壤的Cr 含量。本研究以土壤Cr 含量為目標變量。

圖1 研究區(qū)域及樣點分布Fig.1 Study area and sampling point

目前已有多項研究證實地形(例如海拔高度、坡度等)、土壤類型、土地利用類型是Cr 分布的重要影響因子。海拔高度通過土壤侵蝕來影響重金屬的遷移,土壤類型和土地利用類型影響重金屬的累積狀況。因此,本文選取土壤類型、土地利用類型、海拔高度、坡度等因素作為土壤重金屬含量插值研究的輔助變量。各因子的空間分布情況如圖2所示。

圖2 研究區(qū)輔助變量土壤類型(a)、土地利用類型(b)、海拔(c)和坡度(d)的映射Fig.2 Mapping of auxiliary variables in the study area.a:soil type;b:land use type;c:altitude;d:slope.

2.2 試驗方法

對于研究區(qū)采樣數(shù)據(jù),分別采用本文提出的SRFsei 方法、RF、RFsp、普通克里金(Ordinary Kriging,OK)和泛克里金(Universal Kriging,UK)方法進行插值,通過交叉檢驗對比不同方法的插值精度。所有對比試驗均采用R 語言實現(xiàn),空間半變異函數(shù)、OK和UK 使用gstat 實現(xiàn),隨機森林訓練和預測基于ranger 包實現(xiàn)。為了比較不同插值方法的結(jié)果,采用平均絕對誤差(MAE)和均方根誤差(RMSE)表征插值精度,兩個指標可以避免正負偏差互相抵消的情況,用來量化預測結(jié)果的平均誤差。MAE 和RMSE的計算方法如下:

為了對比SRFsei 方法相比其他插值方差的插值精度改進程度,本文定義了改進程度指數(shù),如公式(8):

式中:imp為SRFsei 方法相比第種方法的插值精度改進程度指數(shù);為SRFsei 方法的插值誤差,可以是MAE或RMSE;E為作為對比的第種方法的插值誤差,可以采用MAE或RMSE,但需要與 E對應一致。

處理覆蓋整個研究區(qū)域的輔助變量數(shù)據(jù),在此基礎上,使用全部采樣點,采用SRFsei 方法、RF、RFsp、普通克里金(Ordinary Kriging,OK)和泛克里金(Universal Kriging,UK)方法對整個研究區(qū)域的Cr 含量進行插值制圖,獲得插值制圖結(jié)果和不確定性估計結(jié)果。

2.3 結(jié)果與分析

原始樣本的Cr 含量的變異函數(shù)圖如圖3所示,利用研究區(qū)域數(shù)據(jù)和R 語言工具繪制出圖。使用如公式(9)所示的球狀模型擬合:

圖3 研究區(qū)Cr 含量半變異函數(shù)Fig.3 Variation function of Cr content in the study area

式中:為塊金效應常數(shù),為拱高,+為基臺值,為空間依賴范圍,即變程值。

其中,研究區(qū)域內(nèi)土壤中Cr 含量的變差函數(shù)的變程值為6421.279,塊金值為77.3,基臺值為369.7。塊金值與基臺值之比(基底比)是表示空間自相關(guān)性程度的指標,一般認為該值為25%和75%是區(qū)分相關(guān)性程度等級的兩個分界值,如果該比值小于25%,則空間相關(guān)性程度很高;如果比值在25%～75%,則空間相關(guān)性程度為中等;如果比值大于75%,則由隨機變化引起的空間變異性程度較大。本文土壤中Cr含量的變差函數(shù)的塊金值與基臺值比率為20.9%,說明研究區(qū)域內(nèi)土壤中Cr 含量存在著較強的空間自相關(guān)性。

本文提出的SRFsei 與RF、RFsp、OK 和UK 的插值交叉檢驗結(jié)果如圖4所示。由圖4 中MAE 結(jié)果可知,SRFsei 具有最小的空間插值誤差,OK 具有最大的空間插值誤差,UK 由于采用了輔助變量,插值精度較OK 有所改善。RF 和RFsp 比OK 和UK有更好的插值精度,在RF 和RFsp 中,RFsp 由于考慮了空間位置關(guān)系,插值精度更高。圖4 中RMSE結(jié)果與MAE 基本一致,SRFsei 具有最優(yōu)的插值精度,OK 具有最大的RMSE。與MAE 不同之處是RFsp方法的RMSE 大于RF 方法。

圖4 不同插值方法交叉驗證結(jié)果的平均絕對誤差(MAE)和均方根誤差(RMSE)Fig.4 Mean absolute error(MAE)and root mean square error(RMSE)of cross-validation results of different interpolation methods

SRFsei 相對傳統(tǒng)的OK 和UK 方法在MAE 和RMSE 上的改進均大于10%。相比RF,SRFsei 在MAE上的改進為6.52%,在RMSE 上的改進為5.34%。而相比RFsp,在MAE 上的改進為5.97%,在RMSE 上的改進為6.49%。從交叉檢驗結(jié)果可知,SRFsei 相比傳統(tǒng)的地統(tǒng)計方法有較大的插值精度提升,而相對于新型的機器學習方法精度提升有限。

采用全部688 個樣點的數(shù)據(jù)基于SRFsei 方法對研究區(qū)空間插值制圖如圖5所示。其中圖5a 為插值制圖結(jié)果,圖5b 為插值不確定性結(jié)果。由圖5a 可知,研究區(qū)域中東部污染較為嚴重,西北部污染較輕。由圖5b 可知,研究區(qū)域絕大部分的插值誤差標準差小于5,誤差較大的區(qū)域主要集中在污染較嚴重的中部區(qū)域。圖中黑色方框所示部分與白色方框所示部分具有相同的污染濃度,但是從誤差圖中可以得出,黑色方框所示部分具有較小的插值誤差,因此可以確定該部分污染較為嚴重,而白色方框部分誤差較高,因此該部分的插值結(jié)果具有較大的不確定性。如需更精確地掌握污染情況,還需要進一步展開加密采樣調(diào)查。

圖5 融合半變異函數(shù)的空間隨機森林插值法的插值結(jié)果圖(a)與誤差標準差圖(b)Fig.5 Interpolation results(a)and standard deviation variance of error(b)of spatial random forest with semivariogram interpolation

RF、RFsp、UK 和OK 對研究區(qū)域的插值制圖結(jié)果及其不確定性估計圖如圖6所示。土壤Cr 含量插值結(jié)果總體相似,但在局部上差異明顯:其中OK和UK 的插值結(jié)果,東南部區(qū)域由于協(xié)方差矩陣是奇異矩陣,沒有計算出插值結(jié)果。相比較OK 方法,其他方法由于使用了輔助數(shù)據(jù),因此具有更豐富的細節(jié)信息。而RF 方法,由于沒有利用空間自相關(guān)性,插值結(jié)果存在不自然的非平滑變化。RFsp 插值結(jié)果中,存在條帶化現(xiàn)象。相比其他方法,SRFsei 插值結(jié)果包含了更豐富的細節(jié)信息,空間變化也更加合理。

圖6 不同插值方法的插值結(jié)果圖(a)及誤差標準差圖(b)Fig.6 Interpolation results(a)and standard deviation variance of error(b)of different interpolation methods

3 結(jié)論

本文基于隨機森林空間預測框架,通過結(jié)合隨機森林方法與空間變異函數(shù),提出了融合半變異函數(shù)的空間隨機森林方法(SRFsei)。該方法在空間插值建模和預測中能夠?qū)⒎从惩寥拉h(huán)境污染源及影響因素的多維輔助變量信息與空間變異信息有效結(jié)合,同時考慮屬性相似性信息和空間位置關(guān)系信息,因而可以有效提高插值進度,得到更為可靠的插值結(jié)果。湖南省湘潭縣北部鄉(xiāng)鎮(zhèn)土壤重金屬數(shù)據(jù)插值案例結(jié)果表明,本文所提出的方法相較于其他傳統(tǒng)地統(tǒng)計的OK 和UK 方法插值精度提升10%以上,相較于新型的機器學習空間插值方法的插值精度提升5%以上。同時,其插值制圖結(jié)果具有更豐富的細節(jié)信息和合理空間變化。因此,本文所提出方法較適用于具有復雜空間異質(zhì)性的土壤環(huán)境變量下的空間插值。但是其插值效率還有待進一步在其他區(qū)域應用中檢驗。