朱洪芬,南 鋒,徐占軍,荊耀棟,段永紅,畢如田

山西農(nóng)業(yè)大學(xué)資源環(huán)境學(xué)院,太谷 030801

土壤有機質(zhì)含量不僅是肥力的重要指標[1],也是全球碳循環(huán)中重要的源和匯[2- 3]。前人已對不同尺度下土壤有機質(zhì)的空間變異性進行了大量研究,證明了在特定尺度下,土壤的異質(zhì)性與外部環(huán)境因素共同影響有機質(zhì)的空間分布[3- 5]。因此,建立土壤有機質(zhì)與影響因子的準確關(guān)系成為間接理解有機質(zhì)空間變異性的方法之一。然而,由于不同尺度下的土壤有機質(zhì)含量與影響因子的空間關(guān)系具有較大差異性[5- 6],當某一因子在不同尺度同時影響土壤有機質(zhì)空間分布時,可能造成該因子與土壤有機質(zhì)的非線性關(guān)系。因此,應(yīng)用傳統(tǒng)統(tǒng)計方法建立土壤有機質(zhì)預(yù)測模型時,模型的精確度與可靠性可能會受到較大影響。

當前,在環(huán)境因子與土壤屬性的空間多尺度關(guān)系的研究中,多采用Pearson相關(guān)分析[7]、小波分析[8- 9]和多維分析[10]等方法。應(yīng)用該類方法的前提是假設(shè)相關(guān)環(huán)境因子對土壤屬性空間分布的影響為線性且平穩(wěn),但此類假設(shè)可能與土壤屬性的實際空間分布不相符。多元經(jīng)驗?zāi)B(tài)分解(multivariate empirical mode decomposition, MEMD)由Huang等首次提出,是與小波分析相對應(yīng)的另一種多尺度分析方法,可根據(jù)數(shù)據(jù)本身的尺度特征經(jīng)驗性的將空間數(shù)據(jù)分解在多個表征尺度上[11],避免了線性和平穩(wěn)性的假設(shè)。盡管MEMD法具有明顯優(yōu)勢,但尚未廣泛應(yīng)用于土壤屬性尺度問題的相關(guān)研究中。

黃土高原由第四紀冰期厚層黃土堆積而成,是全球水土流失最嚴重的地區(qū)之一[12]。黃土高原盆地內(nèi)多平原與丘陵,海拔高度相對較低,溫濕環(huán)境適宜,有利于農(nóng)業(yè)生產(chǎn)[13]。由于土地利用形式存在多樣性及植被覆蓋的斑塊化與破碎化,造成了本區(qū)域內(nèi)土壤有機質(zhì)分布的無序性及影響因子對有機質(zhì)作用的非線性[14]。因此,研究該區(qū)域土壤有機質(zhì)與相關(guān)因子間的空間多尺度關(guān)系可為有機質(zhì)管理提供理論依據(jù)。鑒于MEMD法可用于空間數(shù)據(jù)的非線性和非平穩(wěn)的空間多尺度分析,因此本研究假設(shè)該方法也可用于土壤有機質(zhì)與相關(guān)因子間的空間多尺度研究中。本文以山西省太原盆地為研究區(qū)域,分析盆地內(nèi)不同部位土壤有機質(zhì)的表征尺度,以及土壤有機質(zhì)與影響因子的空間多尺度關(guān)系,并實現(xiàn)有機質(zhì)含量在采樣尺度上的預(yù)測。

1 研究材料與方法

1.1 研究區(qū)概況

太原盆地地處黃土高原中東部、山西省中部(36°55′—38°12′N,111°42′—113°02′E),屬典型的黃土集中分布地帶。盆地南北長約150 km,東西寬為12—60 km,區(qū)域包括汾河流域中游,總面積達6159 km2,人口數(shù)量占全省近1/3。該區(qū)氣候?qū)倥瘻貛Т箨懶约撅L(fēng)氣候類型,年日照總時數(shù)為2360—2796 h,年平均降水量為420—457 mm,總降水量由南到北逐漸減少。境內(nèi)平均海拔800—900 m,東、西、北三面環(huán)山,中部為沖積平原,年平均氣溫7.8—10.3℃。中部平原以沖洪積亞砂粘土質(zhì)黃土為主,邊山丘陵以風(fēng)積厚層亞砂黃土為主。境內(nèi)土壤類型主要有潮土、鹽化潮土、石灰性褐土和褐土性土等四類,主要農(nóng)作物為冬小麥和玉米,是山西省主要農(nóng)業(yè)生產(chǎn)區(qū)。

1.2 樣品的采集與處理

根據(jù)海拔與汾河流向,沿北東-南西方向?qū)⑴璧貏澐譃樯稀⒅?、下三部位。結(jié)合野外調(diào)查和相關(guān)圖件分析,避開大型建設(shè)用地后,在盆地不同位置設(shè)置3條采樣帶(帶1海拔：742—894 m；帶2海拔：723—807 m；帶3海拔：707—1002 m),每條樣帶長約42 km,以330 m間隔設(shè)置采樣點。若某一樣點落于建設(shè)用地或道路等非耕地上,則在離該點最近的耕地上取樣,并盡量使所有樣點在一條直線上。采樣帶1、2和3分別包含121、128和134個樣點,共計383個(圖1)。

圖1 研究區(qū)地理位置及采樣點分布Fig.1 Geographical location of study area and sample points distribution

