基于BME-GWR 法的景觀單元土壤有機碳密度空間預測

劉歡瑤，孟 岑，鄒冬生，吳金水

（1.湖南農業大學資源環境學院，湖南 長沙 410128；2.中國科學院亞熱帶生態農業研究所，湖南 長沙 410125）

【研究意義】土壤是陸地生物圈中最大的碳庫，土壤有機碳（SOC）控制全球碳循環，是重要的大氣碳源［1］。作為評價土壤質量的重要指標，土壤有機碳庫的研究對防治土壤退化和發展可持續農業具有重要意義。準確估算、分析土壤有機碳的儲量和空間分布，為有效調控碳源/匯方向、全球氣候變化以及為土壤質量綜合評估提供了理論支持。土壤有機碳密度（SOCD）的估算精確度受樣點數量的影響。隨著3S技術的飛速發展，與SOCD相關的環境變量數據獲取變得簡單、便捷，充分利用多源環境先驗信息作為“軟數據”是克服傳統以克里格為代表的地統計方法缺陷，提高SOCD估算精度的有效技術途徑。

【前人研究進展】當前已有很多插值方法結合環境輔助變量進行土壤屬性的空間預測，包括多元線性回歸模型、協同克里金、回歸克里金、地理加權回歸等，當環境因子與土壤屬性相關性較強時這些方法預測精度均高于不結合環境輔助變量的地統計方法［2］。然而，由于多數插值方法均未考慮環境輔助變量的先驗分布，很難處理非高斯分布變量，因而限制了利用多源輔助數據進行空間預測的能力［3］。貝葉斯最大熵（Bayesian Maximum Entropy,BME）結合了貝葉斯方法論和最大熵原理來模擬時空變量，將精度較高的“硬數據”和多渠道源得到的“軟數據”區分使用，能夠有效地綜合利用各種不同來源和精度的數據［4］。近幾年來，BME模型嘗試被引入到土壤、大氣污染和環境的時空預測中［4-6］。例如，費徐峰等［5］將田間測量的土壤重金屬含量數據作為軟數據，將實驗室測定的相應土壤重金屬含量作為硬數據，利用土壤重金屬在土壤母質類型中的概率密度函數作為軟數據，建立了土壤重金屬含量BME和普通克里格（OK）插值模型，以探討研究區內重金屬污染分布情況和影響因素。

【本研究切入點】軟數據本身的來源多樣化，例如將另一種檢測手段得到的數據集作為軟數據，按環境相關法得到預測屬性的概率分布作為軟數據［7］，用空間預測模型結合輔助變量獲得相對不確定性的數據作為軟數據［4］，但較少有文獻報道軟數據的獲取途徑對BME結果的影響。受復雜的成土因素（如植被、地形）等環境因子影響，隨著研究尺度的變化，土壤屬性的主控因素也會隨之變化［8］。因此，選擇合適范圍建立環境輔助數據和預測變量之間的關系，是有效控制預測模型不確定性產生的途徑。【擬解決的關鍵問題】本研究基于BME方法框架，選取湖南省桃源縣盤塘鎮王家垱村內SOCD實測數據作為“硬數據”，以地理加權回歸模型（GWR）結合地形、土地利用等多源環境數據獲得“軟數據”，比較貝葉斯最大熵結合地理加權回歸模型（BMEGWR）與傳統GWR方法對土壤SOCD空間分布預測精度，進一步分析基于土地利用類型估算得到的“軟數據”所計算的BME-GWR模型模擬精度的變化，為綜合利用環境變量提高土壤屬性的空間預測模型精度提供依據。

1 材料與方法

1.1 研究區概況

