徐占軍，張 媛，張紹良，李樂樂，余明成

（1. 山西農(nóng)業(yè)大學(xué)資源環(huán)境學(xué)院；晉中 030801；2. 中國(guó)礦業(yè)大學(xué)環(huán)境與測(cè)繪學(xué)院；徐州 221008）

0 引 言1

2014年，中國(guó)煤炭產(chǎn)量為2.6×109t，占世界煤炭總產(chǎn)量的 46.9%[1]。煤炭開采過程中的一系列環(huán)節(jié)占用甚至損毀了大量的土地資源，導(dǎo)致地上植被被移除，土壤動(dòng)物和其他動(dòng)物因生境喪失而消亡，礦區(qū)的農(nóng)作物減產(chǎn)，生態(tài)環(huán)境遭到嚴(yán)重破壞[2]。而且，煤礦開采塌陷的土地大部分是平原農(nóng)業(yè)區(qū)的可耕地，即華北平原、東北平原和長(zhǎng)江中下游平原的農(nóng)業(yè)區(qū)[3]。相關(guān)研究表明，中國(guó)年塌陷的耕地面積約為200 km2[4]。人類煤炭開采活動(dòng)對(duì)農(nóng)業(yè)生態(tài)系統(tǒng)危害嚴(yán)重[5]，對(duì)農(nóng)田土壤有機(jī)碳碳庫(kù)擾動(dòng)十分劇烈[6]。

由于農(nóng)田的土壤有機(jī)碳庫(kù)是減少陸地生態(tài)系統(tǒng)碳排放的最大潛在因素[7-8]，中國(guó)以及世界上的其他煤炭開采大國(guó)必須在區(qū)域尺度上就煤炭開采對(duì)區(qū)域農(nóng)田碳儲(chǔ)量的影響進(jìn)行定量的研究，才可以更好地對(duì)煤炭開采區(qū)的土壤有機(jī)碳庫(kù)進(jìn)行科學(xué)管理，同時(shí)也實(shí)現(xiàn)煤炭低碳開采。而適宜于煤炭開采沉陷區(qū)的區(qū)域土壤有機(jī)碳含量空間預(yù)測(cè)模型是必須要解決的科學(xué)問題。但是由于礦區(qū)受到人類影響較大，區(qū)域內(nèi)部的土壤都不同程度地出現(xiàn)了地表破壞、土壤侵蝕、地面沉陷、植被毀損等問題[9-13]，無論哪一方面的變化都會(huì)對(duì)土壤碳庫(kù)產(chǎn)生影響，使礦區(qū)土壤有機(jī)碳庫(kù)通常存在比較強(qiáng)的空間變異性[14-19]，因此，有必要結(jié)合煤炭開采活動(dòng)引起的擾動(dòng)影響因素對(duì)礦區(qū)土壤有機(jī)碳庫(kù)進(jìn)行估算，以更好地了解煤炭開采活動(dòng)帶來的區(qū)域尺度上土壤有機(jī)碳庫(kù)含量的變化。

當(dāng)前關(guān)于區(qū)域土壤含量預(yù)測(cè)主要有 2種思路，一種是基于遙感數(shù)據(jù)對(duì)區(qū)域土壤有機(jī)碳庫(kù)預(yù)測(cè)，另一種是運(yùn)用地統(tǒng)計(jì)學(xué)方法進(jìn)行預(yù)測(cè)。但目前，基于遙感的方法存在著一個(gè)主要的缺陷，即尚未建立航天遙感與土壤有機(jī)碳庫(kù)之間的定量關(guān)系，因此地統(tǒng)計(jì)學(xué)思路就成為預(yù)測(cè)區(qū)域土壤有機(jī)碳庫(kù)的主要方法。

Kriging空間插值法是應(yīng)用最廣泛的土壤屬性空間預(yù)測(cè)方法[20-23]。但是普通Kriging比較適用于地質(zhì)及人文環(huán)境等各方面條件相對(duì)比較均質(zhì)的區(qū)域，在空間區(qū)域較小但地形復(fù)雜，土壤屬性受人類活動(dòng)擾動(dòng)、變化特別強(qiáng)烈的地區(qū)，其應(yīng)用效果并不理想。

隨著計(jì)算機(jī)技術(shù)的不斷進(jìn)步，結(jié)合輔助信息的Kriging方法已經(jīng)得到了廣泛的應(yīng)用，該方法可以提高預(yù)測(cè)精度。已有研究表明，許多變量可以作為Kriging插值法的輔助變量，提高土壤有機(jī)碳庫(kù)空間預(yù)測(cè)的精度[24]。例如Chai等[25]以北京平谷區(qū)為研究區(qū)，分別利用基于約束性極大似然模型的最佳線性無偏估計(jì)法和回歸Kriging法，對(duì)區(qū)域土壤有機(jī)質(zhì)含量進(jìn)行了空間預(yù)測(cè)，并比較了二者的預(yù)測(cè)精度；Mishra 等[20]以印第安納州為研究區(qū)，分別使用剖面深度分布函數(shù)和普通Kriging模型，預(yù)測(cè)了表層土壤有機(jī)碳的含量，并比較了 2種預(yù)測(cè)模型的預(yù)測(cè)精度。Aladamat R.等[26]建立了全球土壤有機(jī)碳建模系統(tǒng)，并結(jié)合研究區(qū)土地利用的變化情況，預(yù)測(cè)了約旦全國(guó)尺度上2000年到2030年的土壤有機(jī)碳的儲(chǔ)量。

