牛文鵬 ，李青圃，李鋮，程炯，吉冬青*，蔣新宇，劉通

1.廣東省科學(xué)院生態(tài)環(huán)境與土壤研究所/華南土壤污染控制與修復(fù)國家地方聯(lián)合工程研究中心/廣東省農(nóng)業(yè)環(huán)境綜合治理重點實驗室，廣東 廣州 510650；2.珠江水利委員會珠江水利科學(xué)研究院，廣東 廣州 510611

土壤養(yǎng)分，如氮、磷和鉀是植物生長必需的基本元素，對糧食生產(chǎn)和安全有重要影響（Li et al.，2020a；劉茂秀等，2018）。準(zhǔn)確了解土壤養(yǎng)分時空分布規(guī)律是進行農(nóng)業(yè)綜合管理以及生態(tài)環(huán)境管控的重要前提（Yadav et al.，2000）。由于土壤是非均質(zhì)復(fù)雜體，土壤養(yǎng)分質(zhì)量分?jǐn)?shù)在空間上呈現(xiàn)出異質(zhì)性（Liu et al.，2013；Chatterjee et al.，2015）。這種空間異質(zhì)性是不同尺度上自然和人為活動共同作用的結(jié)果（Yadav et al.，2000；Liu et al.，2013；Liu et al.，2017；劉曉林等，2012；崔貝等，2013；崔旭輝等，2016）。相關(guān)研究表明，大氣氮沉降與有機質(zhì)分解是土壤氮的重要自然來源（李一凡等，2019）。土壤磷更易與鈣、鋁、鐵和有機質(zhì)結(jié)合，其空間分布受降雨、土地利用變化、施肥等因素的綜合影響（Werner et al.，2015）。土壤鉀的固定與釋放主要受到土壤質(zhì)地、粘土礦物質(zhì)量分?jǐn)?shù)、土壤酸堿度和濕度等影響（Bruckert et al.，1992）。不同的影響因子具有不同的作用范圍，有些因子（如土壤理化性質(zhì)、養(yǎng)分添加）可能在局部尺度起決定性作用，而有些因子（如氣候和成土母質(zhì)）則在較大的尺度起決定性作用（Mage et al.，2013；Lv et al.，2014；王現(xiàn)潔等，2020）。不同影響因子在不同尺度作用強度不同且交互作用，導(dǎo)致我們觀察到的土壤養(yǎng)分空間格局混合了多個尺度的信息。剝離出不同尺度土壤養(yǎng)分的空間變異規(guī)律及其關(guān)鍵驅(qū)動因子對于土壤養(yǎng)分管理至關(guān)重要（崔貝等，2013）。

由于人力、物力資源的限制，在區(qū)域尺度開展詳細土壤養(yǎng)分調(diào)查費時費力。地理信息系統(tǒng)、地統(tǒng)計學(xué)（普通克里格、指示克里格、協(xié)同克里格等）（Eldeiry et al.，2010；劉曉林等，2012；崔旭輝等，2016）和傳統(tǒng)統(tǒng)計學(xué)方法（如相關(guān)分析和回歸分析）（Huang et al.，2006；Mage et al.，2013）的發(fā)展為揭示區(qū)域尺度上土壤養(yǎng)分空間變異及其驅(qū)動力提供了便利。但是，以上方法多著眼于單一尺度，不能有效區(qū)分或剝離研究要素在特定尺度的空間分異規(guī)律（Lin et al.，2006；Shan et al.，2013；Lv et al.，2014）。因子克里格（Factorial Kriging Analysis，F(xiàn)KA），整合了地統(tǒng)計學(xué)（協(xié)同克里格）、傳統(tǒng)統(tǒng)計學(xué)方法（主成分分析）和協(xié)同區(qū)域化線性模型，可以識別出空間變量的空間結(jié)構(gòu)變異特征，有效過濾掉不同尺度信息的疊加效應(yīng)，在揭示土壤屬性和元素多尺度空間變異和驅(qū)動力方面具有很大潛力（Goovaerts，1992；Imrie et al.，2008；呂真真等，2014）。因子克里格方法被廣泛應(yīng)用于探究土壤性質(zhì)的多尺度空間變異特征及其驅(qū)動力分析（Liu et al.，2013；Zhou et al.，2016；Du et al.，2017）。

