丁亞鵬,張俊華,3,*,劉玉寒,盧翠玲,王爍騫,秦靜婷,丁圣彥,3

1 黃河中下游數字地理技術教育部重點實驗室,河南大學,開封 475004 2 河南大學環(huán)境與規(guī)劃學院,開封 475004 3 河南省大氣污染綜合防治與生態(tài)安全重點實驗室,開封 475004

土壤作為陸地生態(tài)系統(tǒng)中最大的碳庫[1-2],在全球碳循環(huán)中占據重要作用[3-5]。土壤中的碳素主要以有機碳的形式存在[6],土壤有機碳含量及其儲量的變化都會對生態(tài)系統(tǒng)的穩(wěn)定及功能發(fā)揮產生重要影響。土壤有機碳含量受土地利用方式[7-8]、管理措施[9-10]、群落結構[11-12]、地形地勢[13-14]等諸多因素的綜合影響,在空間上具有明顯的差異性[13,15]。土壤有機碳含量的微小變化會對生態(tài)系統(tǒng)和全球氣候變化產生重要影響[16-17]。對土壤有機碳空間分布特征及影響因素的研究有助于認清土壤有機碳變化過程及區(qū)域化反應。同時,能夠加強對土壤碳庫的科學認識,進一步了解土壤碳庫"源"和"匯"的關系,在土壤有機碳的提升和土地資源管理方面有重大意義。

傳統(tǒng)方法采用空間地統(tǒng)計分析對土壤有機碳空間分布進行研究,該方法只考慮樣點之間距離的空間關聯性,而忽視了各影響因素的貢獻率[18-19]。探究土壤有機碳影響因素的研究多采用最小二乘法(Ordinary least squares，OLS)和逐步回歸等方法[20-21]。這些線性回歸模型都是基于全局回歸模型,來反映區(qū)域總體狀況,其認為各要素對空間上所有區(qū)域的影響是相同的,忽略了各影響因子的局部性,導致模擬結果與實際情況相差較大[22-23]。地理加權回歸模型(Geographically weighted regression,GWR)是一種基于樣點地理位置的局域空間分析方法,能夠將局部范圍的因變量和解釋變量進行合并,從而得到每個樣點各解釋變量的回歸系數。與傳統(tǒng)回歸模型相比,該模型在研究土壤有機碳的空間建模中具有更高的精度。因此該方法在研究土壤有機碳分布及與環(huán)境因素、土壤性質的關系方面得到了廣泛應用[13,24-25]。受研究區(qū)范圍影響,在空間大尺度上側重環(huán)境因素和土壤有機碳的關系[13],在小尺度上側重土壤性質對有機碳空間特征的影響[26],綜合考慮環(huán)境因素和土壤理化因子對有機碳的共同作用的成果較少。

流域是社會-經濟-自然要素綜合作用的復合體,其組成要素的多樣性決定了流域的復雜性。在流域尺度上探討各因素對土壤有機碳的影響,是對開展區(qū)域生態(tài)學問題研究、維持各生態(tài)系統(tǒng)良性發(fā)展、綜合考慮各類要素間利用和協調發(fā)展問題進行的有益探索。本文以伊河流域土壤為例,利用GWR模型探索土壤有機碳與環(huán)境因子和其它土壤性質之間的關系,揭示伊河流域土壤有機碳的空間分布特征及影響因素,以期為伊河流域土地利用和管理提供一定的科學依據。

1 材料與方法

1.1 研究區(qū)概況

伊河流域位于河南省西部山區(qū)(33°39′—34°41′N、111°19′—112°54′E),地處我國二、三階梯的過渡地帶。伊河發(fā)源于熊耳山,流域總面積約6100 km2[27],海拔88—2128 m。地勢西南高,東北低,地貌類型豐富,主要包括山地、丘陵和平原。該區(qū)域處于北亞熱帶向暖溫帶過渡區(qū),屬季風性氣候,年均溫12.4—15.1℃,年均降水量700—900 mm[28],四季變化明顯。伊河流域雨熱同期,降水多集中于夏季,且降雨強度較大,土壤易受侵蝕；流域內上游降水較為充足,其年降水量約是下游的2倍。伊河流域土壤類型較多,分布最廣的是褐土,同時還存在棕壤、紅黏土、沙土等。該流域植被類型以暖溫帶落葉闊葉林為主,在海拔較高地區(qū)存在針葉、落葉闊葉混交林,平原區(qū)和低山丘陵區(qū)是人類活動的主要場所,大面積的低山丘陵被開發(fā)為農田,自然植被分布較少。