本文提出在煤炭開采劇烈擾動(dòng)的開采沉陷區(qū)，利用結(jié)合輔助變量的 Kriging法來預(yù)測(cè)土壤有機(jī)碳含量的思路。在高潛水位煤礦區(qū)，煤炭開采沉陷形成了不同的積水區(qū)，根據(jù)不同的積水情況人們采取了不同的土地利用措施，導(dǎo)致不同區(qū)域之間土壤有機(jī)碳含量呈現(xiàn)出很大的差異。本文以徐州九里煤炭開采沉陷區(qū)為研究區(qū)，采用分區(qū)Kriging法，以沉陷積水情況為分區(qū)輔助變量，進(jìn)行了土壤有機(jī)碳含量預(yù)測(cè)；同時(shí)也利用普通Kriging法對(duì)沉陷區(qū)土壤有機(jī)碳庫(kù)含量進(jìn)行了預(yù)測(cè)。通過對(duì) 2種方法的預(yù)測(cè)精度的比較，選出精度更好的礦區(qū)土壤有機(jī)碳庫(kù)儲(chǔ)量估算模型，作為礦區(qū)煤炭開采沉陷區(qū)內(nèi)土壤有機(jī)碳庫(kù)含量空間預(yù)測(cè)的方法。

1 材料與方法

1.1 研究區(qū)概況

九里采煤塌陷區(qū)（34°13'39''–34°26'16''N, 117°06'21''–117°12'16''E）位于江蘇省徐州市九里區(qū)，面積約為42.15 km2，是徐州市張小樓煤礦、龐莊煤礦、夾河煤礦的煤炭開采沉陷區(qū)。3礦區(qū)的井田邊界以及沉陷區(qū)域如圖1所示（紅線范圍內(nèi)）。研究區(qū)煤礦的開采方式為井工開采，選取長(zhǎng)壁采煤的技術(shù)。這種采煤技術(shù)造成了地面塌陷，并伴隨著一系列的生態(tài)環(huán)境問題，主要包括：（1）煤炭開采導(dǎo)致農(nóng)田大面積塌陷；（2）煤炭開采活動(dòng)導(dǎo)致一些位于高潛水位的農(nóng)田變成了沉陷積水水域。

圖1 研究區(qū)位置以及采樣點(diǎn)分布圖Fig.1 Location of study area and distribution of soil sampling sites

研究區(qū)屬于中國(guó)華東平原區(qū)，平均海拔在 35～44 m之間，地下水位較淺，屬于溫暖季風(fēng)氣候區(qū)，年平均氣溫14 ℃，年平均無霜期為200～220 d，年平均降水量為800～-930 mm。該區(qū)是古黃河沖積平原，水稻土、砂姜土為該區(qū)主要土壤類型。農(nóng)業(yè)生產(chǎn)按照冬小麥/夏玉米的方式進(jìn)行輪作，屬于中國(guó)東部的糧食生產(chǎn)區(qū)。

1.2 土壤樣品采集與處理

采樣點(diǎn)包括預(yù)測(cè)樣點(diǎn)和驗(yàn)證樣點(diǎn) 2部分，預(yù)測(cè)樣點(diǎn)采樣采用的方法是正方形格網(wǎng)法。根據(jù)研究區(qū)的實(shí)際面積，格網(wǎng)大小為1000 m 1000 m×，共采集了54個(gè)預(yù)測(cè)樣點(diǎn)。本研究只采用了每個(gè)土樣的0～20 cm土壤有機(jī)碳含量數(shù)據(jù)。在采樣過程中，用手持GPS機(jī)記錄下每一個(gè)采樣點(diǎn)的地理坐標(biāo)，其中沉陷濕地區(qū)內(nèi)采集8個(gè)土壤樣點(diǎn)，季節(jié)性積水區(qū)內(nèi)采集13個(gè)土壤樣點(diǎn)，沉陷未積水區(qū)內(nèi)采集33個(gè)土壤樣點(diǎn)。另外，在整個(gè)沉陷區(qū)內(nèi)再隨機(jī)采集18個(gè)樣點(diǎn)作為驗(yàn)證樣點(diǎn)，與 2種估算模型的預(yù)測(cè)值對(duì)比，對(duì)模型所得出的結(jié)果進(jìn)行精度評(píng)價(jià)。外業(yè)采樣完成后對(duì)土壤有機(jī)碳含量采用重鉻酸鉀氧化法進(jìn)行測(cè)定[27]。本次采樣工作是在2010年10月份研究區(qū)域內(nèi)所有的農(nóng)作物收割完成以后進(jìn)行的。

1.3 煤炭開采沉陷區(qū)土壤有機(jī)碳含量空間預(yù)測(cè)方法

對(duì)區(qū)域土壤理化性質(zhì)進(jìn)行空間預(yù)測(cè)常用的方法是Kriging空間插值法。它是基于變異函數(shù)模型在有限區(qū)域內(nèi)對(duì)區(qū)域化變量的取值進(jìn)行無偏最優(yōu)估計(jì)的空間插值方法。其基本原理和方法在許多文獻(xiàn)中均有詳細(xì)描述[23,28-29]，是基于區(qū)域化變量理論，利用原始數(shù)據(jù)和半方差函數(shù)的結(jié)構(gòu)性，計(jì)算可獲取變量的線性加權(quán)組合，從而對(duì)待估點(diǎn)預(yù)測(cè)值進(jìn)行無偏的最佳估計(jì)，而且這種方法能夠計(jì)算估計(jì)值的可靠性程度（估計(jì)值方差）。