土樣采于2016年3月22—31日。采樣時,首先利用GPS查找樣點位置,并采用“S”形布點法進行樣品采集。使用螺旋取土鉆鉆取0—20 cm耕層土壤樣品5個,混合后作為本采樣點樣品。樣品混合均勻后經(jīng)風(fēng)干、磨碎、過2 mm孔篩后分為兩份,分別用于土壤光譜和土壤有機質(zhì)、質(zhì)地測試。使用環(huán)刀(高5 cm,體積100 cm3)在“S”形樣點中心處采集表層原狀土樣,烘干后(105℃,10 h)測定土壤容重。采用重鉻酸鉀-外加熱法測定土壤有機質(zhì)含量[15],采用沉降法測定土壤質(zhì)地[16],采用地物光譜儀(美國ASD FieldSpec3)測定土壤光譜。其中,光譜范圍為350—2500 nm,數(shù)據(jù)重采樣間隔為1 nm[17]。

1.3 MEMD原理

MEMD是經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition, EMD)的擴展,可直接將獲取的數(shù)據(jù)應(yīng)用于空間域,并根據(jù)空間數(shù)據(jù)特征,將其分解為多個空間序列。MEMD對數(shù)據(jù)的分解基于如下假設(shè),即在給定的空間內(nèi),空間數(shù)據(jù)存在簡單的不同頻率相互疊加的震蕩模式[11]。自然界中,事物整體變異性受多種過程影響,并在不同尺度以不同強度發(fā)生[18]。其中,同一尺度發(fā)生的過程可以分解為相同的本征模函數(shù)(Intrinsic Mode Function, IMF)。定義IMF需滿足以下條件：1)可以是線性或非線性,該IMF在整個數(shù)據(jù)范圍內(nèi),局部極值點和過零點的個數(shù)必須相等或最多相差一個。其中,局部極值點表示一個函數(shù)可以提取的最小值或最大值,過零點表示函數(shù)改變其符號的點。2)該數(shù)據(jù)的震蕩對局部平均值對稱,即在任一數(shù)據(jù)點處由局部極大值和極小值定義的平均包絡(luò)值為0。上包絡(luò)線和下包絡(luò)線由通過分別樣條插值全部局部最大值或最小值產(chǎn)生。根據(jù)該定義,IMF可通過“篩選(sifting)”過程分解任意函數(shù)后獲取,IMF的方差貢獻百分比(%)可通過下列公式計算：

(1)

單個IMF對整體方差的貢獻率決定了每個IMF的相對重要性,因此,也可決定在不同尺度上發(fā)生過程的相對重要性[19],即方差分解法。通過計算每個IMF震蕩的平均次數(shù),可確定對應(yīng)IMF的平均尺度。然而,由于土壤屬性與環(huán)境因子之間存在非線性關(guān)系,使得局部范圍內(nèi)無法預(yù)測其周期或尺度。瞬時尺度可提供每個IMF所代表的尺度范圍,并通過Hilbert變換得到的瞬時頻率來獲取。EMD可將空間數(shù)據(jù)分解為不同的表征尺度,而Hilbert譜分析可計算它們在每個尺度和位置的能量與頻率,便于從不同尺度的數(shù)據(jù)中提取不同位置的頻率。

應(yīng)用MEMD時要求計算空間數(shù)據(jù)序列的局部平均值,而對于多個空間數(shù)據(jù)序列可能無法直接定義最大值和最小值。為解決該問題,Rehman和Mandic[20]建議,在N維空間內(nèi)通過不同方向的投影創(chuàng)建多個N維包絡(luò),取其均值可作為局部平均值。

(1)產(chǎn)生一個合適的方向矢量X數(shù)據(jù)集；

(2)計算空間序列V(s)沿給定方向Xθk的投影pθk(s)；

(5)通過下式計算包絡(luò)曲線的均值M(s)；

(2)

