基于地形分區(qū)IDW 的復(fù)墾排土場(chǎng)土壤養(yǎng)分空間分布

趙艷玲，趙彬程，王 鑫，肖 武，劉慧芳

(1.中國(guó)礦業(yè)大學(xué)(北京)土地復(fù)墾與生態(tài)重建研究所，北京 100083；2.浙江大學(xué) 公共管理學(xué)院，浙江 杭州 310058；3.中國(guó)恩菲工程技術(shù)有限公司，北京 100083)

露天礦的開(kāi)采過(guò)程中，表土層的剝離堆存和回填不可避免地產(chǎn)生土壤養(yǎng)分流失，導(dǎo)致復(fù)墾土壤貧瘠，植被恢復(fù)緩慢，使復(fù)墾區(qū)難以及時(shí)覆蓋植被，加重土壤侵蝕和退化，生態(tài)系統(tǒng)循環(huán)難以建立。王金滿等[1]通過(guò)對(duì)露天煤礦排土場(chǎng)的研究認(rèn)為，土壤養(yǎng)分對(duì)植物生長(zhǎng)起到促進(jìn)作用；肖禮等[2]通過(guò)對(duì)露天煤礦排土場(chǎng)土壤與植被間關(guān)系研究發(fā)現(xiàn)，土壤養(yǎng)分對(duì)于植被恢復(fù)有著重要意義；為準(zhǔn)確掌握土壤養(yǎng)分的空間分布，王黨朝等[3]利用傳統(tǒng)克里金插值法對(duì)勝利露天礦進(jìn)行空間插值，發(fā)現(xiàn)排土場(chǎng)邊坡與平臺(tái)土壤的物理性質(zhì)具有明顯的空間差異；郭凌俐等[4]根據(jù)傳統(tǒng)克里金插值法對(duì)排土場(chǎng)土壤進(jìn)行研究，發(fā)現(xiàn)傳統(tǒng)克里金插值法能直觀揭示復(fù)墾區(qū)不同層次土壤顆粒的空間分布特征；王子良[5]通過(guò)對(duì)銅陵礦區(qū)土壤元素的插值，發(fā)現(xiàn)在最優(yōu)情況下，IDW插值方法具有良好穩(wěn)定的插值效果。但由于研究區(qū)不同，無(wú)法對(duì)插值方法做出明確的優(yōu)劣性判斷，同時(shí)，對(duì)于地形復(fù)雜、人類活動(dòng)干擾較多且變化強(qiáng)烈的地區(qū)(如排土場(chǎng))，一些學(xué)者認(rèn)為全局插值方式并不能直觀有效地對(duì)區(qū)域信息進(jìn)行預(yù)測(cè)，而應(yīng)將研究區(qū)域通過(guò)特征變量進(jìn)行分區(qū)，并對(duì)其分別插值，此種方式能夠提高空間預(yù)測(cè)的精度。胡剛等[6]通過(guò)分區(qū)插值方法，發(fā)現(xiàn)不同地形使用不同的插值方法可以有效提高預(yù)測(cè)精度；U.Mishra 等[7]以印第安納州為研究區(qū)，將剖面深度作為輔助變量和傳統(tǒng)克里金插值結(jié)合使用，有效預(yù)測(cè)了表層土壤有機(jī)碳的含量。

排土場(chǎng)是由人工堆積而成，地形復(fù)雜，且邊坡和平臺(tái)的過(guò)渡較自然地形更為明顯。基于以往研究，本文提出在礦區(qū)排土場(chǎng)中結(jié)合地形將其分為邊坡與平臺(tái)2 類地形，采用分區(qū)反距離權(quán)重法對(duì)排土場(chǎng)土壤養(yǎng)分進(jìn)行插值，同時(shí)與傳統(tǒng)克里金插值法和反距離權(quán)重法進(jìn)行比較。通過(guò)對(duì)3 種插值方式預(yù)測(cè)精度的比較，選出最適宜的估算模型，作為排土場(chǎng)土壤養(yǎng)分含量預(yù)測(cè)的方法，為排土場(chǎng)土壤評(píng)價(jià)奠定基礎(chǔ)。

1 材料與方法

1.1 研究區(qū)概況