研究區(qū)包括多個(gè)煤礦，且處于高潛水位區(qū)。由于每個(gè)礦區(qū)對(duì)煤炭資源的開采程度不同，導(dǎo)致各個(gè)礦區(qū)內(nèi)部的地表沉陷程度不同，進(jìn)而形成了不同程度的沉陷積水區(qū)（包括沉陷未積水區(qū)、季節(jié)性積水區(qū)和常年積水濕地區(qū)），人們針對(duì)礦區(qū)不同的積水情況采取了不同的土地利用方式，以及農(nóng)業(yè)耕作措施，因而土壤有機(jī)碳含量受到影響也會(huì)產(chǎn)生一定的差異[30]。在高潛水位礦區(qū)由開采沉陷導(dǎo)致的沉陷積水情況，是影響沉陷區(qū)土壤有機(jī)碳庫(kù)的主導(dǎo)因素，因此根據(jù)不同的積水情況來分區(qū)進(jìn)行Kriging插值，將樣本中土壤有機(jī)碳庫(kù)的含量分為 2部分，即各個(gè)積水區(qū)的均值和殘差。其計(jì)算公式如下：

式中 Z （ Xkj）表示各個(gè)采樣點(diǎn)的有機(jī)碳含量， Xkj表示的是土樣所在的位置， u（ tk）表示各個(gè)積水區(qū)有機(jī)碳含量的均值， tk表示的是土樣所在沉陷區(qū)的積水情況：主要分為常年積水濕地區(qū)、季節(jié)性積水區(qū)以及未積水區(qū) 3個(gè)部分。其中，k表示不同的區(qū)域積水情況，j表示不同的土壤樣點(diǎn)。r(xkj)表示殘差，即各個(gè)采樣點(diǎn)的有機(jī)碳減去有機(jī)碳均值的差。其中，均值間的差異反映的是沉陷積水情況不同區(qū)域中有機(jī)碳含量的變異性，而殘差之間的差異反映的是沉陷積水情況相同的地區(qū)有機(jī)碳含量的變異性。在這種方法中，將殘差看作是新的區(qū)域變量來進(jìn)行Kriging插值，其變異函數(shù) γr（ h）及待估點(diǎn)Xkj的預(yù)測(cè)公式分別為式（2）和式（3）。每一個(gè)待估點(diǎn)的土壤有機(jī)碳庫(kù)含量預(yù)測(cè)值 Z*（xkj）為所屬沉陷積水情況區(qū)的均值u(tk)與

所估計(jì)的殘差r*（xkj）之和（式4）。

式中h為2個(gè)樣本間的點(diǎn)空間分隔距離，N(h)是分隔距離為h時(shí)的樣本點(diǎn)對(duì)總數(shù)。

式中kjλ為權(quán)重系數(shù)。

在本次研究中，首先用SPSS來對(duì)土壤有機(jī)碳庫(kù)數(shù)據(jù)及其殘差進(jìn)行統(tǒng)計(jì)分析，然后利用ArcGIS 10.0軟件中的Geostatistical Analyst 模塊，檢驗(yàn)數(shù)據(jù)的質(zhì)量和分布，本次的采樣數(shù)據(jù)滿足正態(tài)分布和隨機(jī)性檢驗(yàn)，用該軟件進(jìn)行空間插值并生成空間插值圖。

1.4 預(yù)測(cè)結(jié)果準(zhǔn)確性檢驗(yàn)

本次對(duì)于研究結(jié)果準(zhǔn)確性的檢驗(yàn)是通過比較18個(gè)驗(yàn)證樣點(diǎn)位置處的預(yù)測(cè)值與采樣實(shí)測(cè)值之間的相關(guān)系數(shù)（r）[31]，以及Isaaks和Shrivastava定義的3個(gè)驗(yàn)證指標(biāo)，即平均誤差（mean error，ME）、平均絕對(duì)誤差（mean absolute error，MAE）和均方根誤差（root mean square error，RMSE）來評(píng)估不同插值方法的性能[32]。其中，相關(guān)系數(shù)（r）這一指標(biāo)是用以反映變量之間相關(guān)關(guān)系密切程度的。均方根誤差RMSE實(shí)際上是觀測(cè)值與真值偏差的平方和再與觀測(cè)次數(shù)N來求比值，然后再求平方根來計(jì)算得到的。

式（5）～（7）中，N表示驗(yàn)證點(diǎn)的總數(shù)量，本次研究中N為18，xoi表示驗(yàn)證點(diǎn)的實(shí)測(cè)值，xpi表示根據(jù)估算模型得到的驗(yàn)證點(diǎn)預(yù)測(cè)值。

在實(shí)際的測(cè)量過程中，由于對(duì)樣點(diǎn)進(jìn)行的觀測(cè)次數(shù)總是有限的，因此只能用最可信賴的值來代替真值，而對(duì)一組測(cè)量中的特大或特小誤差，標(biāo)準(zhǔn)誤差的反映非常敏感，所以，可以用標(biāo)準(zhǔn)誤差來很好地反映出預(yù)測(cè)的精確度。r值越大、RMSE值越小，說明所預(yù)測(cè)精度越高, 反之則說明預(yù)測(cè)精度越低。

2 結(jié)果與分析

2.1 土壤有機(jī)碳的統(tǒng)計(jì)分析