伊河流域總人口237.2萬,城鎮(zhèn)化率42.68%,人口集中分布在河流沿岸,尤其是中下游的丘陵和平原區(qū)[28]。伊河流域土地利用類型以耕地和林地為主(圖2),城鎮(zhèn)和村落等建設用地面積約325.10 km2,主要分布在河流沿岸的丘陵和平原區(qū)；耕地面積約2458.87 km2,主要分布在丘陵和平原區(qū),在低山區(qū)也有少量分布；林地面積約2714.57 km2,主要分布在低山和中山區(qū),在丘陵和平原區(qū)有零星分布；草地面積約360.82 km2,主要分布在低山區(qū),在其它地區(qū)有零星分布。

圖2 研究區(qū)土地利用類型圖Fig.2 Land use type map of study area

1.2 樣品采集和指標測定

根據研究區(qū)面積和交通可達性,對研究區(qū)采用6 km×6 km網格取樣,采集表層(0—20 cm)土壤樣品141個,其中農田樣品76個,森林樣品59個,草地樣品6個。將采集的土樣帶回室內自然風干,取適量土樣研磨過篩(2 mm、0.25 mm、0.15 mm)備用。土壤理化性質測定方法參照《土壤農業(yè)化學分析方法》[29],具體方法如下(表1)：

表1 土壤理化性質測定方法

1.3 研究方法

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