研究區(qū)位于內(nèi)蒙古自治區(qū)錫林浩特市北郊，地理坐標(biāo)為東經(jīng)115°30′～116°26′，北緯 43°57′～44°14′。研究區(qū)屬半干旱草原氣候，年溫差較大，1 月平均氣溫-18.8℃，7 月平均氣溫21.2℃，年平均氣溫2.6℃。受季風(fēng)影響，降雨多集中在夏季。年平均降水量294.74 mm，年平均蒸發(fā)量1 794.64 mm。由于研究區(qū)降雨少蒸發(fā)強(qiáng)，導(dǎo)致植被多為耐旱性植物且植被覆蓋率低。隨著礦區(qū)不斷開(kāi)發(fā)，煤炭大量開(kāi)采，礦區(qū)植被受到人為破壞，生態(tài)環(huán)境不斷惡化，因此，需進(jìn)行土地復(fù)墾與植被重建。勝利煤田一號(hào)露天礦北排土場(chǎng)于2006 年開(kāi)始治理，治理時(shí)間為3 a，治理面積107 萬(wàn)m2，排土場(chǎng)高度為60 m，共分為4 個(gè)平臺(tái)(平臺(tái)1—平臺(tái)4)，3 個(gè)邊坡(邊坡1—邊坡3)，綠化面積為101 萬(wàn)m2，其中坡面35.5 萬(wàn)m2、平臺(tái)面積為65.5 萬(wàn)m2。研究區(qū)位置及地形如圖1 所示。

1.2 樣品采集與測(cè)定

圖1 研究區(qū)位置及地形分區(qū)Fig.1 Study area location and terrain division

2017 年6 月，對(duì)勝利一號(hào)露天礦北排土場(chǎng)進(jìn)行實(shí)地調(diào)查及土壤采樣，采樣點(diǎn)分布如圖1 所示。利用網(wǎng)格法布設(shè)采樣點(diǎn)，并考慮北排土場(chǎng)的地形因素，增加邊坡采樣點(diǎn)，共獲得117 個(gè)土壤樣品，舍棄17個(gè)不合格樣品，最終采用100 個(gè)樣品。所有樣品均來(lái)自0～10 cm 的表土層，利用環(huán)刀采集土樣后放入鋁盒，裝入自封袋，采樣同時(shí)記錄樣點(diǎn)編號(hào)、取樣坐標(biāo)等信息。帶回實(shí)驗(yàn)室后，測(cè)定土壤有機(jī)質(zhì)、速效氮、有效磷、速效鉀含量等。

土壤有機(jī)質(zhì)(SOM)含量通過(guò)在加熱條件下，用一定量的重鉻酸鉀-硫酸溶液氧化土壤中的有機(jī)碳，并使用油浴鍋分析獲得；速效氮含量使用半微量凱氏定氮法獲得；有效磷含量使用Olsen 法測(cè)定；速效鉀含量通過(guò)乙酸銨浸提-火焰光度法測(cè)定[8-10]。

1.3 插值方法

1.3.1 傳統(tǒng)克里金插值法

克里金空間插值法是基于變異函數(shù)模型在有限區(qū)域內(nèi)對(duì)區(qū)域化變量的取值進(jìn)行無(wú)偏最優(yōu)估計(jì)的空間插值方法[11-13]。半方差函數(shù)是地統(tǒng)計(jì)學(xué)中用以研究土壤空間變異性的關(guān)鍵函數(shù)之一，能夠定量反映土壤性質(zhì)在不同距離上觀察值之間的變化。在一維條件下，區(qū)域化變量Z(x)在x軸方向上的變異函數(shù)，被定義為區(qū)域變量Z(x)在點(diǎn)x和x+h處的值Z(x)與Z(x+h)差的方差的一半，表示為y(h)，計(jì)算公式如下：

根據(jù)實(shí)測(cè)數(shù)據(jù)及半方差函數(shù)，計(jì)算獲取變量的線性加權(quán)組合，以此對(duì)預(yù)測(cè)值進(jìn)行無(wú)偏的最佳估計(jì)。本研究區(qū)為排土場(chǎng)，根據(jù)各個(gè)采樣點(diǎn)獲取的實(shí)測(cè)數(shù)據(jù)并利用傳統(tǒng)克里金插值法進(jìn)行插值，在結(jié)合半方差函數(shù)y(h)的基礎(chǔ)上，使用若干實(shí)測(cè)點(diǎn)Z(xi)來(lái)推求出未知點(diǎn)的Z0(x0)，通過(guò)采用線性加權(quán)法表達(dá)未測(cè)點(diǎn)x0的預(yù)測(cè)值，具體見(jiàn)式(2)：