綜合考慮地形地貌、土地利用、土壤等因素，在我國亞熱帶紅壤丘陵區選取具有代表性的農業生態景觀單元作為研究區域［9,11］。研究區域位于湖南省桃源縣盤塘鎮王家垱村，總面積為355.12 hm2，地處中亞熱帶北緣，為季風性濕潤氣候（年平均降雨量為1330 mm，年平均氣溫為16.8 ℃）。該區域屬于低山崗地地貌，高程約為81～112 m，土壤類型以地帶性的紅壤為主。研究區土地利用方式主要為林地（35.7%）、稻田（39.9%）、果園（15.1%）、旱地（9.3%），其中林地以人工林（馬尾松、櫟樹、樟樹）為主，稻田以雙季稻為主，旱地主要作物為棉花、油菜、玉米、苧麻等經濟作物，果園以橘子、柚子、桃子為主［10］。

1.2 研究方法

1.2.1 數據來源與處理 通過GPS準確定位采樣點，并記錄采樣點的經緯度、高程，通過ArcGIS 10.3軟件經投影轉化，得到以米為單位的平面坐標；再將矢量等高線和高程控制點轉化為分辨率為5 m的數字高程模型（圖1），并提取1:10000航拍的土地利用信息，得到分辨率為5 m土地利用類型分布圖，在數字高程模型中提取坡度、地形濕度指數等數據。于2003年4—5月，以樣區中心為橫縱坐標軸的交叉點，劃定1.0 cm×1.0 cm的網格，按采樣密度耕地（稻田和旱地）5個土樣/ hm2、果園2～3個土樣/ hm2、林地1個土樣/ hm2的原則隨機布點，并采用“S”形路線多點取樣法（不少于15點）采集表層土樣（0～20 cm），共采集土壤樣品523個［11］，其中，稻田、旱地、果園、林地樣品數分別為249、66、192、16個。

土壤容重采用環刀法測定，土壤有機碳含量采用碳氮元素分析儀（Vario-MAX C/N，德國）測定。SOCD指單位面積內一定厚度土層的土壤有機碳儲量（kg/m2），計算公式為：

式中，SOCDi為第i層土壤有機碳密度（kg/m2）；SOCi為第i層土壤有機碳含量（g/kg）；ρi為第i層土壤容重（g/c m3）；Di為土壤厚度（cm），本研究取值20 cm。

圖1 研究區樣點分布、高程、坡度、地形濕度指數、土地利用類型Fig.1 Sample plot distribution,elevation,slope,topographic wetness index（TWI）and land use types in the research area

1.2.2 地理加權回歸（GWR）GWR模型是普通線性回歸的拓展，可被視作局部的加權最小二乘回歸模型，允許局部參數估計，將觀測數據空間特性的變化參數嵌入到模型中，并假定線性回歸模型中的回歸參數是觀測點處的任意函數，得到GWR模型，反映樣本對回歸方程貢獻在空間上的分異［13］，其計算公式為：

式中，y（u）為位置u的因變量值；β0（u）為截距值；βk（u）為第K個協變量在位置u上的回歸系數，回歸系數由每個觀測點周圍數據估算而得；xk（u）為在位置u上的第k個協變量的值；ε（u）為在位置u上的隨機誤差值。GWR回歸模型中的參數采用加權最小二乘法計算，本研究把SOCD測定數據作為因變量，高程、坡度、地形濕度指數、土地利用類型數據作為自變量進行GWR模型擬合。

1.2.3 貝葉斯最大熵模型（BME）BME模型結合了最大熵原理和廣義貝葉斯條件公式，可以利用不同精度、質量、表示方式（如區間、分類等）的多源數據，這些數據總體上可以分為廣義知識（General knowledge，G）和特定知識（Special knowledge，S）兩類。廣義知識用來描述自然規律、物理法則和硬數據的統計規律（如數學期望、半方差）等。特定知識包括硬數據和軟數據，硬數據是特定觀測點上的誤差可以忽略的數據，如實測數據、時間序列數據等，本研究以樣點SOCD實測數據作為硬數據；而軟數據是相對模糊的、有誤差的數據，包括歷史數據、粗測量數據、模型模擬數據等。

