基于指示克里格的青銅峽灌區土壤鹽漬化風險與潛水埋深和礦化度關系研究

趙志鵬，黃小琴，方 磊

（寧夏回族自治區水文環境地質調查院，寧夏 銀川 750011）

0 引 言

土壤是世界上最復雜多樣的生態系統［1］。土壤鹽漬化是可溶性鹽分在土壤中積累致使土壤特性變化和質量下降的過程［2］。我國鹽漬土面積達到3.69×107hm2，不僅是制約農業可持續發展的主要障礙，也可能造成經濟損失、次生危害和生態環境脆弱［3，4］。在任何引水灌溉系統的干旱半干旱地區都會存在因鹽漬化引起的土地退化問題［5］，與這些區域地勢低平，排水不暢，地下水埋藏較淺關系很大［6］。地下水作為鹽分傳輸、累積和排泄的主要載體，與土壤鹽分關系密切，特別是潛水位埋深和礦化度對土壤鹽分影響顯著［7］。因此，在綜合考慮土壤鹽漬化與地下水關系的基礎上，分析評價土壤鹽漬化風險，對農業生產、水土資源利用和環境保護等方面均有重要意義。但在進行區域性土壤、地下水等特性的分析中經常會出現樣本數據中存在比算術平均值或中位數要大（或小）得多的特異值，從而使樣本偏離正態（或對數正態）分布，影響變異函數的穩健性。指示克立格是一種非參數估計方法，在處理有偏數據的明顯優勢，使其在生態學［8，9］、地質學［10-12］、環境科學［13-15］等領域被廣泛應用。采用指示克里格方法開展土壤鹽漬化風險評價的研究較多，周在明在環渤海低平原區的研究表明地下水礦化度大于2 g/L，地下水位埋深小于3 m 以及土壤鹽分含量大于1 g/kg 的概率分布呈一致性［16］；姚榮江等的研究證明在黃河三角洲典型地區地下水埋深與土壤鹽分的概率分布存在空間上的規律性與相似性［17］；李仙岳等利用指示克里格法比較了永濟灌域春灌前和生育期不同閾值條件下土壤表層含鹽量、地下水埋深和礦化度的概率分布，從概率空間分布的角度提出了不同時期防治土壤鹽漬化的地下水臨界埋深和礦化度［18］。相對其他土壤鹽漬化發育地區，青銅峽灌區近年來土壤鹽漬化的研究主要集中于利用多源遙感、實測高光譜數據開展鹽分反演、評價及動態監測［19］，草甸濕地等特殊土地利用類型鹽分特征分析［20］，灌區鹽漬土空間分布特征［21］等方面，針對整個灌區尺度的鹽漬化風險評價研究較少。本文以青銅峽灌區為研究對象，運用指示克里格法分析不同深度土壤鹽分、潛水位埋深和礦化度的空間概率分布規律，分析土壤鹽漬化與潛水的空間分布關系，以期為地下水資源調控、土地資源高效利用、鹽漬土防治等提供參考。

1 研究區概況

青銅峽灌區南起青銅峽水利樞紐，北至石嘴山，西抵賀蘭山，東鄰鄂爾多斯臺地西緣，可分為銀北、銀川和銀南3 個片區（圖1）。灌區多年平均降水量178 mm，年蒸發量接近年降水量的10 倍，日照時數2 870～3 080 h。灌區分布的土壤以灌於土、灰鈣土為主。糧食作物以水稻、玉米、小麥為主，經濟作物以蔬菜和瓜果類為主。青銅峽灌區主要由黃河沖積平原構成，海拔1 100～1 200 m，地勢由南向北緩緩傾斜，坡降約1/4 000。灌區地下水類型為第四系松散巖類孔隙水，分為單一潛水區和多層結構區，前者位于灌區最南端，系黃河出青銅峽峽口所形成的沖積扇，從西南向東北，含水層巖性由卵礫石逐漸變為含礫粉細砂，含水層厚度變厚；多層結構區分布在單一潛水區以外地區，在大約250 m 深度以上的范圍內，從上向下劃分為潛水、第一承壓、第二承壓含水巖組，各含水巖組間通常具有相對較連續的弱透水層，其中潛水含水層厚度一般為20～60 m，含水巖組巖性表現出自灌區邊緣向中心及由南向北逐漸變細、淤泥質含量增多的特點。引黃灌溉渠系滲漏補給和田間入滲補給占灌區地下水補給量的80%左右，潛水蒸發占地下水排泄量的47%左右［22］。1998-2017 年青銅峽灌區潛水位埋深下降了0.69 m，增加速率為0.038 m/a，年內地下水埋深呈雙峰雙谷特征［23］。2019 年潛水溶解性總固體為0.21～21.95 g/L。土壤鹽漬化一直是制約青銅峽灌區農業發展的重要因素，而鹽漬化的形成演化過程與平緩的地形、較淺的潛水位埋深、強烈的垂向蒸發排泄和局部地下水的高礦化度，以及引黃灌溉導致的鹽分積累等因素密切相關。