式中：y(h)為半方差函數(shù)；h為樣本間距；N(h)表示間距為矢量h的所有觀測(cè)點(diǎn)個(gè)數(shù)；Z(xi)、Z(xi+h)為Z(x)在xi和(xi+h)位置上的實(shí)測(cè)值；Z0(x0)為未知點(diǎn)處x0的線性、無(wú)偏、最優(yōu)估計(jì)值；n為實(shí)測(cè)值個(gè)數(shù)；Wi為各實(shí)測(cè)點(diǎn)權(quán)重。

1.3.2 反距離權(quán)重插值法(IDW)

IDW 即使用若干個(gè)實(shí)測(cè)點(diǎn)Z(xi,yi)來(lái)推求未知點(diǎn)Z0(x0,y0)的一種方法，一般來(lái)說(shuō)，首先要計(jì)算出預(yù)測(cè)點(diǎn)周圍每個(gè)實(shí)測(cè)點(diǎn)的權(quán)重，然后根據(jù)線性加權(quán)法計(jì)算出Z0的預(yù)測(cè)值[14]，具體見(jiàn)式(3)、式(4)：

式中：di為預(yù)測(cè)點(diǎn)周圍各實(shí)測(cè)點(diǎn)到預(yù)測(cè)點(diǎn)之間的距離。

1.3.3 分區(qū)反距離權(quán)重插值法(PIDW)

結(jié)合排土場(chǎng)的地形特征，將排土場(chǎng)分為7 個(gè)子區(qū)，利用反距離權(quán)重法在每個(gè)分區(qū)分別進(jìn)行插值，獲取各個(gè)分區(qū)的空間分布圖，并整合為整個(gè)研究區(qū)的空間分布預(yù)測(cè)圖。

1.3.4 模型驗(yàn)證方案

為檢驗(yàn)插值方法預(yù)測(cè)土壤養(yǎng)分空間分布的效果，以平均誤差ME、相對(duì)誤差RE、均方根誤差RMSE、平均絕對(duì)誤差MAE 以及R2這5 個(gè)指標(biāo)來(lái)比較不同插值方法的插值精度。多種誤差指標(biāo)的綜合使用能更好地反映預(yù)測(cè)值誤差的實(shí)際情況[15-16]。

Shanghai-Kunming high-speed rail XU Sheng LI Jia-yu ZHANG Hui-gang et al.(35)

2 結(jié)果分析

2.1 土壤養(yǎng)分指標(biāo)的描述性統(tǒng)計(jì)分析

研究區(qū)土壤養(yǎng)分的描述性統(tǒng)計(jì)結(jié)果見(jiàn)表1。根據(jù)變異系數(shù)(CV)等級(jí)劃分，CV＜10%屬于弱變異性，若10%≤CV≤100%則為中等變異，當(dāng)CV＞100%是強(qiáng)變異性。表1 顯示，速效氮、有效磷、速效鉀及有機(jī)質(zhì)的變異系數(shù)均大于10%卻小于100%。由此，4 種養(yǎng)分元素都屬于中等變異程度。

描述性統(tǒng)計(jì)分析僅能表示各種土壤養(yǎng)分含量的特征，不能反映整個(gè)研究區(qū)各土壤養(yǎng)分的空間分布情況及邊坡與平臺(tái)之間的空間分布差異。因此，需繪制土壤養(yǎng)分空間分布圖。

表1 土壤養(yǎng)分基本信息Table 1 Basic information of soil nutrients

2.2 不同插值方法的結(jié)果

在進(jìn)行克里金插值時(shí)，要求插值數(shù)據(jù)符合正態(tài)分布。本文對(duì)不符合正態(tài)分布的4 種養(yǎng)分?jǐn)?shù)據(jù)進(jìn)行對(duì)數(shù)轉(zhuǎn)換后再進(jìn)行K-S 檢驗(yàn)，結(jié)果顯示，有效磷、速效氮、速效鉀、有機(jī)質(zhì)4 種養(yǎng)分均符合正態(tài)分布。