首先對(duì)試驗(yàn)區(qū)所采集的土壤樣本進(jìn)行了描述性統(tǒng)計(jì)分析（表 1）。分析顯示，區(qū)域內(nèi)采樣點(diǎn)的土壤有機(jī)碳含量符合正態(tài)分布規(guī)律，不同積水區(qū)內(nèi)土壤有機(jī)碳庫(kù)含量的差異較為顯著，因此能夠進(jìn)行Kriging插值；同時(shí)也為下一步對(duì)土壤有機(jī)碳庫(kù)含量進(jìn)行半方差函數(shù)結(jié)構(gòu)分析以及相應(yīng)的分區(qū)Kriging插值方案的確定提供了依據(jù)，并且減少了統(tǒng)計(jì)分析過程中所產(chǎn)生的不確定性。根據(jù)表 2可知，3個(gè)區(qū)域內(nèi)土壤有機(jī)碳含量從小到大依次為未積水區(qū)、季節(jié)性積水區(qū)、濕地，這說明土壤有機(jī)碳庫(kù)的含量在不同的積水區(qū)域內(nèi)有很大的不同。而且積水區(qū)內(nèi)土壤有機(jī)碳庫(kù)的變異系數(shù)為 0.07，比未積水區(qū)內(nèi)土壤有機(jī)碳庫(kù)的變異系數(shù)0.13要小，尤其在季節(jié)性積水區(qū)內(nèi)其標(biāo)準(zhǔn)差為 0.86，也較小，說明積水區(qū)內(nèi)的土壤有機(jī)碳庫(kù)含量離散程度較小，比較集中，這也從另一方面驗(yàn)證了將沉陷區(qū)內(nèi)積水量的變化作為輔助變量的分區(qū)Kriging插值法的可行性。

表1 沉陷區(qū)域外業(yè)采樣樣本土壤有機(jī)碳庫(kù)含量描述性統(tǒng)計(jì)結(jié)果Table 1 Descriptive statistical results of soil organic carbon content of samples in subsidence area

2.2 地統(tǒng)計(jì)分析及Kriging插值

利用ArcGIS 10.0，對(duì)采樣點(diǎn)樣本土壤有機(jī)碳含量的原始數(shù)據(jù)、去除區(qū)域積水分區(qū)內(nèi)有機(jī)碳含量均值以后的殘差半方差函數(shù)及其擬合參數(shù)進(jìn)行分析[33-36]，其結(jié)果見表3。

表2 沉陷區(qū)不同積水情況下土壤樣本有機(jī)碳描述性統(tǒng)計(jì)結(jié)果Table 2 Descriptive statistics results of soil organic carbon samples with different water accumulation in subsidence area

表3 土壤有機(jī)碳含量及其殘差半方差函數(shù)模型和參數(shù)Table 3 Soil organic carbon content and its residual semi variance function model and parameters

表 3顯示了樣本有機(jī)碳含量原始數(shù)據(jù)和土壤有機(jī)碳?xì)埐罡髯缘陌敕讲詈瘮?shù)擬合指數(shù)和球狀模型的參數(shù)。根據(jù)塊金值的相關(guān)定義，理論上當(dāng)采樣點(diǎn)的距離為 0時(shí)，半變異函數(shù)值應(yīng)為 0，但由于存在測(cè)量誤差以及空間變異，使得 2個(gè)采樣點(diǎn)非常接近時(shí)，它們的半變異函數(shù)值不為0，即存在塊金值，因此塊金值反映的是最小抽樣尺度以下變量的變異性及測(cè)量誤差。結(jié)合對(duì)表 3的分析，可以看出采樣點(diǎn)土壤有機(jī)碳含量與殘差的塊金值分別為0.91和 0.01，土壤有機(jī)碳的空間異質(zhì)性較大，殘差的空間異質(zhì)性比較小，說明該地區(qū)不同積水沉陷區(qū)域內(nèi)部由于試驗(yàn)誤差和一些隨機(jī)因素引起的空間變異都很小。

C是結(jié)構(gòu)方差，表示系統(tǒng)屬性或區(qū)域變量最大空間變異，是由土壤母質(zhì)、地形、氣候等結(jié)構(gòu)性因素引起的變異。由于該沉陷區(qū)域開采前后氣候并未發(fā)生變化，土壤有機(jī)碳庫(kù)的空間變異性基本是由于煤炭開采造成的地表形變和沉陷積水等結(jié)構(gòu)因素引起的變異，該沉陷區(qū)域土壤有機(jī)碳結(jié)構(gòu)方差C是15.17，這種結(jié)構(gòu)性因素是該區(qū)域土壤有機(jī)碳庫(kù)變異的主導(dǎo)因素。

C/(C0+C)即塊金值與基臺(tái)值的比值，被稱為空間相關(guān)度，表示可度量空間自相關(guān)的變異所占的比例，該值的大小表示系統(tǒng)變量的空間相關(guān)性的程度。由表3可知2種方法得出的 C/(C0+C)分別為 0.998與 0.943，均大于75%，這說明其相關(guān)性很強(qiáng)，尤其是樣本的土壤有機(jī)碳含量的相關(guān)性極強(qiáng)，達(dá)到了0.998。

