基于時空變異的旱地土壤有機碳高效采樣策略研究*

姚彩燕，劉紹貴，喬 婷，龍 軍，于東升，史學正，邢世和，陳瀚閱，張黎明?

（1. 福建農林大學資源與環(huán)境學院，福州 350002；2. 福建農林大學公共管理學院，福州 350002；3. 土壤生態(tài)系統(tǒng)健康與調控福建省高校重點實驗室，福州 350002；4. 揚州市農業(yè)環(huán)境監(jiān)測站，江蘇 揚州 225603；5. 土壤與農業(yè)可持續(xù)發(fā)展國家重點實驗室（中國科學院南京土壤研究所），南京 210008）

土壤有機碳（Soil Organic Carbon，SOC）作為全球陸地生態(tài)系統(tǒng)最大的有機碳庫，不僅是評價土壤肥力的重要指標，而且其含量的輕微變化也可導致大氣二氧化碳（CO2）濃度顯著變化，進而影響全球碳平衡[1]。受氣候、地形、土地利用和土壤類型等多方面因素綜合影響，SOC具有很高的時空變異性，不同時間和空間位置的SOC含量特征復雜且不確定，準確了解SOC時空分布變化，對于制定科學合理的碳匯管理措施具有重要意義[2]。樣點是土壤特性時空變異的信息載體，是土壤調查與空間制圖的基礎。在進行土壤采樣時，任何采樣過程都不可能獲得土壤的完全信息，只能對一定數量的土壤樣本數據通過插值形式實現(xiàn)土壤屬性空間分布的定量表達[3]。然而，樣點數量過小，土壤樣品的偶然性增加，缺乏代表性，空間預測將會產生不可接受的誤差；但盲目增大采樣數量，不僅預測誤差不會持續(xù)減小，反而加大了采樣工作量和樣品分析周期，浪費人力財力物力[4]。因此，明確不同樣點數量與SOC預測精度的響應關系是制定高效采樣策略的基礎。

近年來，國內外學者對不同區(qū)域的合理采樣數目開展了一系列研究[5-7]。但當前最佳土壤采樣數目的確定主要集中在若干個樣點數量的土壤屬性空間預測精度對比上，并沒有量化預測精度與采樣數量之間的關系。此外，目前的大多數研究一般局限在單一時間尺度下，未考慮不同時期土壤屬性差異對適宜采樣數目的影響。眾所周知，土壤是動態(tài)的，經過一定時間，SOC含量通常會出現(xiàn)一定范圍的變化，在某些情況下，SOC時間變異將會超過空間變異的影響[8]。因此，為了準確估計SOC含量，基于樣點的預測制圖必須要考慮空間和時間維度。旱地作為我國主要的農田類型，占全國耕地面積的76.9%，其SOC密度顯著低于世界平均水平，固碳潛力巨大[9]。但由于人類活動以及農田管理措施的影響，該土地利用類型下的土壤碳庫呈現(xiàn)強烈的時空變異性，明確不同時期旱地SOC預測精度與采樣數量之間的關系是合理制定我國農田管理措施的基礎。江蘇省北部地區(qū)85%的耕地為旱地土壤，盡管面積僅占全國旱地土壤總面積的3.7%，卻固定了全國10.2%的旱地SOC，具有較大的“碳匯”空間[10]。此外，該地區(qū)旱地的主要土壤類型為潮土、鹽土和砂姜黑土，分別占全區(qū)旱地土壤總面積的52.67%、24.31%和8.13%，這與黃淮海平原主要旱地土壤類型面積比例相一致，有利于研究成果向大區(qū)域拓展[11]?；诖耍狙芯恳越K省北部3.93×104km2旱地（簡稱“蘇北旱地”）為研究區(qū)，利用1980年全國第二次土壤普查和2008年農業(yè)部測土配方施肥項目的實測樣點數據，通過隨機抽樣的方式設置20組不同數量的樣點數據集，分析這兩個時期研究區(qū)SOC的最佳采樣數量，并在此基礎上揭示土壤屬性時空變異下預測精度與樣點數量之間量化關系，研究結果可為蘇北旱地和未來我國第三次土壤普查制定合理的有機碳野外調查方案提供理論依據。