BME的基本原理是利用廣義知識（G）通過最大熵原理獲得先驗概率密度函數pdf〔fG（ymap）〕，再基于結合了硬數據和軟數據的特定知識（S）、貝葉斯條件概率將研究區未測點變量分布的先驗pdf轉換為后驗pdf〔fT（ymap）〕，并可計算待預測點xk的預測均值或最大值。

式中，I=Imn+1U Imn+2U…Im是軟數據集，ydata=[y1，y2，…ym]是特定知識。為預測結果的平均值；為預測結果的最大值；yk為測量值；f（yk）為后驗pdf；dyk為測量值的變化量。

本研究地理加權回歸結合貝葉斯最大熵模型（GWR-BME）是以GWR模型所預測連續柵格表面為基礎，將在樣點空間位置上提取的被預測的SOCD屬性值作為軟數據；而按土地利用類型估算的貝葉斯最大熵模型（GWR-BMEL），則是先提取不同的土地利用類型，在每個土地利用類型單獨以SOCD實測數據作為因變量，高程、坡度、地形濕度指數作為自變量進行GWR模型預測，提取樣點位置的GWR模型預測值作為軟數據。最后，結合“軟數據”和“硬數據”進行BME預測，GWR-BMEL則要將各土地利用類型的預測結果合并得到整個區域的SOCD空間分布。

2 結果與分析

2.1 土壤有機碳密度與環境變量描述性特征統計

從表1 可以看出，研究區域的SOCD 分布范圍為1.30～4.71 kg/m2，平均值為3.22 kg/m2，標準差為0.77 kg/m2，變異系數為24%。不同土地利用類型下SOCD差異顯著，表現為稻田〔3.84（±0.44）kg/m2〕>旱地〔2.92（±0.57）kg/m2〕>果園〔（2.58（±0.44）kg/m2〕>林地〔2.31（±0.55）kg/m2〕，變異系數為11%～24%，均屬于中等變異（10%～100%）。

通過數字高程圖提取高程、坡度、地形濕度指數作為地形特征，與SOCD 進行相關分析。從圖2 可以看出，研究區域內提取的高程、坡度、地形濕度指數對SOCD 有明顯的影響，可以作為地形因子進行空間分布擬合。其中，SOCD 與高程和坡度均呈極顯著負相關、相關系數分別為-0.68 和-0.62，而與TWI 呈極顯著正相關、相關系數為0.34。

表1 研究區土壤有機碳統計Table 1 Statistical values of soil organic carbon density（SOCD）in the research area

圖2 研究區土壤有機碳密度、高程、坡度、地形濕度指數相關分析Fig.2 Correlation analysis of soil organic carbon density（SOCD）,elevation,slope and topographic wetness index（TWI）in the research area

2.2 半方差函數比較

根據最大決定系數（R2）原則，GWR、GWR-BME、GWR-BMEL3種模型最優擬合模型均為指數模型（圖3）。半方差函數模型中所得的主要參數為塊金值（C0）、基臺值（Sill）、變程（Range）以及塊基比值（C0/Sill），可以反映一定空間范圍內的空間變異性和變量的自相關性。其中，塊金值反映了隨機部分形成變量的變異性和測量誤差，GWR、GWR-BME、GWR-BMEL3種模型塊金值分別為0.32、0.24、0.16 kg/m2，結合BME方法可以減少由隨機部分形成的變異；基臺值是變異函數數值隨著采樣點間距的增大從塊金值增大到最大的常數值，反映變量的最大變異性，體現了非隨機部分形成的變異，GWR、GWR-BME、GWR-BMEL的基臺值分別為0.61、0.66、0.52 kg/m2；塊金值與基臺值的比值分別為52%、36%、31%，為中等程度的空間依賴性（25%～75%），尤其是GWR-BMEL表現出的空間自相關性最強；變程從大到小依次為GWR（269 m）、GWR-BME（207 m）、GWRBMEL（138 m），GWR-BME、GWR-BMEL最優擬合模型變程均小于GWR，進一步驗證了塊金值與基臺值比值的結果，反映出較強的空間自相關性。