根據(jù)表 3可知，與原始數(shù)據(jù)相比，土壤有機(jī)碳?xì)埐顢M合函數(shù)的塊金值很小，基臺(tái)值降低，變程也降低，這主要是因?yàn)?，殘差是土壤有機(jī)碳含量減去各個(gè)積水區(qū)內(nèi)的有機(jī)碳含量均值得到的，而不同積水區(qū)域內(nèi)部的土壤有機(jī)碳含量的變化比較小，相對(duì)比較均勻，這也證明了對(duì)沉陷區(qū)內(nèi)土壤有機(jī)碳含量進(jìn)行預(yù)測(cè)時(shí)考慮不同積水沉陷區(qū)的合理性，間接地說明了以分區(qū)對(duì)沉陷區(qū)土壤碳庫(kù)儲(chǔ)量進(jìn)行Kriging插值預(yù)測(cè)的有效性和可行性。

2.3 兩種方法預(yù)測(cè)精度的對(duì)比分析

采用研究區(qū)18個(gè)驗(yàn)證點(diǎn)處的土壤有機(jī)碳含量實(shí)測(cè)值與模型的預(yù)測(cè)值進(jìn)行對(duì)比，對(duì)分區(qū)Kriging模型的預(yù)測(cè)精度進(jìn)行評(píng)價(jià)。圖2是預(yù)測(cè)值與實(shí)測(cè)值的散點(diǎn)分布圖，由2種方法所得到的線性回歸方程顯示，以區(qū)域積水情況為輔助變量進(jìn)行的Kriging插值得到的預(yù)測(cè)值與實(shí)測(cè)值的決定系數(shù)（0.756 4）明顯高于普通的Kriging插值法（0.508 6）。

表4為2種預(yù)測(cè)方法的驗(yàn)證點(diǎn)處土壤有機(jī)碳含量預(yù)測(cè)的準(zhǔn)確性評(píng)估指標(biāo)。ME表示預(yù)測(cè)的平均偏差，從表4中可以看出，普通Kriging（0.020 2 g/kg）和結(jié)合積水情況的Kriging（0.021 8 g/kg）的ME均接近0，說明2種方法預(yù)測(cè)的平均偏差都比較小，總體預(yù)測(cè)精度都比較高。平均絕對(duì)誤差（MEA）反映了預(yù)測(cè)值誤差的實(shí)際情況，由表4可知，普通Kriging的MAE為1.851 1 g/kg，結(jié)合積水情況的Kriging的MAE為1.287 8 g/kg，因此，結(jié)合積水情況的Kriging預(yù)測(cè)實(shí)際誤差更小。作為平均值和方差偏差的聯(lián)合度量指標(biāo)，如表4所示，直接進(jìn)行Kriging插值對(duì)土壤有機(jī)碳含量的預(yù)測(cè)RMSE為0.55，明顯高于分區(qū)Kriging插值（0.35）。綜上所述，考慮區(qū)域內(nèi)土壤積水情況的Kriging插值法對(duì)土壤有機(jī)碳含量的預(yù)測(cè)值精度更高。

圖2 兩種方法得到的土壤有機(jī)碳預(yù)測(cè)值與實(shí)測(cè)值散點(diǎn)圖Fig.2 Scatter diagram of predicted and measured values of soil organic carbon by 2 methods

表4 兩種方法預(yù)測(cè)土壤有機(jī)碳含量的準(zhǔn)確性評(píng)估指標(biāo)Table 4 Accuracy assessment indices of estimated soil organic carbon content with 2 methods

2.4 土壤有機(jī)碳的空間分布特征

通過 2種方法預(yù)測(cè)得到的土壤碳庫(kù)儲(chǔ)量空間分布如圖3所示，可以發(fā)現(xiàn)沉陷區(qū)2種土壤有機(jī)碳含量空間預(yù)測(cè)方法（普通Kriging和分區(qū)Kriging）得到的土壤有機(jī)碳?xì)埐畹目臻g分布格局，而分區(qū)Kriging預(yù)測(cè)的土壤有機(jī)碳含量圖的空間分布特征及空間遞變規(guī)律更加明顯，普通Kriging得到的土壤有機(jī)碳含量范圍為：9.34～16.252 g/kg，而分區(qū)Kriging插值的結(jié)果為9.338～18.058 g/kg，2種方法的估計(jì)范圍大致相同，總體上看研究區(qū)土壤有機(jī)碳分布的空間格局基本一致，有機(jī)碳含量由中部向四周逐漸減少，中部的土壤有機(jī)碳含量最高，據(jù)實(shí)地調(diào)查該區(qū)域?qū)儆跐竦貐^(qū)，與其他土地利用相比，濕地土壤長(zhǎng)期飽和度較低，其凋落物和有機(jī)質(zhì)分解速度較慢，這增加了濕地土壤中的土壤有機(jī)碳含量[37]。

圖3 兩種方法得到的土壤有機(jī)碳含量空間分布Fig.3 Spatial distribution of soil organic carbon content by 2 methods

西南部的土壤有機(jī)碳含量較少，該區(qū)域?qū)儆谖捶e水區(qū)，由于礦區(qū)煤炭開采對(duì)農(nóng)田土壤和植被擾動(dòng)影響嚴(yán)重，造成植被NPP (net primary productivity) 降低，從而減少了土壤有機(jī)碳的一個(gè)重要補(bǔ)給來源，而且開采沉陷形成很多沉陷坡面，土壤有機(jī)碳由于土壤侵蝕的原因流失嚴(yán)重[38-39]。從整體上來說 2種空間插值方法都能夠反映出整個(gè)研究區(qū)的土壤有機(jī)碳分布情況。