當(dāng)前對于土壤養(yǎng)分分布特征及其驅(qū)動力的研究較多，但多集中在單一尺度（Rodríguez et al.，2011；Song et al.，2019；鐘繼洪等，2009；劉曉林等，2012；安永龍等，2018），多尺度的研究相對較少。小尺度的研究有助于精準(zhǔn)施肥，而大尺度上的研究對于區(qū)域平衡施肥和種植區(qū)劃至關(guān)重要（呂真真等，2014）。在制定土壤養(yǎng)分管控策略時，往往需要考慮不同尺度土壤養(yǎng)分變異特征及其驅(qū)動因素。Liu et al.（2013）使用因子克里格和逐步回歸分析定量探究了土壤養(yǎng)分的多尺度空間變異及其環(huán)境驅(qū)動力，發(fā)現(xiàn)土壤養(yǎng)分的空間變異在局部尺度（12 km）主要受控于土地利用類型或海拔，而在區(qū)域尺度（84 km）主要受控于成土母質(zhì)。但是，該研究在進行驅(qū)動力分析時沒有考慮自變量（即環(huán)境變量）的空間變異。事實上，任何地理變量都存在空間變異，在驅(qū)動力分析時要同時考慮土壤養(yǎng)分和環(huán)境變量的空間變異特征和尺度。鑒于此，本研究以典型農(nóng)業(yè)聚集區(qū)-珠江三角洲為例，基于417個表層土壤樣品和6大類環(huán)境變量（氣候、成土母質(zhì)、土壤理化性質(zhì)、地形/植被覆蓋、距污染源的距離、社會經(jīng)濟），綜合考慮自變量和因變量的空間變異特征，使用因子克里格法探討土壤養(yǎng)分的多尺度空間分異規(guī)律及其環(huán)境驅(qū)動力，以期為區(qū)域土壤養(yǎng)分管理提供參考。

1 研究數(shù)據(jù)與方法

1.1 研究區(qū)概況

珠江三角洲位于廣東省中南部，毗鄰港澳，是粵港澳大灣區(qū)的重要組成部分（圖 1）。該地區(qū)地貌主要為三角洲平原與低山丘陵構(gòu)成，范圍包括廣州、佛山、珠海、深圳、江門、東莞、和中山7個城市以及惠州和肇慶部分區(qū)域。珠江三角洲地區(qū)為典型的亞熱帶季風(fēng)氣候，全年溫暖濕潤，年均降水1900 mm，年均氣溫22 ℃。研究區(qū)成土母質(zhì)主要包括（河/海）沖積物、花崗巖和砂頁巖，分別占研究區(qū)總面積的38%、30%和22%（Li et al.，2020b），主要土壤類型包括赤紅壤、河流沖積土、紅壤和人工堆疊土等。研究區(qū)陸域面積4.12×104km2，土地利用類型以林地為主（占比43%），其次為城市用地和水體，分別占研究區(qū)面積的23%和11%（根據(jù)2017年土地利用/覆被圖統(tǒng)計；圖1B）。珠江三角洲人口密集，是廣東省和全國經(jīng)濟發(fā)展的中心區(qū)域，根據(jù)廣東省統(tǒng)計年鑒，該區(qū)域2018年常住人口為6301萬人，GDP 81048億元，分別占全省總?cè)丝诤虶DP的55.5%和80.2%（廣東省統(tǒng)計局，2019）。珠江三角洲是重要的農(nóng)業(yè)聚集區(qū)，農(nóng)業(yè)用地（主要為耕地和園地）約占17%。隨著工業(yè)化和城市化進程加快，研究區(qū)農(nóng)用地迅速減少，土壤質(zhì)量也呈現(xiàn)下降趨勢（侯紅乾等，2013；Li et al.，2020b）。

圖1 珠江三角洲位置（A）、土地利用和覆被圖（B）和土壤樣點采集（C）Fig.1 The location (A), land use and cover map (B) and soil sampling (C) of Pearl River Delta

1.2 土壤樣點采集與養(yǎng)分質(zhì)量分?jǐn)?shù)分析

2018年，基于雙密度網(wǎng)格隨機取樣法對珠江三角洲地區(qū)表層土壤（0—20 cm）進行采樣。其中，中心城區(qū)每5 km×5 km網(wǎng)格隨機布設(shè)1個樣點，城郊和農(nóng)村每10 km×10 km網(wǎng)格隨機布設(shè)1個樣點，共計417個樣點（圖1C）。每個采樣點，由5—8個土樣混合而成，采用四分法取1 kg帶回實驗室分析。土壤樣品經(jīng)自然風(fēng)干，挑出石礫和植物殘體，分別過10目和100目尼龍篩備用。采用凱氏定氮法測定總氮（N），氫氧化鈉高溫熔融后、分別使用鉬銻抗比色法和火焰光度計法測定總磷（P）和總鉀（K），乙酸銨提取、火焰光度計法測定速效鉀（AK）（NT/T 395—2012）。堿解-擴散法測定堿解氮（AN）（LYT 1229—1999）；鹽酸?氟化銨法提取、鉬銻抗比色法測定有效磷（AP）（NY/T 1121.7—2014）。土壤 pH采用 1∶2.5的土水質(zhì)量比提取、電位法測定，土壤有機質(zhì)和質(zhì)地分別采用重鉻酸鉀氧化?容量法和比重計法測定（NT/T 395—2012）。

1.3 數(shù)據(jù)處理與分析

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