(6)利用D(s)=V(s)-M(s)提取D(s)。若D(s)滿足終止迭代準則,則從步驟(1)開始重復(fù)以上步驟,并計算殘差V(s)-D(s)。否則,從步驟(2)開始重復(fù)計算D(s)。殘差變?yōu)閱握{(diào)函數(shù)且再無法提取IMF時,分解過程終止。該殘差可表明原始數(shù)據(jù)的變化趨勢。

1.4 數(shù)據(jù)預(yù)處理

去除土壤光譜數(shù)據(jù)中噪聲影響較大的350—399和2451—2500 nm邊緣波段,利用MATLAB將可見光-近紅外波段(401—2450 nm)2050個土壤光譜數(shù)據(jù)進行主成分壓縮,并將壓縮后前兩個主份(占整體方差的95%以上)選作綜合影響因子。

將土壤有機質(zhì)與地形因子(高程、坡度和地形濕度指數(shù))、物理性狀(容重、砂粒、壤粒和黏粒含量)、光譜主份數(shù)據(jù)等組成多元數(shù)據(jù)序列,利用MEMD法將各樣帶數(shù)據(jù)序列分解為不同IMF。根據(jù)MEMD算法,通過多元序列數(shù)據(jù)中的N個IMF共同震蕩模式來識別每個“共同尺度”[20]。在不同的數(shù)據(jù)源中,MEMD將相似的尺度歸為一組,代表相關(guān)發(fā)生過程的實際尺度。

計算各樣帶在采樣尺度與單個IMF(表征尺度)和殘差下的土壤有機質(zhì)與地形因子、土壤物理性狀、土壤光譜主份等因子的Pearson相關(guān)系數(shù)。采用逐步多元回歸模型,利用地形、土壤物理性狀和光譜主份等被分解的各IMF和殘差,預(yù)測土壤有機質(zhì)在對應(yīng)尺度和殘差處的含量值。最后,通過單個IMF和殘差預(yù)測出的對應(yīng)土壤有機質(zhì)值,利用如下公式估測采樣尺度上有機質(zhì)含量：

(3)

采用由實測值和預(yù)測值計算出的決定系數(shù)(coefficient of determination,R2)、均方根誤差(root mean squared error, RMSE)、標準化均方根誤差(normalized root mean square error, NRMSE)和相對預(yù)測偏差(residual prediction deviation,RPD)評價土壤有機質(zhì)預(yù)測精度,各參數(shù)計算公式如下：

(4)

(5)

(6)

(7)

2 結(jié)果與分析

2.1 采樣尺度上土壤有機質(zhì)與影響因子的關(guān)系

表1為利用Pearson′s相關(guān)分析法對土壤有機質(zhì)與地形因子、土壤物理性狀、光譜主份等影響因子的線性相關(guān)度進行分析的結(jié)果。在樣帶1上,有機質(zhì)含量與高程、容重和光譜主份1顯著相關(guān)；在樣帶2上,與高程、坡度、容重、砂粒、壤粒、黏粒和光譜主份1顯著相關(guān)；在樣帶3上,與高程、容重、砂粒、壤粒和光譜主份1顯著相關(guān)。在全部樣帶上,有機質(zhì)與高程、容重、砂粒、壤粒和光譜主份1顯著相關(guān),表明在采樣尺度上研究區(qū)內(nèi)土壤有機質(zhì)主要受這五個影響因子的作用,相關(guān)局部尺度的發(fā)生過程(如生物活動、耕作措施等)擾亂了其他影響因子對土壤有機質(zhì)的作用,因此表現(xiàn)為無顯著相關(guān)性。土壤有機質(zhì)與物理性狀的相關(guān)性順序為樣帶2>樣帶3>樣帶1。樣帶2中,地形因子、容重、質(zhì)地對土壤有機質(zhì)含量的影響最強。這些結(jié)果表明,在采樣尺度上,地形因子與土壤物理性狀對盆地中部土壤有機質(zhì)含量的影響最強。

表1 采樣尺度上土壤有機質(zhì)與影響因子的相關(guān)系數(shù)

*顯著性水平為P<0.05；**顯著性水平為P<0.01

如表2所示,本研究采用逐步多元回歸法構(gòu)建了土壤有機質(zhì)預(yù)測模型。各相關(guān)因子(樣帶1中為壤粒、光譜主份1和2,樣帶2中為坡度、容重、壤粒和光譜主份1,樣帶3中為高程、容重、砂粒、黏粒、光譜主份1和2)對樣帶1、2和3中土壤有機質(zhì)含量的預(yù)測能力分別達到了52%、56%和54%。

2.2 土壤有機質(zhì)與影響因子的空間多尺度關(guān)系

利用MEMD法將土壤有機質(zhì)與影響因子的多元數(shù)據(jù)序列分解為不同的IMF。土壤有機質(zhì)的IMF如圖2所示。隨著IMF增大,土壤有機質(zhì)空間序列的震蕩周期變長,表明土壤有機質(zhì)的表征尺度隨IMF變大而增加。同時,各樣帶在IMF1處的震蕩幅度均較大,表明在該尺度處的空間變異性較大。土壤有機質(zhì)的殘差序列表明了原始數(shù)據(jù)的變換趨勢,即樣帶1中土壤有機質(zhì)的變化趨勢比較明顯,且與樣帶2一起呈逐漸減小的趨勢,而樣帶3呈先穩(wěn)定后增大的趨勢。各樣帶上,土壤有機質(zhì)與影響因子均具有相同的IMF個數(shù)與類似的震蕩幅度。因此,本研究利用Hilbert變換識別各樣帶不同因子的IMF對應(yīng)的空間尺度,并將所有因子的尺度取其均值,獲取該樣帶上每個IMF的對應(yīng)尺度(表3)。樣帶中各因子的IMF對應(yīng)尺度的變異系數(shù)隨IMF變大而增大,表明在小尺度處土壤有機質(zhì)和影響因子具有相近的表征尺度。樣帶2和3中,1—4 IMF尺度較接近,即在< 5000 m尺度的發(fā)生過程約在表征尺度960、1500、2600和4400 m處；5—8 IMF尺度差異較大,即樣帶2發(fā)生過程的表征尺度約為8573、10901、11591和23659 m,而樣帶3約為6753、11806和12292 m。該結(jié)果表明,盆地中、下部在小尺度處的表征尺度較接近；隨著尺度的變大,表征尺度的差異性增大。另外,盆地上部的表征尺度與盆地中部、下部的差異較大。

表2 基于影響因子的土壤有機質(zhì)的逐步多元回歸預(yù)測模型

Elevation：高程；Slope：坡度；TWI：地形濕度指數(shù),topographic wetness index；BD：容重,bulk density；Sand：砂粒；Silt：壤粒；Clay：黏粒；PC1：光譜主份1；PC2：光譜主份2；括號中的數(shù)據(jù)表示標準偏回歸系數(shù)

圖2 三條樣帶中土壤有機質(zhì)被分解的IMF和殘差Fig.2 Separated IMFs and residues for SOM at the three transects

樣帶中有機質(zhì)與影響因子的IMF方差百分比如表4所示。樣帶1中,土壤有機質(zhì)方差百分比主要分布在IMF1和2處(18.62%和19.41%),表明盆地上部的土壤有機質(zhì)變異性主要在1011和1725 m尺度處。樣帶2中,其方差百分比主要分布在IMF1和5處(30.61%和16.91%),表明盆地中部的土壤有機質(zhì)變異性主要在982和8573 m尺度處。樣帶3中,其方差百分比主要分布在IMF1、5和6處(17.62%、12.53%和20.46%),表明該處土壤有機質(zhì)變異性主要表現(xiàn)在960、6753和11806 m尺度處。總之,整個區(qū)域內(nèi)尺度約1000 m均表現(xiàn)為土壤有機質(zhì)的主要表征尺度。同時,盆地上部土壤有機質(zhì)的主要表征尺度是小尺度,即< 2000 m；盆地中部的是小尺度和中尺度,即1000、8000 m；盆地下部的是小尺度、中尺度和大尺度,即1000、7000、12000 m。這些結(jié)果表明,該區(qū)域內(nèi)垂直于汾河流向的土壤有機質(zhì)主要表征尺度沿河流方向表現(xiàn)分散。

表3 三個樣帶中土壤有機質(zhì)IMF的對應(yīng)尺度/m

括號中的數(shù)據(jù)表示土壤有機質(zhì)和影響因子IMF對應(yīng)尺度的變異系數(shù)值(%)

表4 三條樣帶中土壤有機質(zhì)的每個IMF和殘差占原始數(shù)據(jù)方差的百分比

通過觀察各影響因子的IMF方差百分比發(fā)現(xiàn),高程除在第三條樣帶IMF6處值(21.07%)較大外,其余主要分布在殘差部分(3條樣帶分別為75.93%、72.83%和70.04%),表明高程的主要表征尺度比該實驗提取的尺度大。在樣帶1中,除高程外的其他影響因子的表征尺度主要在IMF1和2處,即小尺度1011和1725 m處。樣帶2中,坡度、砂粒和壤粒含量的表征尺度主要在IMF1和5處(尺度982 和8573 m),其余影響因子主要在IMF1和2處(尺度982 和1530 m)。樣帶3中,除坡度(IMF1、3和6處)、黏粒(IMF1和3處)、光譜主份1和2(IMF1、5和6處)外,其余影響因子主要在IMF1和2處(960和1504 m)。上述結(jié)果表明,除高程外相關(guān)因子的變異性主要發(fā)生在小尺度處。

與采樣尺度上的相關(guān)性不同,土壤有機質(zhì)與影響因子多尺度相關(guān)性見表5。除樣帶3中IMF2外,高程與土壤有機質(zhì)的相關(guān)性主要表現(xiàn)在大尺度,即樣帶1為IMF5、6和殘差；樣帶2為IMF6、7、8和殘差；樣帶3為IMF6、7和殘差。這些結(jié)果表明,盆地內(nèi)高程對土壤有機質(zhì)的影響主要表現(xiàn)在大尺度約10000、12000、23000 m處。坡度與有機質(zhì)關(guān)系在樣帶1中較弱,主要表現(xiàn)在大尺度IMF6和殘差處,在樣帶3中主要表現(xiàn)在IMF2、5、6、7和殘差處；而與樣帶2中有機質(zhì)關(guān)系最為顯著,即除IMF1和3外,其余均顯著相關(guān)。地形濕度指數(shù)與土壤有機質(zhì)的空間多尺度關(guān)系表明,在盆地中部兩者的關(guān)系最明顯,主要表現(xiàn)在尺度982、2588、4487、8573、10901、11591 m處。盆地下部兩者的關(guān)系較明顯,這與采樣尺度上兩者的關(guān)系不一致。土壤容重與有機質(zhì)的多尺度關(guān)系表明,在盆地上、中部小尺度約1000 和1500 m處,中、下部尺度約3000 m處,下部尺度約4500 和6700 m處,上部尺度約5400和9500 m處及在整個盆地內(nèi)大尺度> 10000 m處,兩者顯著相關(guān)。土壤質(zhì)地(包括砂粒、壤粒和黏粒)與有機質(zhì)的多尺度關(guān)系表明,壤粒含量對土壤有機質(zhì)的多尺度關(guān)系最顯著,即除樣帶1中IMF4外,二者在所有表征尺度均顯著相關(guān)。光譜主份與有機質(zhì)的多尺度關(guān)系表明,在研究樣帶中的所有表征尺度上,光譜主份1均與有機質(zhì)含量顯著相關(guān),而光譜主份2的相關(guān)性明顯弱于光譜主份1?？傊?相關(guān)因子對土壤有機質(zhì)的影響隨尺度增大而呈增大趨勢。另外,對于殘差部分,土壤有機質(zhì)與影響因子的相關(guān)性均達到0.05的顯著性水平,表明被選定的影響因子與土壤有機質(zhì)存在一定關(guān)系。

表5 基于MEMD的不同IMF和殘差的土壤有機質(zhì)與影響因子的相關(guān)系數(shù)

*顯著性水平為P<0.05；**顯著性水平為P<0.01

2.3 基于與影響因子多尺度關(guān)系的土壤有機質(zhì)預(yù)測

如表6所示,本研究采用逐步多元回歸獲取了每個IMF土壤有機質(zhì)的預(yù)測模型。從中可以看出,除樣帶3中IMF2外,所有預(yù)測模型的可調(diào)整R2均在0.44—1.00之間變化,且隨IMF增大而增大。可調(diào)整R2的變化趨勢表明,尺度越大,對土壤有機質(zhì)的預(yù)測精度越高,且模型中所選因子對有機質(zhì)的影響更強烈?；谒蠭MF和殘差的有機質(zhì)預(yù)測結(jié)果,本研究采用逐步多元回歸法進行采樣尺度上土壤有機質(zhì)的預(yù)測,結(jié)果列于表7。結(jié)果表明：采樣尺度上各樣帶中土壤有機質(zhì)預(yù)測值和實測值的R2分別為0.90、0.86和0.91,顯著高于采樣尺度上直接利用逐步多元回歸擬合的結(jié)果(0.52、0.56和0.54)。通過MEMD方法獲取的RMSE和NRMSE顯著低于直接逐步多元回歸的預(yù)測結(jié)果,而RPD顯著高于直接逐步多元回歸的預(yù)測結(jié)果。基于上述評價參數(shù)的比較,可以得出采用MEMD對土壤有機質(zhì)的預(yù)測精度高于基于原始數(shù)據(jù)的簡單逐步多元回歸預(yù)測的結(jié)論。

表6基于經(jīng)驗?zāi)B(tài)分解的每個IMF和殘差的土壤有機質(zhì)逐步多元回歸預(yù)測模型和回歸統(tǒng)計特征(F值和可調(diào)整R2)