2.3 預測精度檢驗

圖3 半方差模型對比Fig.3 Comparison of semivariance model

其中，MAE從大到小依次為GWR（0.25）>GWR-BME（0.21）>GWR-BMEL（0.19），RMSE從大到小依次為GWR（0.40）>GWR-BME（0.35）>GWR-BMEL（0.33）。GWR-BMEL模型的MAE和RMSE比GWR-BME模型分別減少9.53%、5.71%，比GWR模型分別減少24%、17.50%。R2以GWR模型最小（0.72），GWRBME（0.79）比GWR提高9.72%，GWR-BMEL的R2最大、為0.81。表明在軟數據的輔助下，GWRBME和 GWR-BMEL的模型精度均高于GWR模型，而GWR-BMEL模型在3種方法中插值精度最高，GWR-BME次之，GWR的模型精度最低（表2）。說明結合了土地利用類型劃分估算單元所得到的軟數據可以較好地體現土地利用類型內部的空間差異性，克服整體估算結果的平滑效應，提高了估算結果的精度。

表2 各土壤有機碳密度模型精度比較Table 2 Comparison of accuracy among different soil organic carbon density（SOCD）models

2.4 土壤有機碳密度空間分布

圖4 土壤有機碳密度（SOCD）空間分布Fig.4 Spatial distribution of soil organic carbon density（SOCD）

利用GWR-BMEL方法預測研究區域的土壤SOCD的空間分布（圖4），得到分辨率為5 m的空間預測結果。從圖4可以看出，GWR-BMEL法擬合的SOCD空間分布特征的整體趨勢特征明顯，反映細節的能力強，對實測點鄰域范圍的空間結構體現較合理。農業生態景觀單元內SOCD的預測結果具有明顯的區域分布特征，SOCD的空間分布格局與土地利用、地形因子的相關性明顯：在研究區域中部高程和坡度較小的坡度、坡下的稻田以及稻田周圍的旱地的SOCD較高；而研究區域東部高程和坡度較大的坡腰、坡頂區域，主要分布林地和果園的區域是SOCD的低值區。

3 討論

我國全國尺度土壤有機碳密度的平均值為3.35 kg/m2［14］，東部地區表層土壤的SOCD為6.8～21.4 kg/m2［15-16］，其中浙江和江西兩省典型森林類型土壤平均SOCD為7.12～15.69 kg/m2［17］。而本研究區域的土壤有機碳密度比其他研究結果較低，這可能是由于農業生態景觀單元地表的作物成熟后多被收割，直接歸還到土壤部分的凋落物相對減少導致土壤有機碳含量較低［18］，且研究區域的林地以人工經濟林為主，林分結構單一，植被凈生產力小于天然林［16］，同時人工林受人為擾動較大，土壤微生物的分解作用加強，限制了當地土壤有機質的積累［19］。另外，研究區以低山崗地地貌為主，耕地多位于低山的坡中和坡下位置，水土流失較嚴重，導致土壤相對貧瘠［20］。

本研究區SOCD的GWR方法最優擬合模型為指數模型，這與很多相關的土壤有機碳研究結果一致［21-22］。研究區土壤有機碳密度的塊金值和基臺值分別為0.32、0.61 kg/m2，塊金值與基臺值的比值為52%。塊金值與基臺值的比值受結構和隨機因素的影響，較高的比值說明空間變化是由施肥、耕作、其他人類活動等隨機因素引起，而較低的比值說明由采樣和實驗分析本身的系統誤差、土壤性質、地形、其他自然因素等結構因素在空間變化中起重要作用。一般用<0.25、0.25～0.75、>0.75分別表示土壤屬性的強、中等、弱空間自相關性，本研究區的土壤有機碳密度為中等程度空間變異［23］。土壤有機碳密度半方差模型的變程為269 m，是土壤有機碳密度的最大相關距離，表明本研究SOCD采樣的網格間距布設合理，推導出的空間結構關系與地面的空間結構關系一致［24］。與此相應，SOCD的GWR-BME、GWR-BMEL兩種擬合模型的塊金值、塊金值與基臺值的比值、變程均小于GWR，表明BME法能更強地解釋SOCD的空間異質性。同時，GWRBMEL的基臺值最小，體現了同一類型的景觀約束下SOCD的變異性會降低［25］。