基于SPSS 20統(tǒng)計軟件，采用平均值±3倍標(biāo)準(zhǔn)差法去除特異值，最終獲得土壤養(yǎng)分樣本395個。前人研究發(fā)現(xiàn)，氣候（Li et al.，2021）、成土母質(zhì)（Liu et al.，2013；崔旭輝等，2016）、地形（李超等，2014）、土壤理化性質(zhì)（Baxter et al.，2003）、植被（Xu et al.，2020）、礦物組成（Li et al.，2021）、肥料添加（胡艷霞等，2013；鄭昊楠等，2019）等是影響土壤養(yǎng)分空間變異的主要因素?；谝陨衔墨I和數(shù)據(jù)可獲得性，本研究選取氣候、成土母質(zhì)、土壤理化性質(zhì)、地形/植被覆蓋、距污染源距離和社會經(jīng)濟6大類共計16個環(huán)境變量作為自變量（表1）。這里，我們使用距污染源距離和社會經(jīng)濟變量來替代肥料添加等人類干擾強度數(shù)據(jù)，因為該數(shù)據(jù)較難獲得。由于土壤養(yǎng)分和環(huán)境變量多不符合正態(tài)分布，考慮到異常值和數(shù)據(jù)分布會影響克里格插值精度，首先將原始數(shù)據(jù)進行高斯變換，使其滿足正態(tài)分布（平均值為0，標(biāo)準(zhǔn)差為1）。分別采用Pearson相關(guān)分析法和因子克里格法探討不同環(huán)境變量對土壤養(yǎng)分空間分布的影響。

表1 環(huán)境變量指標(biāo)及數(shù)據(jù)來源Table 1 Environmental variables and data sources

1.3.2 因子克里格法

尺度是研究某一現(xiàn)象或過程所采用的空間或時間單位。不同學(xué)科對尺度的定義不同，表達也不同，如地圖學(xué)的制圖比例尺、遙感學(xué)的空間分辨率、景觀生態(tài)學(xué)的粒度（最小辨識單元）和幅度（研究范圍大?。╓u，2004）。本研究中的尺度對應(yīng)于空間分析生態(tài)學(xué)的研究步長，類似于景觀生態(tài)學(xué)的粒度，為最小數(shù)據(jù)分析單元，如10—20 km往往對應(yīng)于田塊或街區(qū)尺度，而100 km甚至更大步長則對應(yīng)于城市或區(qū)域尺度。

因子克里格，可用來揭示變量的多尺度空間變異特征及其來源（Lv et al，2013；Lv et al.，2014）。其主要包括以下步驟：（1）基于方差和交叉方差模型識別變量的空間變異尺度；（2）區(qū)域化變量的線性模型（Linear Model of Co-regionalization，LMC），即把變量的總體分解為不同空間尺度上兩個及以上協(xié)變量的加權(quán)線性函數(shù)（公式1）；（3）基于主成分分析探討變量在不同尺度間的相關(guān)關(guān)系；（4）使用協(xié)同克里格法估算變量的區(qū)域化因子及空間組分。

線性協(xié)同區(qū)域化模型通常由 3個空間要素組成，分別是：塊金尺度（Nugget effect）、局部尺度（短變程尺度）和區(qū)域尺度（長變程尺度）。具體公式如下（Imrie et al.，2008）：

式中，Γ(h)為m×m階變異函數(shù)矩陣；Ns為空間尺度的個數(shù)；Bu為協(xié)同區(qū)域化矩陣，用于描述變量與特定空間尺度的關(guān)系。gu(h)為方差函數(shù)，h為步長（km）。

對于特定尺度，土壤污染指數(shù)與環(huán)境變量的相關(guān)系數(shù)ρij可用如下公式表示（Lv et al.，2014）：

式中，qij、λj和σi分別代表特定尺度主成分分析的特征向量、特征值和第i個變量的方差。

為了驗證克里格插值結(jié)果的準(zhǔn)確性，采用留一法進行交叉檢驗，即每次把某一樣點作為測試數(shù)據(jù)，其余樣點作為訓(xùn)練數(shù)據(jù)，將模擬結(jié)果與真實值比較。當(dāng)模擬結(jié)果的平均誤差（Mean Error）越接近于0，平均標(biāo)準(zhǔn)誤方差（Mean of Squared Standardized Error）的值越接近于1時，表明模型擬合效果越好。使用Geovariances Isatis軟件，對有效樣點（n=395）的養(yǎng)分質(zhì)量分?jǐn)?shù)進行因子克里格分析及交叉檢驗。

2 結(jié)果與分析

2.1 土壤養(yǎng)分質(zhì)量分?jǐn)?shù)特征及其與環(huán)境變量的關(guān)系