Table6Predictiveequationsandregressionstatistics(F-value and adjustedR2)forSOMforeachIMFandresidueusingstepwisemultiplelinearregressionbasedonmultivariateempiricalmodedecomposition

樣帶TransectsIMF公式FunctionFR2樣帶1IMF10．03+21．31(0．43)Silt+3．77(0．68)PC1-4．94(0．12)PC242．00．52Transect1IMF20．01+1．06(0．17)TWI-12．01(0．19)BD-27．28(0．46)Sand-9．66(0．26)Clay+3．63(0．53)PC1-21．01(0．42)PC230．00．61IMF30．05+0．07(0．16)Elevation+0．26(0．13)Slope+1．15(0．24)TWI+9．74(0．44)BD+22．34(0．37)Silt-18．42(0．31)Clay+3．63(0．81)PC153．20．77IMF4-0．06+0．03(0．09)Elevation-0．34(0．17)Slope+4．96(0．07)Silt+0．69(0．90)PC1142．60．83IMF5-0．03+0．23(0．48)Elevation+1．46(0．70)Slope+1．55(0．27)TWI-7．08(0．48)BD+58．61(0．59)Silt+6．17(0．10)Clay+4．54(0．72)PC1-0．21(0．24)PC21808．60．99IMF60．01+0．46(0．98)Elevation+3．62(0．75)Slope+6．44(0．18)TWI-3．75(0．25)BD-157．71(2．27)Sand-204．82(1．49)Clay+7．41(1．38)PC1+7．49(0．07)PC261428．31．00殘差-1007．72+11．61(2．02)Slope+138．15(1．38)BD+19．51(2．44)PC1+251．95(4．02)PC2394573．81．00樣帶2Transect2IMF10．07+0．49(0．14)Slope-1．41(0．27)TWI-10．56(0．26)BD+6．77(0．12)Silt+3．32(0．51)PC121．40．47IMF20．03-0．16(0．11)Elevation-1．03(0．41)Slope-7．91(0．20)BD-3．33(0．37)Sand-29．89(0．59)Clay+1．95(0．4)PC1+4．66(0．15)PC224．80．59IMF30．08-0．58(0．30)Elevation-1．12(0．29)TWI-5．67(0．13)BD-15．07(0．26)Clay+4．62(0．92)PC1-17．81(0．44)PC248．80．71IMF4-0．07+0．71(0．28)Slope+1．71(0．36)TWI-2．56(0．06)BD+1．42(0．16)Silt+4．30(0．79)PC1-13．93(0．34)PC2165．50．89IMF5-0．01-0．55(0．51)Elevation-0．52(0．18)Slope+2．67(0．33)TWI+30．02(0．72)Silt-28．44(0．51)PC2418．70．95IMF60．33(0．85)Elevation+0．57(0．52)Slope-2．42(1．22)TWI-37．26(0．68)BD+18．02(0．71)Clay+9．25(2．32)PC1-16．24(0．41)PC210392．01．00IMF70．16(0．30)Elevation-2．73(0．72)TWI-11．86(0．34)Sand-16．14(0．17)Silt+11．10(1．39)PC1-33．83(0．25)PC212571476．71．00IMF8-0．03(0．18)Elevation+1．50(0．10)TWI-46．05(1．10)BD+11．81(0．36)PC21525552．41．00殘差-32．33+1．41(0．42)Slope+32．04(0．69)Clay+17．29(0．17)PC2816066．51．00樣帶3IMF1-0．06-23．23(0．43)Sand-19．10(0．27)Clay+3．81(0．61)PC1-11．97(0．27)PC225．40．44Transect3IMF2-0．52(0．16)TWI+13．43(0．27)Silt+3．26(0．50)PC1-14．41(0．31)PC219．50．38IMF30．01-0．11(0．27)Elevation+0．44(0．18)Slope+1．03(0．33)TWI-9．36(0．21)BD+11．81(0．28)Silt+2．54(0．50)PC1-6．04(0．18)PC217．20．49IMF40．20+2．98(0．66)TWI+12．61(0．26)BD+4．53(0．31)Silt-17．50(0．19)Clay+5．77(1．01)PC1216．30．89IMF50．05-0．16(0．45)Elevation+2．91(0．44)Slope+1．01(0．12)TWI-15．34(0．19)BD-8．10(0．07)Silt-47．86(0．35)Clay+2．52(0．53)PC1-32．90(0．64)PC2323．90．95IMF6-0．4-0．24(2．08)Elevation+6．90(1．14)Slope-13．72(1．36)TWI-18．46(0．08)BD-48．27(0．22)Silt+243．41(0．72)Clay2766．60．99IMF77．91E-5-0．03(0．15)Elevation-2．31(0．21)Slope-1．16(0．03)TWI-78．23(0．46)BD+3．64(0．49)PC1-16．36(0．23)PC21487634911．00殘差117．97-0．73(0．26)Slope+46．93(0．84)Silt-43．50(1．72)PC214528671．00