1 材料與方法

1.1 研究區(qū)概況

蘇北旱地（116°21′～120°54′E，32°43′～35°7′ N）地處黃淮平原與江淮平原的過渡地帶，是中國沿海經濟帶的重要組成部分（圖1）。該地區(qū)屬于亞熱帶向暖溫帶過渡的季風性氣候，年平均降水量800～1 200 mm，年平均氣溫13～16 ℃，全年平均日照時數2 000～2 600 h，無霜期約220 d。境內地勢平坦，呈西北高、東南低的趨勢，地形主要以平原為主，海拔高度在0～50 m之間。研究區(qū)面積最大的旱地土壤類型為潮土和鹽土，其次為褐土、砂姜黑土和棕壤；而石質土、紫色土和石灰土面積較小。成土母質主要包括湖相沖積物、黃泛沖積物、河海相沖積物、其他各河流沖積物以及下蜀黃土等。蘇北是我國重要的糧食產區(qū)，種植制度主要是夏玉米和冬小麥輪作[12]。

1.2 數據來源與處理

1980年和2008年表層樣點數據（0～20 cm）分別來自全國第二次土壤普查中蘇北地區(qū)29個縣（市、區(qū)）在《縣級土種志》、《地市級土種志》、《省級土種志》和《中國土種志》中記錄的旱地剖面資料及農業(yè)部測土配方施肥項目。所有土壤樣點信息均包括地理位置描述和理化性質等，其中SOC含量可由有機質含量乘以Bemmelen轉換系數（0.58）計算?？紤]到樣點采集與化學分析過程的不確定性，對兩期樣點的SOC含量進行正態(tài)分布檢驗后進行異常值剔除處理，所采用的剔除方式是結合數據的常規(guī)統(tǒng)計學特征和空間位置，將每個樣點的SOC含量與總體及其鄰近8個樣點的均值和標準差進行比較，如果樣點SOC含量在總體均值的五倍標準差之外，且大于或小于鄰近樣點均值的三倍標準差，則視為異常值剔除，否則予以保留[13]。最后，1980年和2008年可利用的數據集分別為884個和1 081個（圖2）。

鑒于矢量數據庫可以準確表征SOC預測結果的空間分布特征，本研究使用1︰50 000土壤矢量數據庫實現(xiàn)柵格預測結果的輸出。1︰50 000土壤矢量數據庫是目前中國大區(qū)域尺度最為詳盡的土壤數據庫，由研究區(qū)第二次土壤普查各個縣（市）圖件資料GIS數字化而成。該數據庫旱地圖斑數目為17 024個，基本制圖單元為土種，圖例單元命名采用中國土壤發(fā)生分類系統(tǒng)。

1.3 樣點數量設計

不同樣點子集的創(chuàng)建采用ArcGIS 10.2地統(tǒng)計分析Create Subsets模塊。首先從2008年可利用樣點數據集中隨機剔除197個，剩余884個，保證1980年和2008年樣點總數一致。然后分別從兩期的884個樣點中隨機抽取85%的樣點（751個）作為訓練數據集參與實際的空間插值過程，剩余的15%（133個）作為驗證數據集用于插值預測精度的獨立驗證和不確定性評價。對751個訓練數據集按5%的遞減方式進行隨機抽取，形成其余19組訓練樣點子集通過SPSS軟件對各樣點數據集的SOC進行最大值、最小值、標準差（SD）、均值、變異系數（CV）、偏度、峰度等描述性統(tǒng)計和正態(tài)分布檢驗。

1.4 空間預測方法及其精度評價