珠江三角洲地區(qū)土壤總氮、總磷、總鉀、堿解氮和速效鉀質(zhì)量分?jǐn)?shù)分別為 (1.2±0.5) g·kg?1（平均值±標(biāo)準(zhǔn)差）、(1.1±0.8) g·kg?1、(16.8±7.8) g·kg?1、(115.3±41.8) mg·kg?1和 (184.6±144.1) mg·kg?1，變異系數(shù)介于 44%—78%；土壤有效磷的質(zhì)量分?jǐn)?shù)為(161.4±174.1) mg·kg?1，變異系數(shù)為 108%。通常認(rèn)為，變異系數(shù)＜10%為弱強度變異，介于10%—100%之間為中等強度變異，＞100%為強度變異（表2）。由此可知，土壤養(yǎng)分除有效磷為強度變異，其他元素均為中等強度變異。對比全國第二次土壤普查養(yǎng)分分級標(biāo)準(zhǔn)，59%—80%的土壤樣品土壤養(yǎng)分等級為中等以上；土壤磷和有效磷尤為豐富，處于很豐富等級的樣品分別占總樣品數(shù)的46%和63%。從不同土地利用類型來看，土壤養(yǎng)分質(zhì)量分?jǐn)?shù)普遍呈現(xiàn)出耕地＞園地＞林地/城市綠地，且耕地土壤養(yǎng)分（總鉀除外）質(zhì)量分?jǐn)?shù)顯著高于其他三類土壤（P＜0.05）。而園地土壤的堿解氮、有效磷和速效鉀顯著高于林地和城市綠地土壤，林地和城市綠地的養(yǎng)分質(zhì)量分?jǐn)?shù)無顯著差異（P＜0.05）（圖2）。

表2 土壤養(yǎng)分質(zhì)量分?jǐn)?shù)的統(tǒng)計值Table 2 Statistics of soil nutrients

圖2 不同土地利用類型土壤養(yǎng)分質(zhì)量分?jǐn)?shù)的比較Fig.2 Comparisons of soil nutrient contents among different land use types

土壤養(yǎng)分質(zhì)量分?jǐn)?shù)與環(huán)境變量的相關(guān)關(guān)系分析顯示（表 3），土壤養(yǎng)分質(zhì)量分?jǐn)?shù)普遍與土壤理化性質(zhì)（pH、有機質(zhì)質(zhì)量分?jǐn)?shù) SOM、土壤質(zhì)地Sand*、Silt、Clay）、成土母質(zhì)（沉積物Depo、花崗巖Gran*和砂頁巖Sandsh*）、氣候因子（降雨P(guān)rec和溫度 Temp）、距離水體的距離（Wt*）和海拔（DEM*）呈顯著相關(guān)（P＜0.05；表3）。鉀和有效鉀的質(zhì)量分?jǐn)?shù)與社會經(jīng)濟變量（GDP*、Pop*）和距離工廠距離（Ind*）呈顯著相關(guān)（P＜0.05）。但是，總體上土壤養(yǎng)分與環(huán)境變量的相關(guān)關(guān)系相對較弱（相關(guān)系數(shù)＜0.6，SOM除外）。

表3 土壤養(yǎng)分質(zhì)量分?jǐn)?shù)與環(huán)境變量的Pearson相關(guān)系數(shù)Table 3 Pearson correlation coefficients between soil nutrient contents and environmental variables

2.2 土壤養(yǎng)分多尺度空間變異及其環(huán)境驅(qū)動力

2.2.1 空間變異尺度識別及模型精度檢驗

土壤養(yǎng)分的方差及交叉方差變異圖顯示（圖3），土壤養(yǎng)分的空間變異均存在塊金效應(yīng)，即存在采樣和測量誤差，主要由于實驗誤差和小于實際取樣尺度引起的變異。隨著距離增加，方差值增加或減少并在26 km附近趨于平緩，之后緩慢變化并在90 km處相對穩(wěn)定，接近基臺值。經(jīng)試驗，雙球面模型能夠較好地擬合土壤養(yǎng)分的空間變異（衡量標(biāo)準(zhǔn)：殘差平方和與Akaike Criterion較小）。總體上，本研究的協(xié)同區(qū)域化模型可分解為3部分：塊金值、局部尺度的區(qū)域化因子（球面模型；26 km）和區(qū)域尺度的區(qū)域化因子（球面模型；90 km）。經(jīng)交叉檢驗，土壤養(yǎng)分線性協(xié)同區(qū)域化模擬的平均誤差介于?0.011—0.027，接近于理想值 0，而平均標(biāo)準(zhǔn)誤方差均接近于理想值1（0.99—1.24），說明雙球面模型能夠較好地模擬土壤養(yǎng)分的空間變異特征。

圖3 土壤養(yǎng)分方差函數(shù)和交叉方差變異曲線（土壤養(yǎng)分與環(huán)境變量的方差函數(shù)未顯示）Fig.3 Variogram and cross-variogram of soil nutrients (not shown for those between soil nutrients and environmental variables)