Elevation：高程；Slope：坡度；TWI：地形濕度指數(shù),topographic wetness index；BD：容重,bulk density；Sand：砂粒；Silt：壤粒；Clay：黏粒；PC1：光譜主份1；PC2：光譜主份2；括號中的數(shù)據(jù)表示標準偏回歸系數(shù)

表7 評價土壤有機質(zhì)預(yù)測精度的相關(guān)參數(shù)

RMSE：均方根誤差,root mean squared error；NRMSE：標準化均方根誤差,normalized root mean square error；RPD：相對預(yù)測偏差,residual prediction deviation

采樣尺度上土壤有機質(zhì)實際值與單個IMF或殘差處預(yù)測值及變量預(yù)測有機質(zhì)值的相關(guān)系數(shù)如圖3所示。該相關(guān)系數(shù)的比較可反映每一IMF相對于采樣尺度上有機質(zhì)預(yù)測結(jié)果的重要性程度。結(jié)果表明,樣帶1中有機質(zhì)預(yù)測的主要貢獻者是IMF6；樣帶2中主要貢獻者是IMF5；樣帶3中主要貢獻者是IMF6。即盆地上部尺度12375 m處、中部尺度8573 m處和下部尺度11806 m處對采樣尺度上有機質(zhì)的預(yù)測起主要作用。另外,可能代表更大尺度的殘差在樣帶1中的比重較大,其次是樣帶2,最后是樣帶3,表明被分解的IMF能夠很好解釋土壤有機質(zhì)變異性的順序為盆地的下部>中部>上部。