本研究采用當前土壤學使用最多和效果最穩(wěn)健的普通克里金（Ordinary Kriging，OK）方法實現(xiàn)蘇北旱地SOC含量的空間預測?？臻g預測結果輸出為30 m×30 m的柵格，然后將每個土壤圖斑所包含柵格像元（即SOC含量預測值）平均值提取至矢量土壤圖圖斑。

本研究選擇常用的均方根誤差（Root Mean Square Error，RMSE）和相關系數（r）評價不同樣點數量下SOC空間預測的準確性。由于1980年和2008年SOC數據方差存在一定的波動性，RMSE和r不能直接用于不同時期預測精度比較。為了消除該影響，本研究也采用標準均方根誤差（Normalized Root Mean Square Error，NRMSE）比較不同時期插值預測精度，其中：

式中，n為驗證樣點數量（本研究n＝133），Xmi和Xpi分別表示第i個驗證樣點中SOC含量的實測值和預測值。RMSE或NRMSE越小，r越大，預測精度越高。

1.5 經典統(tǒng)計學確定合理采樣數目的計算公式

本研究所獲取的最佳采樣數量與經典統(tǒng)計學確定的合理采樣數量進行對比，以分析其合理性。在目前大多數研究中，Cochran[14]針對與區(qū)域范圍內純隨機采樣而構造的計算公式通常被經典統(tǒng)計學用于最佳采樣數量的確定：

式中，n為所需要的最佳采樣數量，t為與顯著水平對應的t式分布值，Std為樣本標準差，d為樣本平均值與允許誤差的乘積。

2 結果與討論

2.1 不同樣點數量下1980年和2008年SOC含量統(tǒng)計學特征及空間變異結構

基于訓練樣點總集（751個樣點）的SOC含量統(tǒng)計結果表明（表1），1980年研究區(qū)SOC含量均值為5.70 g·kg-1，而2008年含量均值明顯提高，達到10.41 g·kg-1。兩個時期其余19組樣點子集SOC含量均值、變異系數與訓練樣點總集統(tǒng)計結果相似，表明隨機抽取的不同樣點數量均具有較好代表性，降低采樣數量對研究區(qū)SOC基本統(tǒng)計信息影響較小，這與程道全等[15]和李凱等[16]的研究結果一致。進一步的Kolmogorov-Smirnov檢驗表明，1980年所有樣點子集SOC含量全部符合正態(tài)分布，2008年樣點子集SOC含量經對數轉換后，也均符合正態(tài)分布，滿足地統(tǒng)計插值的要求。

表1 不同樣點數量下研究區(qū)SOC含量（g·kg-1）描述性統(tǒng)計結果及半方差模型擬合參數Table 1 Descriptive statistical and semivariance models fitting parameters of SOC contents (g·kg-1) in the study area relative to number of sampling sites

表征SOC空間變異的指示參數（塊金值、基臺值、結構性比例和變程）在不同時期或樣點數量下具有較大差異，表明SOC空間變異結構隨著樣點數量以及時間的推移發(fā)生明顯變化。1980年研究區(qū)SOC的結構性比例較小（3.43%～14.19%），其變異性受結構性因素影響在空間上強烈自相關；而2008年SOC結構性比例在14.29%～50.00%之間，絕大部分樣點子集都超過了25%，土壤變異的空間相關性減弱為中等水平，這是因為1980—2008年研究區(qū)的耕作、施肥等土壤管理措施經歷了較大改變，對SOC空間變異結構造成了一定影響。

2.2 不同樣點數量下1980年和2008年SOC預測結果的空間分布特征