本文以伊河流域DEM數據為基礎,通過ArcGIS 10.2空間分析模塊對DEM數據進行地形指標的提取計算。共計算了11個地形因子,包括：海拔、坡度、坡向、曲率、平面曲率、剖面曲率、地形起伏度、地表粗糙度、復合地形指數(CTI)、匯流動力指數(SPI)和沉積物運移指數(STI)。其中CTI、SPI和STI分別由公式(1)、(2)、和(3)計算獲得。研究區(qū)DEM來源于91衛(wèi)圖的GoogleEarth高程數據,空間分辨率為15 m。選取的2個氣象因子(年均溫、年均降雨量)通過逐日氣象數據計算獲得,數據來源于中國氣象局(http://data.cma.cn/)。另外選取與土壤樣品采集同時期的歸一化植被指數(NDVI),來反映研究區(qū)植被覆蓋情況,數據來源于美國NASA網站(https://www.nasa.gov/),空間分辨率為250 m。將上述14個因子作為探討影響土壤有機碳的環(huán)境因子。

(1)

SPI=ln(AC×tanβ×100)

(2)

(3)

式中,AC為垂直于水流方向的特定匯流面積,β為坡度。

1.3.2GWR模型

最小二乘法模型(OLS)是隨機變量(y)與確定性變量(x1,x2、xi)的多元線性函數,這是一種基于全局回歸的函數(式4)。GWR模型[30]是對OLS模型的拓展,是一種局部回歸模型,將數據的地理位置嵌入到回歸參數之中,可以實現對參數的局部估計(式5)。

(4)

(5)

式中,yi為樣點i的因變量,xik為第i個點上第k個變量的觀測值,(μi,νi)為第i個點的位置坐標,β0(μi,νi)為截距,βk(μi,νi)為第i個的回歸系數,εi為誤差項。

GWR模型中參數的估算關鍵在于空間權重函數的選取,本文比較4種常用空間權重函數,以確定適合本區(qū)域的最優(yōu)空間函數,分別是固定高斯函數(Fixed Gaussian)(式6)、自適應高斯函數(Adaptive Gaussian)(式7)、固定截尾型函數(Fixed bi-square)(式8)和自適應截尾型函數(Adaptive bi-square)(式9)。

(6)

(7)

(8)

(9)

式中,ij代表給定研究區(qū)域的任意一點,θ>0代表窗寬或者光滑參數。

帶寬的選擇是影響GWR模型分析結果的關鍵因素,帶寬的大小會直接影響回歸參數的估計。本文采用黃金分割搜索來選擇帶寬,以AICc作為帶寬選擇準則。

AIC=2K-2ln(L)

(10)

(11)

其中,K是參數數量,L是似然函數。

1.4 數據處理

實驗數據在Excel中記錄整理；地形因子通過ArcGIS 10.2提取、計算；GWR模型在GWR 4.0中運行；文中空間分布圖均采用WGS_1984地理坐標系；相關分析、多元線性回歸和共線性診斷在R軟件中實現,相關分析采用斯皮爾曼(Spearman)法(P<0.05)。

2 結果與分析

2.1 伊河流域土壤有機碳的描述性統(tǒng)計特征及與環(huán)境變量和土壤性質的相關性

伊河流域表層土壤有機碳范圍在3.37—38.34 g/kg,平均含量為12.23 g/kg,變異系數為0.47,屬于中等變異水平(表2)。流域上、中、下游有機碳平均含量分別為12.43 g/kg、12.48 g/kg和11.72 g/kg,變異系數分別為0.59、0.46、0.29,均屬于中等變異水平,變異系數自上游到下游逐漸降低,表明伊河流域土壤有機碳存在空間分布差異,其中上游差異最大,下游差異最小。

表2 伊河流域表層土壤有機碳含量特征/(g/kg)

將土壤有機碳與環(huán)境因子和其它土壤性質進行相關分析(圖3)。環(huán)境因子中的年均氣溫與土壤有機碳呈顯著負相關(P<0.05),其它環(huán)境因子與有機碳的相關性較弱；土壤因子中的容重、銨態(tài)氮與土壤有機碳表現出顯著負相關(P<0.01)的趨勢；50—250 μm和250—1000 μm粒級含量、活性有機碳、全氮、硝態(tài)氮和全磷與有機碳呈正相關,且相關性顯著(P<0.01)；pH和速效磷與有機碳相關性不顯著。

圖3 土壤有機碳與環(huán)境因子和土壤性質的相關分析Fig.3 Correlation analysis of soil organic carbon with environmental factors and soil propertiesTOC: 有機碳; CTI: 復合地形指數; SPI: 匯流動車指數; STI: 沉積物運移指數; NDVI: 歸一化植被指數; AOC: 活性有機碳; TN: 全氮; TP: 全磷; AP: 速效磷

2.2 GWR模型構建

2.2.1模型變量選擇

相關性作為模型解釋變量選擇的初步條件,通過相關分析,初步篩選出年均氣溫、活性有機碳、全氮、硝態(tài)氮、銨態(tài)氮、容重、全磷、50—250 μm和250—1000 μm粒級含量9個指標作為土壤有機碳的解釋變量。為減少解釋變量存在共線性引起GWR預測產生偏差,需對解釋變量進行共線性診斷,并篩選剔除存在共線性的變量。本文選取方差因子(Variance inflation factor,VIF)、特征值、條件索引和方差比例來判定各解釋變量之間是否存在共線性。由表3發(fā)現,各解釋變量VIF值均小于10,在所有主成分中,只有9和10的條件索引大于30,且每一主成分各解釋變量的方差比例只存在單一變量的方差比例大于0.5,加之最后一個主成分不能提取出來信息,特征值接近0,條件索引相應會很大。綜上所述,年均氣溫、活性有機碳、全氮、硝態(tài)氮、銨態(tài)氮、容重、全磷、50—250 μm和250—1000 μm粒級含量均為土壤有機碳的解釋變量。

表3 土壤有機碳解釋變量間的共線性診斷

2.2.2模型診斷與選擇

以土壤有機碳為因變量,年均氣溫、活性有機碳、全氮、硝態(tài)氮、銨態(tài)氮、容重、全磷、50—250 μm和250—1000 μm粒級含量為自變量,選擇不同的空間權函數進行建模。在模型運行時需對模型的擬合效果進行診斷和選擇,選擇最優(yōu)模型。診斷參數包括殘差平方和、 AICc、R2和調整后R2[31]。不同權函數GWR模型殘差平方和和AICc值均小于OLS,R2和調整后R2均高于OLS(表4),表明對于本研究區(qū)來說,GWR模型優(yōu)于OLS模型。高斯空間權函數的殘差平方和、AICc值均高于截尾型函數,R2和調整后R2均低于截尾型函數。自適應截尾型空間權函數的殘差平方和和AICc值最低,分別為1945.76和825.20,R2與調整后R2最大,分別為0.58和0.49,自適應截尾型空間權函數是本研究區(qū)的最優(yōu)空間權函數,其模型精準度最高。

表4 不同權函數GWR模型與OLS模型對比

2.3 GWR模型土壤有機碳解釋因子系數空間分布

GWR模型預測伊河流域土壤有機碳空間分布整體結果較好,局部決定系數在0.49—0.64之間,自下游到上游,決定系數逐步提高,對上游的預測精度更高(圖4)。土壤有機碳各解釋因子回歸系數具有一定差異?；貧w系數最大的是AOC,范圍在1.28—2.98,表明土壤有機碳受AOC影響最大,容重的回歸系數絕對值最小,在-0.45—-0.01之間,表明有機碳受容重影響較小。各解釋因子空間分布具有局部性和非均勻性,表明同一因子對不同地點有機碳的影響程度不同。年平均氣溫對土壤有機碳的影響整體均為負效應,在伊河流域年平均氣溫越高,土壤有機碳含量越低,其系數絕對值中、上游最大,下游最小,年平均氣溫對中、上游有機碳限制性影響最大,下游影響微弱。全氮對有機碳的影響為正效應,全氮含量越高,有機碳含量越高,全氮系數自上游到下游逐步增加,下游有機碳受全氮影響最大,上游全氮系數在0.173—0.607,上游有機碳受全氮影響較弱。活性有機碳對有機碳是正效應,系數自上游到下游逐步減小,活性有機碳對上游有機碳的影響最大,中游次之,下游最小。硝態(tài)氮的系數在0.18—1.61之間,自下游到上游逐步增加,上游系數較小,硝態(tài)氮對上游有機碳影響較弱,對下游影響較大。銨態(tài)氮系數在-1.13—0.58之間,銨態(tài)氮對有機碳是一個負效應,在中游地區(qū)最為明顯,上游和下游系數較小,影響較弱。容重整體系數較小,對有機碳影響較弱,也存在空間差異,自上游到下游,對有機碳的負效應逐步加大。全磷對有機碳的影響在空間上有較大差異,系數范圍在-0.18—1.29之間,在下游地區(qū),全磷系數為負值,全磷對有機碳影響為負效應,中游和上游其系數為正值,對有機碳影響為正效應,在上游地區(qū)全磷系數較大,上游全磷對有機碳的正效應影響最大。土壤顆粒組成對土壤有機碳的影響在空間上分異較大。其中50—250 μm粒級系數在中上游為負值,表明對有機碳的影響為負效應,但絕對值較小,說明在對中上游有機碳影響較弱,其系數下游為正值,表明對下游有機碳的影響為正效應。250—1000 μm粒級系數范圍在0.83—1.68,下游系數較小,中游最大,表明250—1000 μm粒級對下游有機碳影響較小,對中游和上游影響較大。

圖4 土壤有機碳解釋變量回歸系數分布圖Fig.4 Distribution chart of regression coefficient of soil organic carbon explanatory variables

2.4 土壤有機碳空間分布特征

GWR模型預測伊河流域土壤有機碳含量在4.31—26.93 g/kg之間(圖5),土壤有機碳含量的最大值和最小值都在上游地區(qū),且低值區(qū)占的面積較大,說明上游地區(qū)有機碳含量差別較大,空間變異性最大。下游地區(qū)有機碳含量較接近,大多在10.95—22.03 g/kg之間,空間變異性較小。中游地區(qū)高值區(qū)、中值區(qū)和低值區(qū)范圍接近,空間變異性居中。總體來說,土壤有機碳空間分布的高值區(qū)大多集中在海拔較高的中山區(qū)和地勢平坦的平原區(qū),低值區(qū)主要分布在低山、丘陵和河流沿岸。

圖5 土壤有機碳空間分布圖Fig.5 Spatial distribution map of soil organic carbon

3 討論

環(huán)境因子影響土壤有機碳的空間分布[13-14,24,32],研究發(fā)現,地形因子、NDVI和氣象因子中的年均降水量與有機碳相關性不顯著,與羅梅等[13]研究年均降水量、歸一化植被指數、高程以及地形粗糙指數與有機碳含量呈極顯著相關關系有差別,可能是因為研究區(qū)人類活動劇烈,從而干擾了環(huán)境因子的作用,加之研究區(qū)范圍較小,環(huán)境因子作用不明顯。從土壤有機碳的空間分布來看,與地形有較為一致的規(guī)律,有機碳分布的高值區(qū)大多在海拔較高地區(qū)。高海拔地區(qū)氣溫較低,植被覆蓋以灌木和草本為主,大量枯枝落葉堆積,對降水的緩沖截留,避免了土壤有機碳的流失。但從相關分析來看,有機碳的空間分布與海拔關系較弱,這方面的原因可能是由于采樣造成的,海拔較高的地方難以到達,采樣點較少,從而影響了海拔與有機碳的相關關系。研究發(fā)現,年平均溫度對土壤有機碳是負相關,這與吳春生等[33]研究一致,氣溫降低,土壤有機碳釋放速率降低,碳氮礦化速率減慢,同時環(huán)境微生物活性降低,分解動植物殘體速度減慢,致使土壤有機碳和全氮含量積累量顯著[34-35]。

自然因素和人類活動綜合影響土壤有機碳的空間分布特征,且具有明顯的地域性。在海拔較高的中上游區(qū)域,土壤有機碳含量主要受年平均氣溫、活性有機碳和全磷含量的影響?；钚杂袡C碳主要來源于土壤有機碳的水解,植物根系的分泌物、及其凋落物的分解,和土壤微生物本身及其代謝產物[36-37]；而全磷主要來源于成土母質[38]。這表明在自然狀態(tài)下,土壤有機碳主要受立地環(huán)境、成土母質和地表覆蓋的影響。在中上游的低山丘陵區(qū),土壤有機碳主要受容重、250 μm以上粒級含量和速效氮的影響。低山丘陵區(qū)大多被改造為耕地,受季風性氣候的影響,降水集中,多暴雨,土壤侵蝕嚴重,細顆粒物、速效氮易被淋溶、沖刷,是引起土壤中碳含量減少的原因,加之低山丘陵區(qū)耕作交通不便,化肥和有機肥難以大量施用,收獲的秸稈主要被用來當作燃料,使得有機碳歸還較少,也就導致該區(qū)有機碳含量較低。這表明在低山丘陵區(qū),人類活動和環(huán)境因素,尤其是在地形和氣候的共同作用影響了土壤有機碳的含量。平原區(qū)有機碳含量較高,主要受全氮和硝態(tài)氮的影響,可能是因為平原區(qū)主要是耕地,地勢平坦,交通便利,大量化肥和有機肥的投入,從而使得土壤有機碳含量較高[9]。

GWR模型預測土壤有機碳可以恰當地反映出因子在局部區(qū)域的影響,并且其預測效果在自然狀態(tài)下或人為干擾較弱的地區(qū)更好。通過對比四種空間權重矩陣函數結果,GWR模型在空間權重矩陣函數的選擇上,自適應型函數可以根據采樣點之間的距離調整權重,更具有合理性[39]。該模型在土壤方面應用的最大的問題之一就是解釋變量的局部共線性[40],土壤性質間相互作用,有較強的相關性[41],很容易產生局部共線性問題,這是現階段GWR模型在土壤方面應該被重視的問題之一,也是今后該模型需要解決的主要問題。

4 結論

伊河流域表層土壤有機碳范圍在3.37—38.34 g/kg,流域上、中、下游有機碳平均含量分別為12.43 g/kg、12.48 g/kg和11.72 g/kg,伊河流域土壤有機碳存在空間分布差異,其中上游差異最大,下游差異最小。相關分析表明,土壤有機碳含量受年平均氣溫以外的環(huán)境因子影響較小,與土壤理化性質相關性顯著。

GWR模型預測土壤有機碳可以恰當地反映出各因素在局部區(qū)域的影響。結果表明,在海拔較高的中上游區(qū)域,土壤有機碳含量較高,主要受立地環(huán)境、成土母質和地表覆蓋的影響。在中上游的低山丘陵區(qū),人類活動和環(huán)境因素共同影響下土壤有機碳含量較低。中下游平原區(qū)農業(yè)活動和化肥投入是影響土壤有機碳含量較高的原因。這為伊河流域各亞區(qū)生態(tài)系統(tǒng)的合理發(fā)展和管理提供了依據。

GWR模型預測伊河流域土壤有機碳空間分布整體結果較好,局部決定系數在0.49—0.64之間,自下游到上游,決定系數逐步提高,對上游的預測精度最高。該模型能夠較好預測伊河流域土壤有機碳的空間分布特征,對人為干擾較弱的區(qū)域預測效果更好。