圖3 采樣尺度上土壤有機質(zhì)實際值與單個IMF或殘差處預(yù)測值及變量預(yù)測有機質(zhì)值的相關(guān)系數(shù)Fig.3 Coefficient of determination between SOM at the measurement scale and predicted IMFs (residue) or total SOM predicted by each variable

土壤有機質(zhì)的變異性主要是土壤本身異質(zhì)性、地形、植被、人類活動等綜合因素影響的結(jié)果。不同尺度的某一因子對有機質(zhì)預(yù)測值與采樣尺度上實測值的相關(guān)系數(shù)如圖3所示。該結(jié)果可表示每一因子對采樣尺度上土壤有機質(zhì)預(yù)測的相對重要性。3條樣帶中光譜主份1(綜合因子)是有機質(zhì)預(yù)測中最重要的影響因子。此外,樣帶1中容重、樣帶2中砂粒和壤粒含量、樣帶3中地形濕度指數(shù)、坡度和容重對其有機質(zhì)預(yù)測也起了重要作用?？傊?地形因子、土壤容重和光譜主份對樣帶3中有機質(zhì)的影響最為明顯,土壤質(zhì)地對樣帶2中有機質(zhì)的影響最為明顯,且影響順序為中部>下部>上部,而樣帶1中土壤有機質(zhì)受地形因子、土壤容重和質(zhì)地的影響程度最弱。

與傳統(tǒng)回歸分析相比(表2),MEMD可捕捉到樣帶1中容重、樣帶2中砂粒及樣帶3中地形濕度指數(shù)和坡度對有機質(zhì)的影響。同時,回歸方程中標準回歸系數(shù)也可表示這些影響因子在每一尺度上對有機質(zhì)預(yù)測的相對重要性(表6)。在樣帶1中,光譜主份1對所有表征尺度上土壤有機質(zhì)預(yù)測起顯著影響。土壤質(zhì)地在小尺度(≤3078 m)和大尺度(≥9535 m)處對有機質(zhì)預(yù)測起顯著影響。同時,地形因子對有機質(zhì)預(yù)測的影響隨尺度增大而增強,如坡度因子在IMF3、4、5和6處系數(shù)分為0.13、0.17、0.70和0.75。土壤容重對IMF2、3、5、6和殘差處有機質(zhì)預(yù)測都有顯著影響。因此,容重對樣帶1采樣尺度上有機質(zhì)預(yù)測的貢獻主要體現(xiàn)在這幾個尺度上的綜合作用。

在樣帶2中,由于IMF7和8中有機質(zhì)預(yù)測值與采樣尺度上實測值的相關(guān)系數(shù)不顯著,故可忽略(圖3)。從IMF1到6處,地形因子對有機質(zhì)預(yù)測的影響逐漸增強。土壤容重和質(zhì)地對不同尺度有機質(zhì)預(yù)測的影響無明顯規(guī)律。砂粒含量對采樣尺度上有機質(zhì)預(yù)測的影響主要體現(xiàn)在IMF2處。在樣帶3中,地形因子對有機質(zhì)預(yù)測的影響隨尺度的增大而增強,且在IMF6處達到最大值,然后隨尺度增大而減小,且濕度指數(shù)對樣帶3中有機質(zhì)預(yù)測的影響最為明顯。土壤質(zhì)地在其小尺度IMF1和大尺度IMF6處影響最強。容重在IMF3、4和7處對有機質(zhì)預(yù)測有穩(wěn)定影響,其余尺度不明顯。綜上,與傳統(tǒng)回歸分析相比,MEMD方法在某些表征尺度可捕捉到相關(guān)因子對有機質(zhì)分布的影響,故其預(yù)測精度高于采用傳統(tǒng)回歸分析方法。