隨著樣點數量的減少，蘇北旱地SOC的空間分布特征發(fā)生明顯改變（圖3和圖4）。訓練樣點總集下，SOC空間分布既能反映整體格局又能體現(xiàn)變異細節(jié)，但隨著土壤樣點數量的降低，表征SOC局部變異的能力減弱，并逐漸改變整體空間分布格局，其中SOC含量高值和低值區(qū)域面積變化最為劇烈。造成不同數據集極值區(qū)域估算結果明顯差異的原因可能是因為兩端的極值具有較高空間異質性，樣點數量減少導致表征土壤空間異質性的能力減弱。一方面，隨著樣點數量的降低，樣點數據集的SOC含量最大值呈現(xiàn)減小趨勢，1980年和2008年的變幅分別由751個樣點的13.17和20.65 g·kg-1減小至38個樣點的11.37和15.78 g·kg-1；而最小值呈增大趨勢，1980年和2008年的變幅分別由751個樣點的0.70和5.86 g·kg-1增大至38個樣點的1.68和6.09 g·kg-1，表明樣點數量減少會使得所探測到的SOC含量范圍降低，進而影響其表征研究區(qū)的空間完整性。另一方面，巫振富等[17]研究指出，當最能凸顯土壤空間分布格局的高、低值“關鍵區(qū)域”樣點數量足夠多時，才可保證預測結果最大程度地反映研究區(qū)土壤屬性的實際空間分布特征，但樣點數量的減少會導致“關鍵區(qū)域”能否滿足樣點需求成為一種不確定性事件。

續(xù)表

以訓練樣點總集空間預測結果的高、低值區(qū)為基準，劃定研究區(qū)SOC空間分布的“關鍵區(qū)域”。經統(tǒng)計，751個訓練樣點總集下1980年和2008年SOC高值區(qū)的樣點數量分別為147和40個，低值區(qū)分別為130和372個。當樣點數量依次減少為638（85%）、563（75%）、451（65%）個時，位于“關鍵區(qū)域”的土壤樣點數量仍然較多，其中1980年SOC高值區(qū)樣點數量分別為121、116和97個，低值區(qū)樣點數量分別為111、100和81個；2008年SOC高值區(qū)樣點數量分別為33、28和29個，低值區(qū)樣點數量分別為320、287和242個，盡管局部的一些細節(jié)信息逐漸被隱藏，但此時的SOC的整體分布格局和趨勢仍然比較清晰，樣點數量降低并未導致預測制圖結果差異明顯。然而，當樣點數量進一步減少至113（15%）、38（5%）個時，位于SOC高值區(qū)的數量1980年分別僅有28和3個，2008年分別僅有8和2個；位于低值區(qū)的數量1980年分別僅有12和7個，2008年分別僅有48和17個，樣點數量無法滿足“關鍵區(qū)域”的制圖需求，導致SOC在空間分布上的細節(jié)信息基本被平滑，呈現(xiàn)出單一的分布格局，不能很好地反映蘇北旱地SOC的真實分布狀況，特別是5%樣點數量下空間結構出現(xiàn)失真畸變，表現(xiàn)的空間分布格局與訓練樣點總集相異。

從圖2也可以看出，雖然2008年的采樣點空間分布較不均勻，但該時期各個樣點集的空間分布差異整體較小，這是由于其分布于“關鍵區(qū)域”的樣點數較多，能更穩(wěn)定地反映“關鍵區(qū)域”空間分布的真實情況，隨著樣點數量的降低，樣點集統(tǒng)計特征的變異系數變化程度也明顯小于1980年。這表明空間插值精度受采樣數量和空間分布的共同影響，位于“關鍵區(qū)域”的樣點數量越多，其估算結果越穩(wěn)定。楊琳等[3]也認為，與單純地增加樣點數量相比，基于樣點代表性的土壤采樣在土壤空間制圖成本和準確性上具有更為明顯的優(yōu)勢。因此，在SOC高、低值等“關鍵區(qū)域”設置足夠的樣點數量，從而提高采樣點所反映的空間異質性，是提高土壤空間預測效果的重要手段之一。

2.3 不同樣點數量下1980年和2008年SOC空間預測精度