但是在局部分布上，2種預(yù)測(cè)方法之間存在著一定的差異。對(duì)比分析發(fā)現(xiàn)，直接進(jìn)行Kriging插值得到的土壤有機(jī)碳空間分布的圖斑比較連續(xù)規(guī)整，這明顯跟研究區(qū)的實(shí)際情況不太相符，說明平滑效應(yīng)使不同積水情況區(qū)域之間的差異性降低，其差值結(jié)果只能反映研究區(qū)土壤有機(jī)碳的大致分布情況，卻不能更為精確地反映研究區(qū)的實(shí)際情況，比如圖3a局部S1區(qū)以及S2區(qū)域由于個(gè)別預(yù)測(cè)樣點(diǎn)的土壤有機(jī)碳含量偏高，造成區(qū)域內(nèi)土壤有機(jī)碳含量空間插值結(jié)果異常偏高。普通Kriging插值是基于采樣點(diǎn)的空間分布特征以及樣點(diǎn)與樣點(diǎn)之間的空間位置關(guān)系來進(jìn)行空間插值，事先未對(duì)插值樣點(diǎn)進(jìn)行分類，所以無法消除鄰近的不同類別空間預(yù)測(cè)樣點(diǎn)對(duì)空間插值結(jié)果的影響。而分區(qū)Kriging由于按照影響沉陷區(qū)土壤有機(jī)碳含量的最大因素——沉陷積水情況進(jìn)行分區(qū)，就可以有效消除這方面的誤差影響。而且，與分區(qū)Kriging的模擬結(jié)果相比，直接Kriging法模擬的礦區(qū)土壤有機(jī)碳分布特征規(guī)律性不明顯，忽略了土壤有機(jī)碳空間平穩(wěn)過渡。因此，分區(qū)Kriging模擬更能反應(yīng)土壤有機(jī)碳的空間遞變特點(diǎn)，有利于分析不同因素對(duì)土壤有機(jī)碳空間分布的影響，能更好地反映研究區(qū)域內(nèi)部土壤有機(jī)碳含量的分布情況。

分區(qū)Kriging空間插值的思想在相關(guān)領(lǐng)域也得到了學(xué)者的認(rèn)可，Zhang等[31]通過對(duì)比在土壤現(xiàn)場(chǎng)采樣過程中分配樣點(diǎn)的 4種不同模式，即不分類的網(wǎng)格模式，基于土壤類型的分配模式，基于土地利用類型的分配模式以及基于土地利用模式-土壤類型的分配模式下得到的土壤有機(jī)碳的變異系數(shù)，得出評(píng)估土壤有機(jī)碳空間分布的最有效方法是基于土地利用模式-土壤類型分配模式。同樣地，Sandeep等[40]在研究美國(guó)賓夕法尼亞州的土壤有機(jī)碳時(shí)，采用結(jié)合環(huán)境變量的地理加權(quán)回歸克里金法，與原來的回歸克里金法進(jìn)行對(duì)比發(fā)現(xiàn)，因?yàn)榍罢呖紤]了空間非平穩(wěn)性以及殘差的空間自相關(guān)性，其精度會(huì)提高。所以在受人類煤炭開采活動(dòng)擾動(dòng)劇烈的沉陷區(qū)土壤有機(jī)碳含量的預(yù)測(cè)應(yīng)用中，其預(yù)測(cè)精度也高于回歸克里金法。

3 結(jié) 論

本文提出了利用分區(qū) Kriging法來對(duì)煤炭開采沉陷區(qū)土壤有機(jī)碳含量進(jìn)行空間插值預(yù)測(cè)的方法，并將該方法與傳統(tǒng)的直接Kriging方法得到的預(yù)測(cè)結(jié)果進(jìn)行對(duì)比，發(fā)現(xiàn)：

1）結(jié)合區(qū)域內(nèi)部積水情況來進(jìn)行的分區(qū) Kriging法得到的研究區(qū)土壤有機(jī)碳含量的范圍為：9.338～18.058 g/kg，而直接Kriging方法得到的范圍為：9.34～16.252 g/kg。

2）從區(qū)域預(yù)測(cè)精度上分析，與直接進(jìn)行 Kriging插值相比，分區(qū)Kriging預(yù)測(cè)方法精度更高。研究區(qū)內(nèi)由于煤炭開采造成了不同的沉陷積水區(qū)，土地利用方式也各不相同，導(dǎo)致區(qū)域內(nèi)土壤有機(jī)碳含量出現(xiàn)差異。結(jié)合區(qū)域積水情況的分區(qū)Kriging，可以消除鄰近不同積水情況區(qū)的采樣點(diǎn)土壤有機(jī)碳含量對(duì)其空間插值預(yù)測(cè)精度的影響，提高土壤有機(jī)碳含量的預(yù)測(cè)精度。

綜上所述，可以選擇結(jié)合區(qū)域積水情況的分區(qū)Kriging空間插值作為煤炭開采沉陷區(qū)土壤有機(jī)碳含量空間預(yù)測(cè)的預(yù)測(cè)方法。在煤礦區(qū)，優(yōu)化的土壤有機(jī)碳空間預(yù)測(cè)方法，也有利于更精確地進(jìn)行土壤有機(jī)碳動(dòng)態(tài)模擬和理解本地土壤有機(jī)碳庫(kù)的時(shí)空演化，為礦區(qū)土地低碳復(fù)墾乃至區(qū)域內(nèi)的土地資源的低碳利用提供更加科學(xué)的依據(jù)。

[1] BP. The BP Statistical Review of World Energy 2014[EB/OL]. http://www.bp.com/statisticalreview, 2018-03-12.