2.2.2 土壤養(yǎng)分與環(huán)境變量的多尺度環(huán)境驅(qū)動力

不同空間尺度相關(guān)分析結(jié)果顯示（表 4），土壤養(yǎng)分與環(huán)境變量的關(guān)系存在相互加強和抵消的效應(yīng)。多數(shù)土壤養(yǎng)分質(zhì)量分?jǐn)?shù)與土壤有機質(zhì)（SOM）、土壤質(zhì)地（Silt、Clay）、沉積物（Depo）在3個空間尺度上均呈正相關(guān)，與土壤沙粒質(zhì)量分?jǐn)?shù)（與Sand*正相關(guān)）、距離污染源的距離（與Wt*正相關(guān)）和海拔（與DEM*正相關(guān)）在3個空間尺度上均呈負(fù)相關(guān)，而與其他環(huán)境變量的關(guān)系則存在尺度效應(yīng)。多數(shù)土壤養(yǎng)分質(zhì)量分?jǐn)?shù)與土壤pH、花崗巖（Gran*）、溫度（Temp）和NDVI在塊金尺度呈負(fù)相關(guān)，在短變程和長變程尺度呈正相關(guān)。類似地，土壤養(yǎng)分質(zhì)量分?jǐn)?shù)與降雨（Prec）在塊金和短變程尺度呈負(fù)相關(guān)，在長變程尺度呈正相關(guān)；土壤養(yǎng)分質(zhì)量分?jǐn)?shù)與砂頁巖（Sandsh*）在塊金尺度和短變程尺度呈正相關(guān)，在長變程尺度呈負(fù)相關(guān)；而與GDP在短變程尺度呈負(fù)相關(guān)，在塊金尺度和長變程尺度呈正相關(guān)。

表4 不同空間尺度土壤養(yǎng)分質(zhì)量分?jǐn)?shù)與環(huán)境變量的相關(guān)系數(shù)Table 4 Correlation coefficients between soil nutrient contents and environmental variables at different spatial scales

為了更好地顯示多空間尺度變量的相互關(guān)系，我們把土壤養(yǎng)分空間組分與主成分得分的相關(guān)系數(shù)投影在單元環(huán)中（圖 4）。塊金尺度，前兩個主成分可解釋該尺度變量空間變異的45.1%，其中第一主成分（PC1）和第二主成分（PC2）分別揭示33.5%和 11.6%的變異。PC1與土壤質(zhì)地（SOM、Silt）和土壤養(yǎng)分質(zhì)量分?jǐn)?shù)（K除外）相關(guān)系數(shù)較大，表明該主成分主要受到這些因子的影響。PC2與Sand*和Clay相關(guān)關(guān)系較密切，主要受到這兩種土壤理化性質(zhì)的影響（表4，圖4A）。從變量與變量之間的距離來看，土壤養(yǎng)分質(zhì)量分?jǐn)?shù)與有機質(zhì)質(zhì)量分?jǐn)?shù)關(guān)系尤為密切。短變程尺度，前兩個主成分可解釋該尺度變量空間變異的59.2%。根據(jù)變量與主成分（圖 4）以及因變量與自變量之間的相關(guān)系數(shù)（表4，圖4B），土壤養(yǎng)分可劃分為4組：（1）土壤總氮，主要受控于土壤理化性質(zhì)（SOM、Sand*、Silt和Clay），其次為距離水體的距離和溫度；（2）土壤總鉀，主要受成土母質(zhì)（Sandsh*）的影響；（3）土壤有效磷，主要受到人口密度（Pop*）的影響；（4）土壤總磷、堿解氮和有效鉀，主要受控于海拔（DEM）、距離工廠的距離（Ind*）和社會經(jīng)濟因素（GDP*或Pop*）。此外，總磷質(zhì)量分?jǐn)?shù)還與土壤pH關(guān)系密切。類似地，在長變程尺度，前兩個主成分可解釋該尺度變量空間變異的67.4%。根據(jù)環(huán)境驅(qū)動的不同，土壤養(yǎng)分可劃分為3組（表4，圖4C）：（1）土壤總氮、總磷、堿解氮和有效磷，主要受控于氣候因素（Prec、Temp），其次為距離水體的距離和土壤有機質(zhì)質(zhì)量分?jǐn)?shù)（圖4C；表4）?？偟杂胁煌?，其主要受控于土壤理化性質(zhì)（SOM、Sand*和Clay）的影響，氣候因素的作用位于其次。此外，總磷還受到成土母質(zhì)（Depo）的影響。（2）總鉀，主要受成土母質(zhì)（Sandsh*）與距污染源距離（Ind*）影響，其次為社會經(jīng)濟因素（GDP*、Pop*）。（3）有效鉀，其主要受到氣候因素的影響（表4）。