1980年和2008年研究區(qū)SOC含量空間預測精度隨樣點數量變化如圖5所示。在1980年預測制圖結果中，r和RMSE隨樣點數量變化的特征可劃分為：100%～75%、70%～20%和15%～5%三個階段。在第一階段，r和RMSE變化幅度較小，預測精度并未隨著采樣點數量的減少而產生明顯差異，說明當樣點數量≥75%時，通過OK插值手段獲得的SOC含量預測精度較高且已經達到相對穩(wěn)定水平，研究區(qū)的最佳采樣數目在75%（563個）左右，這可能是因為樣點數量增加至一定程度以后，空間上的聚集度較高，導致信息達到一定的飽和水平。在這種情況下，即使采樣點數量增加，也無法提供進一步的信息。在第二階段，隨著樣點數量的減少，SOC含量的空間預測精度逐漸降低，r和RMSE變化幅度分別介于0.39～0.56和2.09～2.35 g·kg-1之間，原因在于預測精度受到樣點代表性和密度顯著影響。根據極值區(qū)域的樣點統(tǒng)計，第二階段的數量由216個減少至63個，SOC空間數據的代表性顯著降低。此外，ArcGIS 10.2的空間統(tǒng)計分析表明，第二階段樣點之間的平均最近鄰距離為4～9 km，明顯大于前一階段。在第三階段，r和RMSE隨著樣點數量的減少急劇變化，預測精度較第二階段降低得更快。半方差模型擬合結果顯示，15%～5%樣點數量下的擬合優(yōu)度僅在0.08左右，SOC空間變異結構的推斷存在較大誤差，這可能是預測精度急劇降低的主要原因。然而，與1980年相比，使用20組相同樣點梯度的訓練子集數據對2008年研究區(qū)SOC含量預測時，預測精度隨樣點數量變化特征存在一定差異。2008年預測精度的變化趨勢同樣表現(xiàn)為先緩慢再劇烈下降，但階段劃分有所不同，主要分為：100%～70%、65%～20% 和15%～5%三個。在每個階段，預測精度特征與1980年也不盡相同。第一階段，預測精度并未隨著樣點數量的降低表現(xiàn)出顯著變化，即70%的樣點數時預測效果能達到相對穩(wěn)定水平。可見，2008年研究區(qū)的最佳采樣數目在526個左右，較1980年少5%，說明隨著樣點數量的增加，其預測精度更快趨于穩(wěn)定。第二階段，預測精度隨樣點數量降低呈現(xiàn)出明顯的下降趨勢，r和RMSE變化幅度分別處于0.42～0.63和2.11～2.56 g·kg-1之間。在該階段內，2008年各個樣點梯度間的預測精度差距明顯大于1980年，表明樣點數量的增加對其預測效果的提高更為顯著。第三階段，隨著樣點數量減少，樣點梯度間預測精度差距進一步增大，即對樣點數量的敏感性增強。但是，與1980年不同的是，2008年預測精度在15%樣點數量下出現(xiàn)了一定回升，原因可能與樣點數量稀少而引起的預測精度不確定性有關。研究表明，預測精度受到樣點數量和空間分布的共同影響，隨著樣點數量的減少，空間分布位置對于預測效果的影響逐漸增大，甚至可以取代樣點數量成為預測效果的主導性因素，但隨機抽樣忽略了樣點的空間布局，使得樣點是否對研究區(qū)SOC空間分布具有較好的全局代表性成為一種概率事件[17]。