[2] 原野, 趙中秋, 白中科, 等. 露天煤礦復(fù)墾生物多樣性恢復(fù)技術(shù)體系與方法：以平朔礦排土場(chǎng)為例[J]. 中國(guó)礦業(yè),2017, 26(8): 93－98.Yuan Ye, Zhao Zhongqiu, Bai Zhongke, et al. Technology system and method of biodiversity restoration for the reclamation of opencast coal mine:a case study from the dumps in Pingshuo mine[J]. China Mining Magazine, 2017,26(8): 93－98. (in Chinese with English abstract)

[3] 羅愛武. 淮北市采煤塌陷區(qū)土地復(fù)墾研究[J]. 安徽師范大學(xué)學(xué)報(bào)（自然科學(xué)版）, 2002, 25(3): 286－289.Luo Aiwu. The land restoration in the sunk areas of coal extraction in HuaiBei City[J]. Journal of AnHui Normal University (Nature Science), 2002, 25(3): 286－289. (in Chinese with English abstract)

[4] 彭蘇萍, 趙建慶. 中國(guó)西部煤礦區(qū)生態(tài)環(huán)境控制及改善研究[C]//中國(guó)科協(xié)2000年學(xué)術(shù)年會(huì). 2000.

[5] 侯湖平, 徐占軍, 張紹良, 等. 煤炭開采對(duì)區(qū)域農(nóng)田植被碳庫(kù)儲(chǔ)量的影響評(píng)價(jià)[J]. 農(nóng)業(yè)工程學(xué)報(bào), 2014, 30(5): 1－9.Hou Huping, Xu Zhanjun, Zhang Shaoliang, et al. Effect evaluation on vegetation carbon pool of region agro-ecosystem by coal mining in mining area[J].Transactions of the Chinese Society of Agricultural Engineering (Transactions of the CSAE), 2014, 30(5): 1－9.(in Chinese with English abstract)

[6] 王坤. 高潛水位采煤塌陷區(qū)充填復(fù)墾土壤碳動(dòng)態(tài)研究[D].徐州：中國(guó)礦業(yè)大學(xué), 2016.Wang kun. Study on Carbon Dynamics of Reclaimed Soil in Coal Mining Subsidence Area with High Groundwater Level[D]. Xuzhou: China University of Mining and Technology, 2016.

[7] Miller B A, Koszinski S, Wehrhan M, et al. Comparison of spatial association approaches for landscape mapping of soil organic carbon stocks[J]. Soil, 2015(1): 217－233.

[8] Zhang G S, Ni Z W. Winter tillage impacts on soil organic carbon, aggregation and CO2emission in a rainfed vegetable cropping system of the mid–Yunnan plateau, China[J]. Soil& Tillage Research, 2017, 165: 294－301.

[9] Lan C Y, Shu W S, Wong M H. Revegetation of lead/zinc mine tailings at shaoguan, guandong province, China:Phytotoxicity of the tailings[J]. Studies in Environmental Science, 1997, 66(97): 119－130.

[10] Yan Chu. The influence of coal mining on environmental quality of mining area and the research progress of repairing technology[J]. Advances in Environmental Protection, 2016,6(1): 1－6.

[11] Anthony B. Restoration of mined land-using natural processes[J]. Ecological Engineering, 1997, 8(4): 255－269.

[12] Zhang Ling, Wang Jinman, Bai Zhongke, et al. Effects of vegetation on runoff and soil erosion on reclaimed land in an opencast coal-mine dump in a loess area[J]. Catena, 2015,128(5): 44－53.

[13] Huang Yi, Tian Feng, Wang Yunjia, et al. Effect of coal mining on vegetation disturbance and associated carbon loss[J]. Environmental Earth Sciences, 2013, 73(5): 1－14.

[14] álvaro-Fuentes J, Easter M, Paustian K. Climate change effects on organic carbon storage in agricultural soils of northeastern Spain[J]. Agriculture Ecosystems &Environment, 2012, 155(155): 87－94.

[15] Yang X M, Drury C F, Wander M M, et al. Evaluating the effect of tillage on carbon sequestration using the minimum detectable difference concept[J]. Pedosphere, 2008, 18(4):421－430.

[16] Mishra U, Ussiri David A N, Lal R. Tillage effects on soil organic carbon storage and dynamics in Corn Belt of Ohio USA[J]. Soil & Tillage Research, 2010, 107(2): 88－96.

[17] Rosemary F, Vitharana U W A, Indraratne S P, et al.Exploring the spatial variability of soil properties in an Alfisol soil catena[J]. Catena, 2017, 150(3): 53－61.

[18] Teng Mingjun, Zeng Lixiong, Xiao Wenfa, et al. Spatial variability of soil organic carbon in Three Gorges Reservoir area, China[J]. Science of the Total Environment, 2017,599(6): 1308－1316.

[19] Wu Lizhi, Li Long, Yao Yunfeng, et al. Spatial distribution of soil organic carbon and its influencing factors at different soil depths in a semiarid region of China[J]. Environmental Earth Sciences, 2017, 76 (19): 654.

[20] Mishra U, Lal R, Slater B, et al. Predicting soil organic carbon stock using profile depth distribution functions and ordinary kriging[J]. Soil Science Society of America Journal,2009, 73(2): 614－621.

[21] Liu Yaolin, Guo Long, Jiang Qinghu, et al. Comparing geospatial techniques to predict SOC stocks[J]. Soil &Tillage Research, 2015, 148(3): 46－58.