圖4 塊金（A）、短變程（B）和長變程（C）尺度土壤養(yǎng)分空間組分與主成分得分相關(guān)關(guān)系的單位投影圓Fig.4 Projections of correlations between spatial components for soil nutrient contents and the principal component scores into unit circles at nugget (A),short-range (B) and long-range (C) scales

2.2.3 土壤養(yǎng)分的多尺度空間結(jié)構(gòu)

圖5和圖6顯示，短變程和長變程尺度上土壤養(yǎng)分空間分異明顯。短變程尺度，土壤養(yǎng)分質(zhì)量分?jǐn)?shù)高的區(qū)域較為分散，主要位于研究區(qū)的中部，邊緣地區(qū)質(zhì)量分?jǐn)?shù)較低（圖 5）。土壤總氮和總磷質(zhì)量分?jǐn)?shù)空間變異范圍較小，數(shù)據(jù)分布范圍主要集中在?0.5—1.0（圖5A、B）。土壤堿解氮和有效鉀的空間分布較為相似，質(zhì)量分?jǐn)?shù)高的區(qū)域主要散布于廣州市、佛山市、珠海、惠州中西部和東莞西部（圖5D、F）?？傗浐陀行Я踪|(zhì)量分?jǐn)?shù)較高的區(qū)域也基本集中在以上城市，但是集聚位置又略有不同（圖5C、E）。在長變程尺度，土壤養(yǎng)分質(zhì)量分?jǐn)?shù)空間變異范圍較短變程尺度減小，數(shù)據(jù)變異集中在?0.5—0.5，沒有明顯的高值區(qū)域。土壤總氮、總磷、總鉀、堿解氮、有效磷和速效鉀空間變異規(guī)律較為接近，質(zhì)量分?jǐn)?shù)較高的區(qū)域（黃色區(qū)域）主要位于廣州北部與中部、惠州中部、珠海北部、佛山南部和江門中部（圖 6）?？傮w來說，劃分為同一組的土壤養(yǎng)分元素空間分異特征較為相似（圖5、6）。長變程尺度，土壤養(yǎng)分質(zhì)量分?jǐn)?shù)的空間異質(zhì)性明顯小于短變程尺度，前者的數(shù)值分布范圍明顯較小。

圖5 短變程尺度土壤養(yǎng)分空間組分的分布Fig.5 Distributions for spatial components of soil nutrients at short-range scale

圖6 長變程尺度土壤養(yǎng)分空間組分的分布Fig.6 Distribution for spatial components of soil nutrients at long-range scale

3 討論

3.1 土壤養(yǎng)分多尺度空間變異

珠江三角洲地區(qū)土壤養(yǎng)分質(zhì)量分?jǐn)?shù)存在空間變異，多數(shù)屬于中等變異（有效磷除外，為強度變異）。根據(jù)《全國第二次土壤普查》養(yǎng)分分級標(biāo)準(zhǔn)，研究區(qū)養(yǎng)分質(zhì)量分?jǐn)?shù)普遍較高，多數(shù)土壤樣品養(yǎng)分狀態(tài)處于中等偏上水平（等級≥豐富），尤其磷和有效磷質(zhì)量分?jǐn)?shù)尤為豐富。研究區(qū)屬南方酸性紅壤，富含鐵、鋁離子與施入土壤的磷肥發(fā)生化學(xué)反應(yīng)并形成沉淀累積在土壤中，不能被植物利用，導(dǎo)致土壤磷質(zhì)量分?jǐn)?shù)普遍較高（王現(xiàn)潔等，2020）。從不同土地利用類型來看，土壤養(yǎng)分質(zhì)量分?jǐn)?shù)普遍呈現(xiàn)出耕地＞園地＞林地/城市綠地。其他研究也發(fā)現(xiàn)類似規(guī)律，耕地和果園過量施肥可能是導(dǎo)致其養(yǎng)分質(zhì)量分?jǐn)?shù)較高的主要原因（侯紅乾等，2013；胡艷霞等，2013）。中國的肥料利用率普遍較低，氮、磷和鉀肥利用率普遍不足30%（于飛等，2015；王榮萍等，2016；王現(xiàn)潔等，2020）。過量的養(yǎng)分可深層淋溶進入地下水，或進入地表水導(dǎo)致富營養(yǎng)化，造成環(huán)境污染風(fēng)險。廣東省蔬菜地土壤氮磷質(zhì)量分?jǐn)?shù)均較高，具有很大的流失風(fēng)險（張永起等，2010）。