本次選取反距離權(quán)重法與傳統(tǒng)克里金插值法進(jìn)行插值，同時(shí)由于邊坡與平臺(tái)土壤養(yǎng)分含量差異較大，再選取反距離權(quán)重法進(jìn)行分區(qū)插值，共3 種空間插值方法。通過(guò)ME、MAE、RMSE、RE、R2共5 個(gè)指標(biāo)的驗(yàn)證結(jié)果，得出最適合研究區(qū)的插值方法。表2 為3 種插值方法交叉驗(yàn)證結(jié)果。

據(jù)表2 可知，4 種養(yǎng)分的空間分布插值結(jié)果中，插值效果最好的方法均為分區(qū)反距離權(quán)重法。對(duì)有效磷和有機(jī)質(zhì)而言，反距離權(quán)重法次之，傳統(tǒng)克里金插值法精度最差。對(duì)速效氮而言，傳統(tǒng)克里金插值法次之，反距離權(quán)重法精度最差。對(duì)速效鉀而言，各種指標(biāo)均為3 種方法的最小值。

根據(jù)圖2 顯示，4 種土壤養(yǎng)分元素采用分區(qū)反距離權(quán)重法插值的R2最高，其中速效氮與速效鉀的R2提升較為明顯，達(dá)到0.30 以上。相較來(lái)說(shuō)，反距離權(quán)重法對(duì)土壤養(yǎng)分的R2整體上也要優(yōu)于傳統(tǒng)克里金插值法。以土壤N 為例，分區(qū)反距離權(quán)重法得到的研究區(qū)土壤速效氮含量為8.09～63.73 mg/kg，傳統(tǒng)克里金方法得到的含量為11.22～55.62 mg/kg，而反距離權(quán)重法得到的含量為11.12～63.68 mg/kg；相較于反距離權(quán)重法和傳統(tǒng)克里金插值法，分區(qū)反距離權(quán)重法對(duì)土壤N 含量的預(yù)測(cè)范圍更加準(zhǔn)確，對(duì)土壤N 的預(yù)測(cè)范圍與實(shí)測(cè)值相差不大，能較好地預(yù)計(jì)土壤養(yǎng)分的極值點(diǎn)(較大值和較小值)，反距離權(quán)重法次之，傳統(tǒng)克里金插值結(jié)果最差，其他土壤元素(P、K 和有機(jī)質(zhì))也有相同的結(jié)果。這是因?yàn)榭死锝鸩逯祵?duì)于預(yù)測(cè)結(jié)果有平滑效果，使其對(duì)極值點(diǎn)的預(yù)測(cè)效果較差，但較小值對(duì)于排土場(chǎng)而言屬于土壤貧瘠區(qū)，不利于植被的恢復(fù)，應(yīng)當(dāng)受到重視。綜合比較不同空間插值方法的交叉驗(yàn)證結(jié)果顯示，針對(duì)本研究區(qū)，分區(qū)反距離權(quán)重法插值的精度相對(duì)較高。

表2 不同插值法交叉驗(yàn)證結(jié)果Table 2 Cross-validation results for different interpolation methods

圖2 不同插值法驗(yàn)證點(diǎn)散點(diǎn)圖Fig.2 Scatter plots of verification points for different interpolation methods

3 討論

圖3 北排土場(chǎng)土壤養(yǎng)分空間分布Fig.3 Spatial distribution of soil nutrients in the north dump

據(jù)圖3 可知，排土場(chǎng)的速效氮含量整體呈現(xiàn)出由北至南逐漸減少的規(guī)律；而有效磷含量均偏低，僅中心地區(qū)含量有所提升；整體來(lái)看排土場(chǎng)的速效鉀含量偏低，但東北地區(qū)含量處于水平值之上；土壤有機(jī)質(zhì)整體含量偏高，中心地區(qū)含量相對(duì)偏低。平臺(tái)與邊坡的養(yǎng)分分布空間差異較大，但這種空間差異并沒(méi)有表現(xiàn)出規(guī)律性[17]。這是由于排土場(chǎng)的土壤是由人為堆積而成，與自然土壤的分布規(guī)律不同，表現(xiàn)出明顯的空間異質(zhì)性。文獻(xiàn)[18]研究發(fā)現(xiàn)邊坡區(qū)域?qū)λ鞯谋３肿饔帽绕脚_(tái)差，且水土流失、土壤受侵蝕狀況與坡度成正比，使邊坡土壤養(yǎng)分低于平臺(tái)。但本研究并無(wú)此類結(jié)論，這可能是由于研究區(qū)排土場(chǎng)復(fù)墾年限較短，且降水較少，多以風(fēng)蝕為主，致使邊坡平臺(tái)之間養(yǎng)分分布并無(wú)規(guī)律。