土壤有機質(zhì)實測值和預(yù)測值的局部小波如圖4所示。與MEMD相比,逐步多元回歸預(yù)測誤差在樣帶1中主要是在位置19.8—33.0 km處局部高估,造成尺度1—2和4 km處局部方差的增大；樣帶2中,主要是在位置19.8—33.0 km處局部低估,造成尺度4—8 km處局部方差的減?。粯訋?中,主要是在位置19.8—36.3 km處局部高估,造成8 km處局部方差的增大。

圖4 土壤有機質(zhì)實測值和逐步多元回歸、MEMD法預(yù)測值的局部小波圖Fig.4 Local wavelet spectrum of measured and predicted SOM by SMLR or MEMD

3 討論

本研究采用MEMD法將太原盆地不同位置的土壤有機質(zhì)空間序列分解為不同的表征尺度,發(fā)現(xiàn)主要表征尺度在1000 m處。因此,在利用柵格數(shù)據(jù)存儲該盆地內(nèi)土壤屬性時,建議將表層土壤有機質(zhì)最佳格網(wǎng)寬度設(shè)置為1000 m,以便捕捉到有機質(zhì)較大的空間變異性。盆地上部殘差的方差百分比較大(24.19%),可能是需設(shè)計更長的采樣帶才可利用MEMD識別的更大尺度。同時,盆地上部的表征尺度與中、下部差異較大,可能是由于盆地上部距太原和晉中市等大、中型城市較近,人為活動造成土地利用破碎。此外,有機質(zhì)空間序列在垂直河流方向上的主要表征尺度沿汾河流域方向表現(xiàn)呈分散狀態(tài),表明沿河流方向的相關(guān)因子對有機質(zhì)空間分布影響的主導(dǎo)性增強。

在采樣尺度和空間多尺度關(guān)系上,盆地中部土壤有機質(zhì)分布與地形因子、容重和土壤質(zhì)地的相關(guān)性最強,這可能是由該區(qū)處于侵蝕平原(上部)與陷落盆地(下部)的交界處的特殊地質(zhì)條件造成[23]。在采樣尺度上,盆地中部有機質(zhì)與光譜主份1的相關(guān)性最弱,可能是由于中部區(qū)域土壤有機質(zhì)的變異性較小,導(dǎo)致有機質(zhì)與土壤光譜之間的相關(guān)性減弱[17]。另外,有機質(zhì)與某些因子在采樣尺度不存在顯著相關(guān),但是二者的空間多尺度關(guān)系在某些表征尺度呈顯著相關(guān),表明土壤有機質(zhì)與影響因子的空間多尺度關(guān)系能獲取更多信息。同時,基于有機質(zhì)與相關(guān)因子空間多尺度關(guān)系的有機質(zhì)預(yù)測精度顯著高于采樣尺度上直接利用逐步多元回歸分析的結(jié)果,進一步表明單一尺度(采樣尺度)的簡單分析難以揭示兩者的復(fù)雜關(guān)系。

MEMD法適用于分析非線性、非平穩(wěn)數(shù)據(jù)序列的空間多尺度關(guān)系[24]。由于該理論建立時間尚短,在土壤屬性的空間多尺度研究中還未廣泛應(yīng)用。同時,MEMD法也存在一定缺陷。Hu和Si研究結(jié)果表明,每條樣帶的所有IMF和殘差的方差百分比之和并不等于100%[25]。此外,與小波分析不同,MEMD法會針對不同采樣帶分解出不同的表征尺度,特定樣帶的研究結(jié)果不能推廣。因此,MEMD法重在分析其發(fā)生過程,在結(jié)果驗證方面存在缺陷。

4 結(jié)論

本文采用MEMD法分析了太原盆地土壤有機質(zhì)與影響因子的空間多尺度關(guān)系,主要結(jié)論如下：

(1)盆地上部主要表征尺度為1011 和1725 m,中部為982 和8573 m,下部為960 、6753 和11806 m。整個盆地內(nèi),尺度約1000 m處是土壤有機質(zhì)的主要表征尺度,且盆地內(nèi)垂直河流方向的有機質(zhì)序列主要表征尺度沿河流方向表現(xiàn)分散。

(2)土壤有機質(zhì)與影響因子在采樣尺度和MEMD空間多尺度的相關(guān)性順序為：盆地中部>下部>上部。

(3)在3種景觀樣帶上,光譜主份1與有機質(zhì)的相關(guān)性均顯著。其次,盆地上部的容重、中部的砂粒和下部的地形濕度指數(shù)對其影響較明顯,而在采樣尺度上盆地下部二者的關(guān)系并不顯著。因此,單一尺度分析不能夠全面揭示土壤有機質(zhì)與相關(guān)因子在所有空間尺度上的復(fù)雜關(guān)系,而MEMD法對有機質(zhì)的預(yù)測精度要顯著高于直接利用逐步多元回歸分析。