因子克里格結(jié)果顯示，土壤養(yǎng)分空間異質(zhì)性隨尺度增大而減小。雖然養(yǎng)分高值區(qū)主要位于研究區(qū)中部（廣州、佛山、珠海、惠州中西部和東莞西部），但是土壤養(yǎng)分的空間變異性在短變程尺度較長變程尺度更強，高值分布區(qū)更為分散。類似地，其他針對土壤性質(zhì)空間分異規(guī)律的研究也發(fā)現(xiàn)小尺度上的異質(zhì)性高于大尺度（Lv et al.，2014；Zhou et al.，2016；Du et al.，2017）。土壤養(yǎng)分來源、利用效率及其在土壤中的遷移是決定土壤養(yǎng)分空間變異的主要因素（安永龍等，2018）。例如，不同類型的有機肥和土壤母質(zhì)對土壤磷的固定能力不盡相同，導(dǎo)致土壤有效磷質(zhì)量分?jǐn)?shù)存在差異（王榮萍等，2014；余煒敏等，2018a；余煒敏等，2018b）。小尺度上土壤養(yǎng)分分布異質(zhì)性高，可能與其更多受到隨機因素（如化肥農(nóng)藥的添加）影響有關(guān)。此外，研究區(qū)土壤有效態(tài)形式的養(yǎng)分空間分異性高于全態(tài)形式，這與很多研究結(jié)果類似（李旭等，2012；安永龍等，2018）。這可能由于不同植物肥料利用效率不同，導(dǎo)致殘留在土壤中的養(yǎng)分在空間上存在較大差異（王榮萍等，2016；王現(xiàn)潔等，2020）。

3.2 土壤養(yǎng)分多尺度驅(qū)動機制

土壤養(yǎng)分的空間變異，是多個尺度上自然（氣候、成土母質(zhì)、地形、土壤理化性質(zhì)等）和人為因素（特別是人為活動，如施肥、耕作方式、種植制度等）共同作用的結(jié)果（曹宏杰等，2013；李超等，2014；崔旭輝等，2016；安永龍等，2018）。土壤養(yǎng)分與環(huán)境變量的關(guān)系在不同尺度存在增強和抵消的情況（表 4），導(dǎo)致傳統(tǒng)的相關(guān)分析等統(tǒng)計方法揭示出來的兩者之間的關(guān)系相對較弱（表 3）。因子克里格法可以有效過濾掉不同尺度的交互影響及干擾，有效揭示土壤元素與環(huán)境變量的相互關(guān)系，在本研究及其他研究中得到證明（Goovaerts，1992；Imrie et al.，2008；Liu et al.，2013；Lv et al.，2014）。