結合BME方法用于模擬SOCD空間分布的精度高于GWR，這與多數研究結果相同。楊勇等［7］以DEM生成的相關地形因子作為環境數據，分別用普通克里格（OK）和BME模型對湖北省京山縣土壤有機質含量的密集樣本和稀疏樣本進行預測，結果表明BME的預測精度高于OK，尤其對稀疏樣本預測的精度提高幅度更大。Xiao等［4］采用GWR-BME模型，選取氣溶膠光學厚度、地形數據、氣象數據和污染排放等輔助信息對我國大區域尺度的PM2.5空間分布進行擬合，GWR-BME預測結果精度較GWR模型更高。與傳統方法相比，GWR-BME、GWR-BMEL方法結合了“硬數據”和“軟數據”［6］，降低了預測結果的不確定性，尤其是通過劃分土地利用類型獲取軟數據的GWR-BMEL模型，通過在空間范圍上對每個土地利用類型單獨模擬，獲得其空間變化的局部特征以及空間對象本身相關性和異質性參數，進一步有效利用“軟數據”的空間自相關性，可以較好地體現研究區SOCD的空間分布。

本研究利用GWR-BMEL對研究區SOCD模擬分布，結果表明在中間高程和坡度較小的稻田以及稻田周圍的旱地SOCD較高。這是由于農田（稻田和旱地）的復種指數高，化肥投入量高，導致表層土壤的養分較高，導致SOCD高于林地和果園，尤其是稻田由于淹水環境會降低土壤有機碳的分解速率，進一步降低土壤有機質的分解速率，有利于土壤有機碳積累［25］。同時，研究區域的水稻和旱地作物經常種植在坡度、海拔較低的區域［26］，因而形成當地SOCD的高值區域集中分布在坡度、海拔較低的稻田和旱地的格局。高程和坡度較大的林地和果園是主要的SOCD低值區域，其中研究區域海拔和坡度較大的林地種植年限較短，以種植馬尾松和油茶等人工林為主，其凈生產力遠低于天然林，有機物質投入量少［26］；當地的果園多種植桔科植物，因而地面凋落物較少，從而減少了土壤有機碳來源。而且，受耕作措施的影響，果園的深埋施肥對表層土壤的有機碳貢獻也不明顯，果園土壤連年深翻加速了土壤有機碳礦化過程［27-28］。另外，與農田地形較平坦相比，果園和林地常分布山地，具有較高的水力傳導性，土壤有機質容易隨著降雨而發生流失；這些原因也可能會導致SOCD的降低［29］。

4 結論

本研究中，由于BME-GWR和BME-GWRL模型結合了“軟數據”作為輔助變量，對研究區SOCD的空間異質性有更強的解釋能力，交叉驗證結果也反映BME-GWR和BME-GWRL模型較GWR模型的擬合精度更高，其中BME-GWRL模型基于土地利用類型的空間范圍對軟數據進行擬合，考慮到了“軟數據”估算單元的不確定性，對比不劃分土地利用類型直接模擬的BME-GWR預測結果，從整體上更加接近于觀測數據。本研究結果表明，考慮軟數據的獲取途徑，是在實測數據和輔助變量有限的條件下提高空間預測精度的有效途徑，為合理利用多源輔助數據提供新思路。