致謝：感謝郭穎、王鵬和徐振對樣品采集和測試提供的幫助。

[1] Bauer A, Black A L. Quantification of the effect of soil organic matter content on soil productivity. Soil Science Society of America Journal, 1994, 58(1): 185- 193.

[2] Lal R. Soil carbon sequestration impacts on global climate change and food security. Science, 2004, 304(5677): 1623- 1627.

[3] Hu K L, Wang S Y, Li H, Huang F, Li B G. Spatial scaling effects on variability of soil organic matter and total nitrogen in suburban Beijing. Geoderma, 2014, 226- 227: 54- 63.

[4] Li D F, Gao G Y, Lü Y H, Fu B J. Multi-scale variability of soil carbon and nitrogen in the middle reaches of the Heihe River basin, northwestern China. Catena, 2016, 137: 328- 339.

[5] Zhou Y, Biswas A, Ma Z Q, Lu Y L, Chen Q X, Shi Z. Revealing the scale-specific controls of soil organic matter at large scale in Northeast and North China Plain. Geoderma, 2016, 271: 71- 79.

[6] Zhu H F, Bi R T, Duan Y H, Xu Z J. Scale-location specific relations between soil nutrients and topographic factors in the Fen River Basin, Chinese Loess Plateau. Frontiers of Earth Science, 2016: 1- 10, doi: 10.1007/s11707-016-0587-y.

[7] She D L, Shao M A, Hu W, Yu S E. Variability of soil water-physical properties in a small catchment of the Loess Plateau, China. African Journal of Agricultural Research, 2010, 5(22): 3041- 3049.

[8] Zhu H F, Hu W, Bi R T, Peak D, Si B C. Scale-and location-specific relationships between soil available micronutrients and environmental factors in the Fen River basin on the Chinese Loess Plateau. Catena, 2016, 147: 764- 772.

[9] Liu Q J, Shi Z H, Fang N F, Zhu H D, Ai L. Modeling the daily suspended sediment concentration in a hyperconcentrated river on the Loess Plateau, China, using the Wavelet-ANN approach. Geomorphology, 2013, 186: 181- 190.

[10] Wang D, Fu B J, Zhao W W, Hu H F, Wang Y F. Multifractal characteristics of soil particle size distribution under different land-use types on the Loess Plateau, China. Catena, 2008, 72(1): 29- 36.

[11] Huang N E, Shen Z, Long S R, Wu M C, Shih H H, Zheng Q A, Yen N C, Tung C C, Liu H H. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proceedings of the Royal Society of London Series A: Mathematical, Physical and Engineering Sciences, 1998, 454(1971): 903- 995.

[12] She D L, Shao M A, Timm L C, Reichardt K. Temporal changes of an alfalfa succession and related soil physical properties on the Loess Plateau, China. Pesquisa Agropecuária Brasileira, 2009, 44(2): 189- 196.

[13] 劉志鵬. 黃土高原地區(qū)土壤養(yǎng)分的空間分布及其影響因素[D]. 北京: 中國科學(xué)院研究生院, 2013.

[14] Hu W, Shao M A, Si B C. Seasonal changes in surface bulk density and saturated hydraulic conductivity of natural landscapes. European Journal of Soil Science, 2012, 63(6): 820- 830.

[15] 鮑士旦. 土壤農(nóng)化分析(第三版). 北京: 中國農(nóng)業(yè)出版社, 2013.

[16] 龐獎勵, 黃春長, 賈耀峰. 不同方法測定黃土和古土壤樣品粒度的比較. 陜西師范大學(xué)學(xué)報: 自然科學(xué)版, 2003, 31(4): 87- 92.

[17] 南鋒, 朱洪芬, 畢如田. 黃土高原煤礦區(qū)復(fù)墾農(nóng)田土壤有機質(zhì)含量的高光譜預(yù)測. 中國農(nóng)業(yè)科學(xué), 2016, 49(11): 2126- 2135.

[18] Goovaerts P. Geostatistical tools for characterizing the spatial variability of microbiological and physico-chemical soil properties. Biology and Fertility of Soils, 1998, 27(4): 315- 334.

[19] Oladosu G. Identifying the oil price-macroeconomy relationship: An empirical mode decomposition analysis of US data. Energy Policy, 2009, 37(12): 5417- 5426.

[20] Rehman N, Mandic D P. Multivariate empirical mode decomposition. Proceedings of the Royal Society A, 2010, 466(2117): 1291- 1302.

[21] Ur Rehman N, Mandic D P. Filter bank property of multivariate empirical mode decomposition. IEEE Transactions on Signal Processing, 2011, 59(5): 2421- 2426.

[22] Zhao X M, Patel T H, Zuo M J. Multivariate EMD and full spectrum based condition monitoring for rotating machinery. Mechanical Systems and Signal Processing, 2012, 27: 712- 728.

[23] 姜佳奇, 莫多聞, 呂建晴, 廖奕楠, 魯鵬, 任小林. 山西太原盆地全新世地貌演化及其對古人類聚落分布的影響. 古地理學(xué)報, 2016, 18(5): 895- 904.

[24] She D L, Zheng J X, Shao M A, Timm L C, Xia Y Q. Multivariate empirical mode decomposition derived multi-scale spatial relationships between saturated hydraulic conductivity and basic soil properties. Clean-Soil, Air, Water, 2015, 43(6): 910- 918.

[25] Hu W, Si B C. Soil water prediction based on its scale-specific control using multivariate empirical mode decomposition. Geoderma, 2013, 193- 194: 180- 188.