本研究結(jié)果顯示，土壤養(yǎng)分的空間變異及其關(guān)鍵環(huán)境驅(qū)動因子因尺度而異。在短變程尺度，土壤養(yǎng)分質(zhì)量分?jǐn)?shù)空間異質(zhì)性較高，高值分布區(qū)較為分散，多數(shù)土壤養(yǎng)分尤其是有效態(tài)養(yǎng)分主要受到人為活動影響（如距離工廠距離、GDP和人口密度）和海拔的影響。土壤養(yǎng)分質(zhì)量分?jǐn)?shù)與DEM呈負(fù)相關(guān)（與DEM*呈正相關(guān)），與Liu et al.（2013）的研究結(jié)果一致。表明在地形、重力和徑流侵蝕等因素的綜合作用下，土壤養(yǎng)分從高海拔區(qū)運移并累積至低海拔區(qū)。土壤總磷、堿解氮和有效鉀與社會經(jīng)濟因素呈負(fù)相關(guān)（與GDP*和Pop*呈正相關(guān)），這可能源于經(jīng)濟發(fā)達和人口密集區(qū)建設(shè)用地比重較大，對土壤干擾大，導(dǎo)致土壤較貧瘠。此外，土壤總氮質(zhì)量分?jǐn)?shù)受控于土壤理化性質(zhì)、距離水體距離和溫度，這可能與粘質(zhì)土（粘粒和粉粒質(zhì)量分?jǐn)?shù)較高）有利于截留土壤養(yǎng)分、面源污染、溫度的升高有利于促進成土母質(zhì)礦化作用和養(yǎng)分循環(huán)有關(guān)（Rodríguez et al.，2011；Sun et al.，2017）。土壤有效磷與土壤有機質(zhì)質(zhì)量分?jǐn)?shù)、粘粒質(zhì)量分?jǐn)?shù)、粉粒質(zhì)量分?jǐn)?shù)和氣溫呈負(fù)相關(guān)，與前人研究結(jié)論一致（Mage et al.，2013；Yu et al.，2016）。長變程尺度，土壤養(yǎng)分質(zhì)量分?jǐn)?shù)更多受控于自然因素（氣候和成土母質(zhì)）的影響，這與前人研究結(jié)果類似（Mazzarino et al.，1998；Tang et al.，2016；王幼奇等，2014）。有學(xué)者發(fā)現(xiàn)，沉積巖發(fā)育的土壤總磷質(zhì)量分?jǐn)?shù)較蛇紋巖發(fā)育的土壤高2—4倍（Kitayama et al.，2000）；碳酸鹽巖、輝長巖和石英閃長巖發(fā)育的土壤固鉀能力較強（Portela et al.，2019）。雖然成土母質(zhì)不是土壤氮的主要來源（大氣是主要來源），但是其會影響土壤中氮素的累積和分解；沉積物、石灰?guī)r、泥巖和玄武巖等母質(zhì)發(fā)育的粘土具有較高的有機氮儲存能力（Zhong et al.，2019）。此外，Li et al.（2021）發(fā)現(xiàn)氣候是決定土壤鉀變異的最重要因素，其次為地形條件、粘土礦物組成和關(guān)鍵土壤理化性質(zhì)。Yu et al.（2016）發(fā)現(xiàn)CO2和溫度升高能夠促進土壤中磷的釋放/流失。Mazzarino et al.（1998）和脫云飛等（2020）發(fā)現(xiàn)，降雨的增加可以加快氮、磷運移速度，促進土壤養(yǎng)分的礦化過程?？傮w上，土壤養(yǎng)分的關(guān)鍵影響因素因尺度而異，這在前人研究中也有發(fā)現(xiàn)。例如，Liu et al.（2013）發(fā)現(xiàn)大尺度上土壤總磷和總氮分布受控于成土母質(zhì)，小尺度上主要受土地利用類型或地形因素影響。這種尺度效應(yīng)意味著管理者需要在局部尺度（田塊或街區(qū)尺度）、城市或區(qū)域尺度制定不同的調(diào)控措施，促進土壤養(yǎng)分的有效利用?；凇疤卣?機制-管理”思路，本研究提出土壤養(yǎng)分多尺度管理概念模型，根據(jù)不同尺度肥力區(qū)劃及其驅(qū)動力制定不同的養(yǎng)分管理策略（圖7）。具體來說，局部尺度，針對土壤高肥力區(qū)，應(yīng)避免過量施肥，可采用減量施肥、測土配方施肥、施用有機替代肥（如生物有機肥、秸稈還田）等方式減少化肥添加（鄭昊楠等，2019；胡貴貴等，2020；譚海燕等，2021；董佳琦等，2021）。在土壤高肥力且海拔與坡度較高的區(qū)域，采取免耕、等高壟作等措施以防治水土流失，減少面源污染（張少良等，2020）。而針對部分土壤低肥力區(qū)，應(yīng)根據(jù)養(yǎng)分狀況調(diào)整耕作方式，制定合理施肥措施來提升土壤肥力（曹宏杰等，2013；Shang et al.，2014；Krauss et al.，2020）。區(qū)域尺度，需重視氣候及其成土母質(zhì)對土壤養(yǎng)分的影響。在高肥力且降雨集中區(qū)，需重點防范氮磷隨降雨遷移造成的面源污染，應(yīng)設(shè)置長期定位觀測點加大對該區(qū)域面源污染的監(jiān)控。在低肥力且花崗巖廣泛分布區(qū)，應(yīng)采取施用土壤改良劑等措施提升土壤養(yǎng)分固存能力；對于其他低肥力區(qū)，應(yīng)根據(jù)養(yǎng)分狀況調(diào)整農(nóng)業(yè)種植結(jié)構(gòu)和耕作方式，改變區(qū)域養(yǎng)分空間分布的不均衡性。

圖7 基于“特征-機制-管控”思路的土壤養(yǎng)分多尺度管理概念模型Fig.7 A multi-scale conceptual model for soil nutrient management based on the principle of “soil characteristics-processes-managements”

4 結(jié)論

本研究以珠江三角洲為例，采用多元因子克里格法分析了研究區(qū)6種土壤養(yǎng)分的多尺度空間結(jié)構(gòu)及其環(huán)境驅(qū)動力。主要結(jié)論如下：

（1）對比全國第二次土壤普查養(yǎng)分分級標(biāo)準(zhǔn)，研究區(qū)土壤養(yǎng)分等級為中等及以上，磷和有效磷尤為豐富。不同土地利用類型土壤養(yǎng)分質(zhì)量分?jǐn)?shù)不盡相同，按照從高到低排序為：耕地＞園地＞林地/城市綠地。

（2）研究區(qū)土壤養(yǎng)分質(zhì)量分?jǐn)?shù)均為中等變異，除了有效磷為強度變異。有效態(tài)形式的土壤養(yǎng)分空間變異性明顯大于全態(tài)形式。

（3）因子克里格法將土壤養(yǎng)分的空間變異分解為塊金效應(yīng)、短變程尺度變異和長變程尺度的變異。不同空間尺度，土壤養(yǎng)分的分布格局不盡一致。短變程尺度養(yǎng)分空間變異性更強，異質(zhì)性更高。此外，土壤養(yǎng)分的關(guān)鍵環(huán)境驅(qū)動因子因尺度而異：塊金尺度和短變程尺度，主要受地形和局部人類活動的影響，而長變程尺度主要受氣候和成土母質(zhì)影響。因此，未來應(yīng)在不同空間尺度制定不同的土壤養(yǎng)分管理策略增加養(yǎng)分利用效率。