圖1 土壤采樣點和潛水監測點分布圖Fig.1 Distribution of soil sampling points and phreatic water monitoring points

2 數據與方法

2.1 數據來源

本研究使用了土壤可溶鹽、潛水位埋深、潛水礦化度3類數據（圖1）。其中:土壤可溶鹽數據來自于2018 年4 月和2019 年4 月春灌前開展的2 次野外采集樣品的測試，通過綜合考慮灌區農業種植結構、土地利用類型和地下水狀況等條件，盡可能均勻地布設了80 個采樣點，分層采集了0～20、20～40、40～60 cm深度土壤可溶鹽樣品；各采樣點使用GPS 定位，并對當地土地耕作、灌溉制度、農作物發育特征等信息進行調查；采集的土壤樣品由寧夏回族自治區地質礦產中心實驗室按照《土工試驗方法標準（GB/T50123-1999）》測試，其中K++ Na+使用火焰光度法測定、SO2-4用比濁法測定、Ca2+和Mg2+使用EDTA 滴定法測定、HCO-3和CO2-3用雙指示劑滴定法測試、Cl-用AgNO3滴定法測試、pH用pH計測定，采用離子加和法計算土壤全鹽量。潛水位埋深數據來自于灌區內117眼國家級和省級監測井自動化監測數據，另外選用了218 個同時期潛水位調查點數據作為補充，數據時段均為2019年4月25-30日。潛水礦化度數據來自于灌區內103 眼國家級監測井水質測試報告，數據時段為2019 年4月。國家級監測點潛水位埋深和水質數據按照國家地下水監測工程要求的標準開展全流程質量控制［24］，補充的潛水位調查點采用水位測量儀測量地下水埋深，使用RTK 獲取點的位置和高程信息。

2.2 研究方法

采用軟件GS+9.0 確定不同閾值下土壤全鹽量、潛水位埋深、潛水礦化度的變異函數模型，并將其模型參數輸入Arc-GIS10.6 進行指示克里格插值，進而繪制概率分布圖。步驟是:①根據研究需要確定土壤全鹽量、潛水位埋深、潛水礦化度評價的閾值。②根據設定閾值利用指示函數對土壤全鹽量、潛水位埋深和潛水礦化度數據進行二態指示變換，將得到的指示變換值（即1 或0）作為評價數據。③在GS+9.0 對獲取的指示變換值進行空間變異分析，獲取最佳變異函數模型及其參數。④將變異函數模型參數輸入ArcGIS10.6 進行指示克里格插值，得到各指標滿足相應閾值的空間概率分布圖。經典統計及相應的圖件制作采用SPSS25.0、GraphPad Prism 9.3.1軟件進行。

3 結果與分析

3.1 統計特征值分析

表1 為統計特征值分析結果。0～20、20～40、40～60 cm 深度土壤全鹽量平均值分別為12.27、4.85、4.75 g/kg，潛水位埋深、潛水礦化度的平均值分別為4.61 m 和2.55 g/L。變異系數Cv反映樣點的離散程度，對于估計結果可以提供一些預警信息，Cv＜0.1為弱變異性、0.1≤Cv≤1為中等變異性，Cv＞1為強變異性并標識存在一些特別大的樣本值［25］。青銅峽灌區不同深度土壤全鹽量、潛水位埋深、潛水礦化度的Cv均大于1，都屬于強變異，造成這種現象的原因首先與研究區范圍較大而產生的地形地貌、地下水賦存條件、土壤及地層特征等結構性因素有關，另外也與灌排方式、種植結構、耕種方式等隨機性因素的作用有關。單樣本K-S 正態檢驗結果表明，均不服從正態分布或對數正態分布。從圖2 可知，0～20、20～40、40～60 cm 土壤全鹽量頻率分布分別集中在小于上四分位數的10.35、4.69、4.98 g/kg 區間，潛水位埋深和潛水礦化度頻率分布也集中在小于上四分位數的4.85 m 和2.44 g/L 的區間，由于都在局部地區存在高值，頻率分布呈現明顯的“高頂”現象。