從圖5也可以看出，1980年和2008年蘇北旱地SOC含量空間預測精度隨樣點數量變化表現(xiàn)出不同的響應特征，每個時期所需的樣點數量有所差異。這主要是由于不同時期的環(huán)境、氣候以及農田管理措施（施肥、灌溉和耕作等）不同，導致有機碳在土壤中的含量及其空間結構不同而產生響應特征差異。由于1980年研究區(qū)SOC空間自相關變程較小、局部變異較多，不僅需要更多的樣點來推測空間分布，并且其預測精度的提高速率也較為緩慢和不穩(wěn)定；相比而言，2008年研究區(qū)SOC空間自相關變程較大、局部變異較少，所以預測精度對于采樣數量的敏感性較強，空間信息達到飽和狀態(tài)所需樣點數量也相對較少。從NRMSE來看，不同時期SOC含量的預測效果差距也十分明顯，20組樣點數量下1980年和2008年的NRMSE分別為0.34～0.43和0.20～0.25 g·kg-1，在樣點數量相同的情況下2008年的空間預測精度優(yōu)于1980年，主要的原因有以下3點：一是空間變異性對OK插值精度具有負面影響，該方法多適用于土壤屬性變異較弱的環(huán)境區(qū)域[18]；二是由于不同時期樣點的空間布局差異所致；三是2008年的土壤空間樣點盡管分布不均，但位于關鍵區(qū)域的樣點數量較多，其半方差模型擬合參數明顯優(yōu)于1980年，采樣點具有更好的代表性。由此可見，研究區(qū)最優(yōu)采樣點數量及其預測精度不是固定的，而是由每個階段土壤屬性的空間變異特征以及樣點的分布位置、空間布局共同決定。

根據經典統(tǒng)計學確定合理采樣數目的計算公式（表2），在95%的置信區(qū)間和5%的誤差條件下，1980年和2008年分別只需要261和100個采樣點即可滿足空間預測需求，而在90%的置信區(qū)間下，分別僅需20和8個采樣點即可滿足15%誤差條件下的空間預測需求。本研究顯示，1980年和2008年的最優(yōu)采樣數量應分別為563和526個，明顯多于利用Cochran計算公式得到的最佳采樣數量。Cochran在描述土壤元素空間變異時假定研究區(qū)土壤元素的含量變化是隨機的，樣本之間服從同一已知的概率分布并且相互獨立，認為樣品均值是反映研究區(qū)土壤特性的最優(yōu)方式，決定土壤最佳采樣數量的影響因素主要包括土壤樣本的變異系數、允許誤差和可接受的置信區(qū)間。然而，地統(tǒng)計研究表明土壤特性的變化并非完全隨機，其在空間上具有一定的結構性和關聯(lián)性，Cochran計算公式雖然能夠顯著減少采樣數量，但只能概括SOC含量變化的宏觀趨勢而不能反映其局部變異特征，如果盲目采用Cochran公式確定采樣數量，可能會造成超過預期的預測誤差和不確定性。本研究采用隨機抽樣與地統(tǒng)計學的克里金插值相結合的方法確定合理采樣數目，充分考慮了SOC的空間變異和結構特征，較Cochran將研究區(qū)作為一個均質的整體來確定最優(yōu)采樣數目的方式更具科學性。

表2 基于Cochran采樣數量公式獲取的蘇北旱地SOC最佳采樣數Table 2 Optimal number of sampling sites for prediction of SOC in the upland area in North Jiangsu based on the Cochran sampling quantity formula

3 結 論

本研究通過隨機抽樣方式設置20組不同數量的樣點數據集，結合1︰50 000高精度土壤矢量數據庫明確樣點數量對1980年和2008年SOC普通克里金插值預測精度的影響，結果表明，1980年蘇北旱地SOC預測值與實測值的相關系數和均方根誤差分別在0.15～0.56和2.09～2.63 g·kg-1之間，樣點數量約563個時可達到相對可靠的預測水平；2008年SOC預測值與實測值的相關系數和均方根誤差分別在0.24～0.63和2.11～2.62 g·kg-1之間，空間信息達到穩(wěn)定預測水平所需樣點數量約526個。這表明相同區(qū)域范圍內的最優(yōu)采樣數量及其預測精度不是固定的，是由每個階段內土壤屬性的空間變異特征以及樣點的分布位置、空間布局共同決定。因此，未來的土壤采樣設計需要根據研究區(qū)域土壤屬性的時空變異特點制定合理的樣點數量以及布局方式，以提高土壤采樣的科學性和高效性。此外，有研究表明系統(tǒng)抽樣能夠較空間隨機抽樣更好地反映土壤屬性的空間變異，未來的研究還應進一步考慮系統(tǒng)抽樣方法對土壤屬性變量空間預測結果的影響[19-20]。