分區(qū)反距離權(quán)重插值方法基于排土場(chǎng)的地形在常規(guī)反距離插值方法上進(jìn)行了改進(jìn)，但是在每個(gè)分區(qū)上的反距離插值方法仍需遵循“地理學(xué)第一定律”[19]，否則插值效果不會(huì)得到改進(jìn)，因此，分區(qū)反距離權(quán)重法仍然具有一定的局限性。交叉驗(yàn)證結(jié)果顯示，速效氮和速效鉀的分區(qū)反距離權(quán)重插值結(jié)果明顯好于反距離權(quán)重法和傳統(tǒng)克里金插值法，但有效磷的插值結(jié)果顯示3 種插值方法沒(méi)有明顯差異，此結(jié)果表明，有效磷在分區(qū)后，插值點(diǎn)仍然是空間不平穩(wěn)的，因此，不適應(yīng)本文采取的所有插值方法。

本次基于地形分區(qū)，使插值點(diǎn)在相同地形條件下進(jìn)行插值，屬于局部插值方法，這一插值方法同時(shí)也被應(yīng)用于類似的研究中?；诰植坎逯邓枷?，徐占軍等[11]根據(jù)積水狀態(tài)對(duì)采煤沉陷區(qū)土壤有機(jī)碳含量進(jìn)行分區(qū)克里金插值并與傳統(tǒng)克里金插值進(jìn)行比較，得出分區(qū)克里金模擬更能反映土壤有機(jī)碳的空間遞變特點(diǎn)，有利于分析不同因素對(duì)土壤有機(jī)碳空間分布的影響。吳子豪等[20]利用3 種結(jié)合土地利用的克里金插值模型對(duì)農(nóng)田土壤有機(jī)碳含量進(jìn)行預(yù)測(cè)，發(fā)現(xiàn)中位數(shù)中心化克里金法結(jié)合土地利用類型可生成精度更高的土壤有機(jī)碳密度空間分布圖。

4 結(jié)論

a.排土場(chǎng)的速效氮含量整體呈現(xiàn)出由北至南逐漸減少的規(guī)律；有效磷含量均偏低，僅中心地區(qū)含量有所提升；排土場(chǎng)的速效鉀含量偏低，東北部含量偏高；土壤有機(jī)質(zhì)整體含量偏高，中心地區(qū)含量相對(duì)偏低。

b.根據(jù)分區(qū)反距離權(quán)重法得到排土場(chǎng)土壤速效氮含量為 8.09～63.73 mg/kg，速效鉀含量為37.19～488.26 mg/kg，有效磷含量為0.89～6.76 mg/kg，有機(jī)質(zhì)含量為2.08～21.3 g/kg，與實(shí)測(cè)數(shù)據(jù)更接近，分區(qū)反距離權(quán)重法對(duì)養(yǎng)分含量的預(yù)測(cè)范圍較為精確，對(duì)極值預(yù)測(cè)精度較高。

c.平臺(tái)與邊坡的土壤養(yǎng)分空間差異較大，表現(xiàn)出明顯的空間異質(zhì)性，但并無(wú)明顯規(guī)律。由于排土場(chǎng)土壤由人為堆積而成，與自然土壤分布規(guī)律不同，且研究區(qū)排土場(chǎng)復(fù)墾年限較短，降雨較少，土壤侵蝕多以風(fēng)蝕為主，致使邊坡平臺(tái)間養(yǎng)分分布并無(wú)規(guī)律。

d.從預(yù)測(cè)精度上來(lái)看，北排土場(chǎng)最優(yōu)插值方法為分區(qū)反距離權(quán)重法。這是由于研究區(qū)地形復(fù)雜，坡度起伏較大且土壤分布受人為影響導(dǎo)致土壤中養(yǎng)分含量分布復(fù)雜，結(jié)合基于地形的分區(qū)反距離權(quán)重法插值，可以有效消除不同地形采樣點(diǎn)土壤養(yǎng)分含量對(duì)空間插值預(yù)測(cè)精度的影響，從而提高土壤養(yǎng)分含量的預(yù)測(cè)精度。

請(qǐng)聽(tīng)作者語(yǔ)音介紹創(chuàng)新技術(shù)成果等信息，歡迎與作者進(jìn)行交流

OSID 碼