表1 青銅峽灌區土壤全鹽量、潛水位埋深和潛水礦化度統計特征值Tab.1 Descriptive statistics of total soil salt、groundwater depth and TDS in Qingtongxia Irrigation Area

圖2 青銅峽灌區土壤全鹽量、潛水位埋深和潛水礦化度小提琴圖Fig.2 Violin plot of total soil salt， groundwater depth and TDS in Qingtongxia Irrigation Area

3.2 指示變異函數模型構建

3.2.1 閾值的選擇

綜合考慮青銅峽灌區地形地貌、水文地質條件以及土壤鹽漬化程度等多方面因素確定各變量閾值。其中根據中國鹽漬土劃分標準［2］和相關研究成果［26］，選取土壤全鹽量1.0、2.0、4.0、10.0 g/kg 分別作為土壤輕度、中度、重度鹽化及鹽土化的閾值；參考Zhang 等人的研究［27］，選取潛水礦化度為1.0、2.0、2.5、3.0 g/L 為水質咸化閾值；參考金曉媚［28］、鄧麗［29］、周在明［16］等人研究，選取1.5、2.0、2.5、3.0 m作為潛水蒸發作用的臨界水位埋深。二態指示變換規則是，土壤全鹽量和潛水礦化度大于相應閾值為1、小于相應閾值為0；潛水位埋深為大于相應閾值為0、小于相應閾值為1。

3.2.2 指示變異函數模型

指示克立格法通過指定閾值對數據進行二態指示變換，可以削弱偏態分布和特異值對變異函數穩健性產生的影響。表2是各變量不同閾值指示變換后的指示變異函數模型和相關參數?？梢钥闯?青銅峽灌區不同深度、不同閾值條件下土壤全鹽量的指示半方差函數符合球狀模型、指數模型、高斯模型等不同類型。潛水位埋深的指示半方差函數符合指數模型。潛水礦化度閾值為1.0 和2.0 g/L 時指示半方差函數符合指數模型，閾值為2.5 和3.0 g/L 時指示半方差函數符合高斯模型。從決定系數和殘差看，理論模型均顯示出較好的擬合程度。

表2 不同閾值下土壤全鹽量、潛水位埋深和潛水礦化度的指示變異函數理論模型Tab.2 Indicator semivariogram models of soil salinity， phreatic water depth and phreatic TDS for different threshold values

3.2.3 不同變量的空間異質性分析