[22] 杜挺, 楊聯(lián)安, 張泉, 等. 縣域土壤養(yǎng)分協(xié)同克里格和普通克里格空間插值預(yù)測(cè)比較——以陜西省藍(lán)田縣為例[J].陜西師范大學(xué)學(xué)報(bào)(自然科學(xué)版), 2013, 41(4): 85－89.Du Ting, Yang Lian'an, Zhang Quan, et al. Spatial prediction comparison of soil nutrient between ordinary Kriging and Cokring at county scale-A case study in Lantian county of Shaanxi province[J]. Journal of Shaanxi Normal University(Natural Science Edition), 2013, 41(4): 85－89. (in Chinese with English abstract)

[23] 施周, 閆杭召, 畢晨, 等. 基于地統(tǒng)計(jì)學(xué)——克里格插值法的村鎮(zhèn)地表水體水質(zhì)監(jiān)測(cè)[J]. 環(huán)境工程學(xué)報(bào), 2017(4):2607－2613.Shi Zhou, Yan Hangzhao, Bi Chen, et al. Surface water quality monitoring in town based on Geotatistics-Kriging interpolation[J]. Chinese Journal of Environmental Engineering, 2017(4): 2607－2613. (in Chinese with English abstract)

[24] Dai Fuqiang, Zhou Qigang, Lv Zhiqiang, et al. Spatial prediction of soil organic matter content integrating artificial neural network and ordinary kriging in Tibetan Plateau[J].Ecological Indicators, 2014, 45(5): 184－194.

[25] Chai Xurong, Shen Chongyang, Yuan Xiaoyong, et al.Spatial prediction of soil organic matter in the presence of different external trends with REML-EBLUP[J]. Geoderma,2008, 148(2): 159－166.

[26] Aladamat R, Rawajfih Z, Easter M, et al. Predicted soil organic carbon stocks and changes in Jordan between 2000 and 2030 made using the GEF SOC Modelling System[J]. Agriculture Ecosystems & Environment, 2007, 122(1): 35－45.

[27] Fang Yu, Yan Zhilei, Chen Jichen, et al. Effect of chemical fertilization and green manure on the abundance and community structure of ammonia oxidizers in a paddy soil[J]. Chilean Journal of Agricultural Research, 2015, 75(4): 487－496.

[28] Wang L X, Okin G S, Caylor K K, et al. Spatial heterogeneity and sources of soil carbon in southern African savannas[J]. Geoderma, 2009, 149(3-4): 402－408.

[29] Sun Zhili, Wang Jian, Li Rui, et al. LIF: A new Kriging based learning function and its application to structural reliability analysis[J]. Reliability Engineering & System Safety, 2017, 157(4): 152－165.

[30] Shrestha R K, Lal R. Changes in physical and chemical properties of soil after surface mining and reclamation[J].Geoderma, 2011, 161(3–4): 168－176.

[31] Zhang Z Q, Yu D S, Shi X Z, et al. Effect of sampling classification patterns on SOC variability in the red soil region, China[J]. Soil & Tillage Research, 2010, 110(1): 2－7.

[32] Isaaks E H, Srivastava R M. An Introduction to Applied Geostatistics[M]. Oxford University Press, 1989, 33(33): 483－485.

[33] Lévesque J, King D J. Airborne digital camera image semivariance for evaluation of forest structural damage at an acid mine site[J]. Remote Sensing of Environment, 1999,68(2): 112－124.

[34] Li Xinrong. Influence of variation of soil spatial heterogeneity on vegetation restoration[J]. China Science:Earth Science, 2005, 48(11): 2020－2031.

[35] 郭凌俐, 王金滿, 白中科, 等. 黃土區(qū)露天煤礦排土場(chǎng)復(fù)墾初期土壤顆粒組成空間變異分析[J]. 中國(guó)礦業(yè), 2015,24(2): 52－59.Guo Lingli, Wang Jinman, Bai Zhongke, et al. Analysis of spatial variability of soil granules in early stage of reclamation at opencast coal minedump in loess area[J].China Mining Magazine, 2015, 24(2): 52－59. (in Chinese with English abstract)

[36] 代富強(qiáng), 周啟剛, 劉剛才. 基于回歸克里格和遙感的紫色土區(qū)土壤有機(jī)質(zhì)含量空間預(yù)測(cè)[J]. 土壤通報(bào), 2014, 45(3):562－567.Dai Fuqiang, Zhou Qigang, Liu Gangcai. Spatial prediction of soil organic matter contents in a purplish soil region with regression Kriging and remote sensing[J]. Chinese Journal of Soil Science, 2014, 45(3): 562－567. (in Chinese with English abstract)

[37] Qualls R G, Richardson C J. Decomposition of Litter and Peat in the Everglades: The Influence of P Concentrations[M].Springer New York, 2008.

[38] Hu Zhenqi, Xiao Wu. Optimization of concurrent mining and reclamation plans for single coal seam: A case study in northern Anhui, China[J]. Environmental Earth Sciences,2013, 68(5): 1247－1254.

[39] Indorante S J, Jansen I J, Boast C W. Surface mining and reclamation: Initial changes in soil character[J]. Journal of Soil & Water Conservation, 1981, 36(6): 347－351.

[40] Sandeep K, Rattan L, Liu Desheng. A geographically weighted regression kriging approach for mapping soil organic carbon stock[J]. Geoderma, 2012, 189–190 (6): 627－634.