塊金值C0由測量誤差和空間變異共同作用引起，反映區域化變量的隨機異質性程度；基臺值Sill表示結構性和隨機性變異構成的總變異，表示變量空間異質性的強弱，C0/Sill為＜25%、25%～75%、＞75%對應空間相關性較強、中等和很弱，在整個尺度上具有恒定變異時其值將接近100%。變程是不同變量空間異質性的尺度函數，用來衡量最大變異程度的空間距離，在變程之內具有空間自相關性［30］。區域化變量的空間變異是由結構性、隨機性因素共同引起的，其中氣候、地形地貌、土壤母質、地質與水文地質條件等非人為的結構性因素是變量具有空間連續性的原因，而灌溉制度、種植結構、耕作措施、地下水開發等隨機因素則會弱化其空間自相關性。根據表2，各模型塊金值C0都較小，說明在本研究尺度上由于采樣誤差、短距離的變異、隨機和固有變異引起土壤全鹽量、潛水位埋深和潛水礦化度的變異都不大。土壤全鹽量C0/Sill的不同反映出在不同深度由空間自相關引起的空間異質性所占的比例有差異，總體表現出隨土層深度增加呈倒U 型特點:0～20 cm 閾值為1.0 g/kg，20～40 cm 閾值為1.0 g/kg 和4.0 g/kg 和40～60 cm 閾值為4.0 g/kg 時，C0/Sill均小于25%，表現出較強的空間自相關，即空間變異主要由結構型因素作用；而20～40 和40～60 cm 閾值為10.0 g/kg 時，C0/Sill為100%，說明在整個研究尺度具有恒定變異；其余條件下C0/Sil處于25%～75%之間，表現出中等的空間自相關，也即空間變異由結構型因素和隨機性因素共同作用。潛水位埋深閾值為最小的1.5 m 和最大的3.0 m 時，C0/Sill分別為5.82%和12.59%，表現出較強的空間自相關，空間變異主要由水文地質條件等結構型因素作用，而閾值為2.0 和2.5 m 時，C0/Sill處于25%～75%之間，表現出中等的空間自相關，說明人類活動等隨機性因素對空間變異的影響程度增加。研究區潛水礦化度的C0/Sill在4種閾值條件下均＜25%，表現出較強的空間自相關，即其空間變異主要由結構性因素控制的，具體表現在由于含水層巖性、透水能力、水力坡度、水循環能力等的不同，造成灌區單一潛水區和多層結構區潛水位埋深和TDS 特征的明顯區別，即單一淺水區潛水埋藏較深而TDS 多小于0.5 g/L；而在多層結構區，潛水的水位埋深和水化學成分在空間上的變化都較大，永寧縣以南地下水水力坡度大、徑流條件好、水交替作用強烈，大部分為低礦化淡水，而在永寧縣以北隨著地下水徑流條件的逐漸變差和水位埋深的逐漸變淺，地下水質逐漸變差。另外，從表2 知，潛水位埋深和潛水礦化度的變程分別在8.70～45.30 和12.99～19.20 km，自相關范圍均不大（約為青銅峽灌區最大直線范圍的1/10 到1/3）。不同深度土壤全鹽量的變程總體隨著閾值的變大而增大，比如0～20 cm 閾值為1.0 和2.0 g/kg 時自相關范圍較小，而閾值為4.0 和10.0 g/kg 時，自相關距離已超過青銅峽灌區最大直線范圍的1/2，反映出當閾值較小時表層土壤全鹽量受自然因素和人類活動的影響較大，而隨著閾值的增大土壤自身性狀對其鹽分特征的影響變大。

3.3 土壤全鹽量與潛水環境的概率空間分布特征

將變異函數模型參數輸人ArcGIS10.6 進行指示克里格插值，繪制各指標滿足相應閾值的概率空間分布圖，并定義概率大于40%為高概率區、反之為低概率區。圖3（a）～（d）、（e）～（h）和（i）～（l）分別是不同深度土壤處于輕度（＞1 g/kg）、中度（＞2 g/kg）、重度鹽化（＞4 g/kg）和鹽土化（＞10 g/kg）的概率分布圖，表3是不同閾值條件下各指標高低概率區面積統計情況。可以看出，春灌前青銅峽灌區0～20 cm 深度土壤發生輕度、中度、重度鹽化和鹽土化的高概率區占灌區總面積的比例為92.38%、73.71%、51.83%、24.61%；20～40 cm深度土壤發生輕度、中度、重度鹽化和鹽土化的高概率區占灌區總面積的比例為92.15%、50.21%、6.87%、0.53%，40～60 cm深度土壤發生輕度、中度、重度鹽化和鹽土化的高概率區占灌區總面積的比例為89.16%、43.06%、6.09%、0.46%；即隨著閾值增大，各深度高概率區明顯減少，這與灌區鹽化表現出的表聚型特點相符［21］。對比圖1 發現，當土壤深度為0～20 cm，閾值為2 g/kg 時，銀北片區的大部分、銀川片區的全部均為高概率區，其面積約占灌區總面積的73.71%；當閾值為4 g/kg 時，高概率區面積減少至3 174.25 km2，主要分布在銀川片區全部區域和銀北片區的局部；當閾值為10 g/kg 時，除銀川片區中部和銀北片區石嘴山南部的小范圍外約1 507.33 km2區域外，其余均演化為低概率區。當土壤深度為20～40和40～60 cm，發生中度鹽化（閾值為2 g/kg）的高概率區分別占灌區總面積的50.21%和43.06%，主要分布在銀北片區東北部、西南部和銀川片區的北部；發生重度鹽化（閾值為4 g/kg）的高概率區面積分別減少至420.74 和372.84 km2，主要集中在銀北片區的西南部；發生鹽土化（閾值為10 g/kg）的高概率區面積進一步減少至32.48 和28.45 km2，占比也降至0.50%左右。總的來看，不同深度大閾值的高概率區均包含在上一級小閾值的高概率區，并在空間分布上具有明顯的相似性，尤其是深部土壤重度鹽化和鹽土化高風險區主要集中銀北片區的西南部一帶。

圖3 青銅峽灌區不同閾值下不同深度土壤全鹽量概率空間分布Fig.3 Probability maps of soil salinity at different depths under different thresholds in Qingtongxia Irrigation Area

圖4（a）～（d）為不同閾值潛水位埋深的概率分布，對比圖1發現，不同閾值概率分布圖的共同特點是:銀北片區是高概率集中分布區，低概率區主要集中在銀川片區中部和銀南片區。當潛水位埋深閾值為1.5 m 時，研究區潛水埋深概率分布區間主要在0～0.2 的低概率之間，僅13.42%的區域為高概率區且主要分布在銀北片區的東北部和銀南片區東南部局部小區域；當閾值為2.0 m 和2.5 m 時時，高概率區面積在原有基礎上分別擴大至約1 852.48、2 564.33 km2，占比增加至30.25%和41.87%，增加部分均主要來源于銀北片區；隨著閾值的增加為3.0 m 時，高概率風險區繼續增加，僅約34.82%的區域為低概率區且主要集中分布于銀川片區中的城市建設區周圍和銀南片區的單一潛水區，其余地區基本為高概率區。分析不同閾值下潛水位埋深概率圖，發現小閾值的高概率區包含在大閾值的高概率區，其分布區域隨閾值的增大而擴散。

圖4 青銅峽灌區不同閾值下潛水位埋深概率空間分布Fig.4 Probability spatial distribution of phreatic water depth under different thresholds in Qingtongxia Irrigation Area

圖5（a）～（d）為不同閾值下潛水礦化度的概率分布，當閾值為1.0、2.0、2.5、3.0 g/L 時高概率區面積分別是3 621.55、3 381.99、2 067.43、1 957.72 km2，占比分別為59.13%、55.22%、33.76%和31.96%。結合圖1發現，不同閾值下銀川片區和銀北片區都是高概率集中分布區，低概率區主要集中于銀南片區，同時表現出大閾值的高概率區包含在小閾值的高概率區，其分布區域隨閾值的增大逐漸縮小的變化特點。

圖5 青銅峽灌區不同閾值下潛水礦化度概率空間分布Fig.5 Probability spatial distribution of phreatic TDS under different thresholds in Qingtongxia Irrigation Area

3.4 土壤全鹽量與潛水環境的概率空間分布關系

分析圖3、圖4 和圖5 發現3 個不同深度土壤發生輕度鹽化（閾值為1 g/kg）高風險區基本覆蓋整個灌區，這與由于本研究土壤可溶鹽樣品采集于4月份的春灌前，此時土壤處于融凍期，返鹽激烈且不能及時脫鹽導致鹽分積累加重的事實有關；另外3個深度發生鹽土化（閾值為10 g/kg）的范圍分布很小。即分析不同深度土壤發生輕度鹽化和鹽土化的概率空間分布，與潛水埋深和礦化度的關系意義不大。以下重點對比深度為20～40 cm 時土壤發生中度（閾值為2 g/kg）、重度鹽化（閾值為4 g/kg）的概率空間分布與潛水位埋深和礦化度之間的關系。使用ArcGIS10.6 的Intersect 功能計算不同閾值土壤全鹽量高概率區與潛水位埋深、礦化度高概率區相交的部分，即同時滿足土壤全鹽量相應閾值和潛水環境相應閾值高概率分布的區域，并統計匹配率發現（表4），青銅峽灌區潛水位埋深越淺，表層土壤含鹽量越高，反之愈輕；而潛水礦化度愈大，土壤鹽漬化愈嚴重，反之也愈輕。這與王金哲等在西北干旱區［31］和樊自立等在塔里木河流域［32］的研究結論相同。可能的原因是潛水埋深越淺，地下水中的鹽分越易通過潛水蒸發由地下水帶至土壤耕層而表聚；特別是當潛水礦化度較大時，越易發生次生鹽漬化。因此可以認為存在臨界埋深和礦化度，即當潛水位埋深小于某臨界值時、礦化度大于某臨界值時，土壤表層發生鹽漬化的風險就大。對比圖3 和圖4 發現，潛水埋深2.0 m 及以上閾值的概率空間分布與土壤全鹽量2.0 g/kg和4.0 g/kg的概率空間分布有較大相似性，且隨著潛水埋深的增大匹配率逐漸增大，并且明顯看出潛水埋深2.0 m 是匹配率突增的臨界值，可以初步判定發生中度和重度鹽化的潛水埋深臨界值為 2.0 m，與黃權中在河套灌區的研究結果相同［33］。再對比圖3 和圖5 可知，土壤全鹽量閾值為2.0 g/kg 的高概率分區范圍與潛水礦化度1.0 g/L 的相似性最大，匹配率達到74.22%；與潛水礦化度2.0 g/L 的相似性也較大，匹配率達到69.34%；隨著潛水礦化度閾值增加至2.5 g/L，匹配率大幅降低至39.02%；土壤全鹽量閾值為4.0 g/kg 的高概率分區范圍與潛水礦化度2.0 g/L 的相似性最大，匹配率達到59.45%；與潛水礦化度1.0、2.5 和3.0 g/L 匹配率都僅為30%左右；即可以初步判定發生中度和重度鹽化的潛水礦化度臨界閾值為2.0 g/L。以上分析表明青銅峽灌區土壤鹽分的空間變異性與潛水位埋深和礦化度的空間分布密切相關，但同時也發現，土壤全鹽量的高概率空間分布并不與潛水位埋深和礦化度的高概率空間分布完全相同，這是因為土壤鹽分特征還受到的氣象、地形地貌、土地利用方式等因素的共同影響［21］，而農業生產活動的頻繁程度、作物類型、施肥措施和灌溉水量及時間也是影響土壤鹽分的重要因素［25，34］。因此在青銅峽灌區影響農田土壤鹽漬化的一個重要因素是灌區潛水埋深和礦化度。因此，有必要進一步健全灌排系統，將大水漫灌用田間畦灌、滴灌、噴灌等節水型灌溉措施來替代，同時結合其他土壤改良方法，多渠道聯合控制土壤鹽漬化。

表4 20～40cm土壤全鹽量與潛水環境匹配率統計表Tab.4 Statistical table of matching rate between total salt content of 20～40 cm soil and phreatic environment

4 結 論

根據對青銅峽灌區2018、2019兩年4月份春灌前0～20、20～40 和40～60 cm 深度土壤全鹽量、潛水水位埋深、潛水礦化度進行空間變異性分析表明。

（1）青銅峽灌區土壤全鹽量、潛水位埋深、潛水礦化度的變異系數均大于1，屬于強變異；頻率分布呈現明顯的“高頂”現象。

（2）青銅峽灌區不同深度不同閾值土壤全鹽量的指示半方差函數符合球狀模型、指數模型、高斯模型等不同類型。潛水位埋深的指示半方差函數符合指數模型。潛水礦化度閾值為1.0 和2.0 g/L 時指示半方差函數符合指數模型，閾值為2.5 和3.0 g/L 時指示半方差函數符合高斯模型。不同閾值條件下的土壤全鹽量和潛水位埋深呈中等—較強的空間自相關，潛水礦化度表現出較強的空間自相關。

（3）青銅峽灌區0～20 cm 深度土壤發生輕度、中度、重度鹽化和鹽土化的高概率區占灌區的92.38%、73.71%、51.83%、24.61%；20～40 cm 深度比例為92.15%、50.21%、6.87%、0.53%；40～60 cm 深度比例為89.16%、43.06%、6.09%、0.46%。潛水位閾 值 為1.5、2.0、2.5、3.0 m 時 高 概 率 區 占 灌 區 的13.42%、30.25%、41.87%、65.18%。潛水礦化度閾值為1.0、2.0、2.5、3.0 g/L 時高概率區面積占灌區的59.13%、55.22%、33.76% 和31.96%。隨著閾值增加，高風險區面積均明顯減少。

（4）青銅峽灌區土壤發生中度和重度鹽漬化時的潛水臨界埋深為2.0 m，潛水臨界礦化度為2.0 g/L。