翟文華，付 宇，曹文庚，李澤巖，任 宇

（1.華北水利水電大學 地球科學與工程學院，河南 鄭州 450046；2.華北水利水電大學 測繪與地理信息學院，河南 鄭州 450046；3.中國地質(zhì)科學院 水文地質(zhì)環(huán)境地質(zhì)研究所，河北 石家莊 050061）

砷（As）是地殼的一種天然元素，無機砷化合物是一種毒性很強的物質(zhì)。部分地區(qū)飲用水和灌溉地下水中存在砷超標現(xiàn)象。當?shù)叵滤猩橘|(zhì)量濃度超過世界衛(wèi)生組織（WHO）設(shè)定的生活飲用水最高容許質(zhì)量濃度10 μg／L時，便可認為是高砷地下水。長期飲用高砷水會損害人體胃腸道系統(tǒng)、皮膚和神經(jīng)系統(tǒng)［1］。我國是受砷影響較為嚴重的國家之一，塔里木盆地、額濟納盆地、黑河盆地、柴達木盆地、東北平原、華北平原等地區(qū)被確定為潛在高砷地區(qū)［2］。黃河下游沖積平原豫北區(qū)位于華北平原南部，淺層地下水水質(zhì)較差，砷超標問題較為嚴重［3］，如新鄉(xiāng)市封丘縣曹崗鄉(xiāng)高砷暴露區(qū)病人癥狀表現(xiàn)為皮膚色素沉著或色素脫失［4］，這表明高砷地下水已經(jīng)對當?shù)赜盟踩腿梭w健康造成威脅。江欣悅等［5－6］研究黃河下游沖積平原豫北區(qū)高砷地下水空間分布及成因表明，該地區(qū)地下水還原性較強，不同沉積環(huán)境條件下生物地球化學和水文地質(zhì)作用是砷富集的驅(qū)動因素。

高砷地下水空間分布研究可以幫助人們快速識別高砷地下水影響范圍。隨著科學技術(shù)的發(fā)展，機器學習為非線性問題的解決提供了較好的方法，并且在地下水砷研究方面取得了豐富的成果。Bindal等［7］采用1 473個地下水砷采樣數(shù)據(jù)，結(jié)合地球化學環(huán)境、含水層條件和地形等20個變量，基于隨機森林和單變量特征選擇的混合集成模型預測了印度北方邦高砷地下水的空間分布。Podgorski等［8］采用23 799個地下水砷實測值和26個環(huán)境參數(shù)，利用隨機森林模型預測了印度高砷地下水的空間分布，研究表明印度有18萬～3 000萬人處于地下水砷質(zhì)量濃度超過10 μg／L的地區(qū)。Liang等［9］采用反向傳播神經(jīng)網(wǎng)絡(luò)方法對中國臺灣蘭陽平原地下水砷濃度進行了預測，與普通克里金法比較發(fā)現(xiàn)，反向傳播神經(jīng)網(wǎng)絡(luò)方法具有更高的預測精確度。當前研究大多采用單一機器學習模型預測高砷地下水的空間分布，其中線性分類器處理非線性數(shù)據(jù)時性能不佳，神經(jīng)網(wǎng)絡(luò)方法則容易陷入局部極小值，且學習收斂速度較慢，因此傳統(tǒng)的機器學習模型預測高砷地下水空間分布難度較大。為了解決該問題，本文采用堆疊（Stacking）集成學習模型，充分發(fā)揮基學習器中每個機器學習算法的優(yōu)勢，降低單一算法泛化能力不佳的風險，提高模型的預測精確度。

為了識別黃河下游豫北區(qū)地下水中潛在高砷區(qū)域，基于研究區(qū)1 081個地下水砷質(zhì)量濃度實測值，選取人類活動、氣候、沉積環(huán)境、土壤理化特征、水文地質(zhì)等作為環(huán)境因子，采用Stacking集成學習模型對黃河下游豫北區(qū)高砷地下水空間分布進行建模，預測該地區(qū)高砷地下水的空間分布，明確在預測過程中環(huán)境變量對地下水砷分布的相對重要性，以期為該地區(qū)地下水資源的合理利用和有效管理提供依據(jù)。

1 研究方法

1.1 研究區(qū)概況

研究區(qū)位于黃河下游河南省北部沖積平原地區(qū)（見圖1，其中ρ（AS）為砷質(zhì)量濃度），面積19 733.75 km2，地勢總體由西南向東北傾斜。研究區(qū)屬于溫帶大陸性季風氣候區(qū)，四季分明，年均氣溫13.3～15.6℃，年降水量496.7～751.3 mm，降水集中在7—9月，年蒸發(fā)能力988.0～1 023.9 mm，5月、6月蒸發(fā)最為強烈。黃河和衛(wèi)河是研究區(qū)兩大主要河流。研究區(qū)為歷史上黃河決口、改道最頻繁的地區(qū)之一，地表仍有黃河河道變遷的遺跡。研究區(qū)地下水總體徑流方向與地勢變化基本一致，在太行山前由西向東流動，最終在華北平原東部入海。研究區(qū)地下水補給來源主要為地表水、大氣降水及灌溉水下滲，排泄以人工開采和蒸發(fā)為主［4］。

1.2 環(huán)境變量選擇及處理

含水層中砷的遷移和流動受地貌、水文地質(zhì)、生物地球化學和人類活動影響［10－16］，例如地形、地貌、沉積物特征、土壤性質(zhì)、土地利用類型、地下水流量和植被覆蓋等，有學者利用這些環(huán)境變量進行了地下水砷分布預測。Tan等［17］基于地面標高、坡度、氣溫、降水量、蒸散量、地貌類型、地下水水位等90個環(huán)境因子，利用提升回歸樹模型對孟加拉國高砷地下水分布進行了預測。綜合分析可能影響研究區(qū)高砷地下水空間分布的環(huán)境變量，選擇氣候、人類活動、沉積環(huán)境、水文地質(zhì)、土壤理化特征等6類共21個初始環(huán)境變量進行研究區(qū)高砷地下水分布預測，見表1（與河流距離為構(gòu)建距離河流間隔為500 m的多環(huán)緩沖區(qū)，地下水水位累計變幅計算年份為1949—2020年，水位年際變化計算年份為2019—2020年，黃河決口情況為歷史累計黃河決口次數(shù)，第四紀地貌類型有沖洪積扇及沖洪積平原、泛濫平原、海積平原、湖泊、洼地、河道帶、黃土狀土、基巖，淺層和深層土壤理化特征包括砂粒含量、粉砂含量、黏土含量、土壤有機碳含量、土壤pH值，土地利用類型包括耕地、建筑用地、林地、水系、草地）。研究區(qū)環(huán)境變量均按照500 m空間分辨率劃分像元。由于環(huán)境變量較多，其中可能含有與研究區(qū)高砷地下水空間分布預測不符的冗余信息和噪音，因此利用以隨機森林為迭代分類器的遞歸特征消除技術(shù)消除冗余信息和噪音［18］。通過反復計算各環(huán)境變量對砷分布重要性影響占比，并刪除重要性影響占比小的環(huán)境變量，最終選出含有18個環(huán)境變量的最佳變量子集，即刪除黏性土層、富水性、土地利用類型3個環(huán)境變量。

表1 模型預測變量及對砷分布重要性影響占比

1.3 模型構(gòu)建與驗證

在利用機器學習算法進行地下水水質(zhì)建模時，線性模型（邏輯回歸、支持向量機SVM和線性判別分析LDA等）、極端梯度提升（XGBoost）、隨機森林（RF）、神經(jīng)網(wǎng)絡(luò)等都有較好的預測性能［19－22］，其中隨機森林模型可以較好處理高維數(shù)據(jù)、異常值、噪聲、過擬合和多重共線性問題；極端梯度提升模型可以很好地處理高維數(shù)據(jù)，不易陷入過擬合，同時通過多線程并行計算提高計算效率；支持向量機可以將非線性問題轉(zhuǎn)化為某個高維空間的線性問題；線性判別分析模型簡單，無需調(diào)參，且可以較好地處理不同類別訓練樣本量差別很大的數(shù)據(jù)。相對于以上獨立機器學習模型，Stacking模型的優(yōu)勢在于將不同類型的機器學習模型進行組合。不同類型機器學習模型對于數(shù)據(jù)空間以及結(jié)構(gòu)的學習存在較大差別，Stacking模型可以從不同角度觀察數(shù)據(jù)特征，更加全面地學習數(shù)據(jù)，從而得到更加準確的結(jié)果。目前Stacking模型在醫(yī)學、大氣污染、金融和工程造價等領(lǐng)域廣泛應用。孫麟［23］通過Stacking模型將多種算法結(jié)合起來，提高了病人是否患有宮頸癌的預測精確度，為醫(yī)學相似問題的預測提供了新思路；趙濱等［24］利用Stacking模型建立地面PM2.5濃度估算模型，提高了大范圍區(qū)域大氣污染預測的精確度。

基于以上論述，根據(jù)Stacking模型“基學習器好而不 同，元 學 習 器 好 而 結(jié) 構(gòu) 簡 單”的 原 則［25］，將XGBoost、RF、SVM作為Stacking模型的基學習器，LDA作為Stacking模型的元學習器進行融合得到集成模型，用于預測高砷地下水的空間分布。Stacking模型構(gòu)建流程見圖2。

圖2 Stacking模型構(gòu)建流程

二元目標變量的預測可以提高模型預測的準確性和有效性，因此將10 μg／L作為閾值，將砷質(zhì)量濃度≤10 μg／L、＞10 μg／L重新編碼為0、1。砷質(zhì)量濃度≤10 μg／L、＞10 μg／L的采樣數(shù)據(jù)在數(shù)據(jù)量上存在差異，需要進行類平衡處理，將經(jīng)過重新編碼和類平衡處理的362個采樣點數(shù)據(jù)按7∶3劃分訓練集（X＿－train，Y＿train）和測試集（X＿test，Y＿test），其中X＿train、X＿test為篩選的18個環(huán)境變量數(shù)據(jù)集，Y＿train、Y＿test為砷質(zhì)量濃度編碼數(shù)據(jù)集。首先，利用訓練集采用RF、XGBoost、SVM模型進行五折交叉驗證，訓練后每種模型得到一組與訓練集數(shù)量相同的數(shù)據(jù)，將這3組數(shù)據(jù)進行組合，作為元學習器的特征數(shù)據(jù)。在對訓練集進行5折交叉驗證的同時，對測試集進行計算，不同的是需要將計算結(jié)果進行平均，將3種模型計算的結(jié)果進行組合，作為元學習器的測試集。然后，將得到的元學習器的特征數(shù)據(jù)用于LDA模型訓練，測試集用于驗證Stacking模型的性能。通過ROC曲線下面積（AUC）、準確率（Accuracy）、特異性（Specificity）和敏感性（Recall）對模型性能進行評估，其中：AUC值越大，模型性能越好；Accuracy是全部預測正確的樣本占所有樣本的比例；Recall指正確分類砷質(zhì)量濃度＞10 μg／L樣品的能力；Specificity指正確分類砷質(zhì)量濃度為≤10 μg／L樣品的能力。最后，利用構(gòu)建好的Stacking模型，對黃河下游沖積平原豫北區(qū)高砷地下水空間分布進行預測，并繪制研究區(qū)高砷地下水概率分布圖。

2 結(jié)果與討論

2.1 地下水砷質(zhì)量濃度統(tǒng)計與分布特征

研究區(qū)地下水采樣數(shù)據(jù)來自2010—2020年中國地質(zhì)科學院水文地質(zhì)環(huán)境地質(zhì)研究所與河南省自然資源監(jiān)測院進行的水文地質(zhì)調(diào)查結(jié)果，調(diào)查采集淺層地下水樣品共1 081組，采樣點位置見圖1。經(jīng)計算，研究區(qū)地下水砷質(zhì)量濃度最大值、最小值、平均值分別為0.1、190、7.06 μg／L，中值、標準差分別為1.20、16.52 μg／L，變異系數(shù)為2.34，具有強變異性，表明研究區(qū)地下水砷質(zhì)量濃度具有較強的空間變異特征。研究區(qū)高砷地下水采樣點共181個，超標率為16.74%。從采樣點分布看，高砷地下水分布呈東北—西南走向，研究區(qū)中部、南部高砷地下水主要分布在太行山前洼地及黃河決口扇地區(qū)；東部以濮陽為中心，高砷地下水呈不均勻分布。

2.2 地下水高砷分布特征

XGBoost、RF、SVM、Stacking模型在測試集上的評估結(jié)果見表2，Stacking模型具有最大的AUC、Accuracy、Specificity和Recall值。從模型評估指標來看，Stacking模型預測精確性最高，Recall、Specificity分別為0.75、0.877 2，表明Stacking模型可以準確預測研究區(qū)砷質(zhì)量濃度≤10 μg／L、＞10 μg／L的區(qū)域。

表2 不同模型性能對比 %

采用XGBoost、RF、SVM、Stacking模型計算研究區(qū)各單元地下水中砷質(zhì)量濃度超過10 μg／L的概率，計算結(jié)果見圖3。XGBoost、RF、SVM、Stacking模型計算的高砷地下水空間分布總體趨勢相似，研究區(qū)地下水砷污染范圍較大，主要集中在中部、南部和東部部分地區(qū)，但不同模型計算的局部區(qū)域差異明顯。XGBoost、RF模型計算結(jié)果中，新鄉(xiāng)市封丘縣、安陽市滑縣北部高砷地下水分布概率較大，不能精確展現(xiàn)高砷地下水分布情況。SVM模型計算結(jié)果不能細致描繪局部地區(qū)高砷地下水分布變化情況。Stacking模型計算結(jié)果可以較好表達局部高污染區(qū)域細部變化，局部沒有采樣點的地區(qū)也能較好地進行預測，且與實際情況相符，原因是，Stacking模型能夠充分結(jié)合不同機器學習算法的優(yōu)勢來消除誤差。Stacking模型計算的高砷概率為0.09～0.88，根據(jù)概率閾值0.5劃定高砷地區(qū)面積為6 673.25 km2，占研究區(qū)總面積的33.82%，高砷地下水集中分布在研究區(qū)太行山前洼地及黃河決口扇地區(qū)，具體分布在新鄉(xiāng)市延津縣、原陽縣和封丘縣北部，安陽市滑縣和內(nèi)黃縣南部，濮陽市范縣和濮陽縣等地。其中高砷地下水分布概率大于0.8的地區(qū)主要分布在新鄉(xiāng)市延津縣、原陽縣和封丘縣北部，面積為1 237.25 km2，占研究區(qū)總面積的6.27%，這些地區(qū)應加強地下水水質(zhì)監(jiān)測與管理。

圖3 地下水砷質(zhì)量濃度超過10 μg／L的概率分布

2.3 驅(qū)動因素分析

特征重要性分析可以用來評估不同環(huán)境變量對高砷地下水分布的影響。基于隨機森林的gini方法計算得到研究區(qū)不同環(huán)境變量對砷分布重要性影響占比，見表1。黃河決口情況、年均氣溫、年降水量、高程、水力梯度是對高砷地下水分布影響排序靠前的環(huán)境變量，根據(jù)環(huán)境變量重要性綜合排序結(jié)果可以得出，沉積環(huán)境（黏砂比、黃河決口情況、第四紀地貌類型、高程）對地下水砷富集具有較大影響。黃河決口情況對地下水砷空間分布影響最大，原因是黃河決口和改道形成了砂與土互層的沉積環(huán)境。研究區(qū)高砷地下水主要分布在太行山前洼地與黃河決口扇所在的沉積環(huán)境中，其中太行山前洼地位于新鄉(xiāng)市北部與安陽市滑縣交界處，屬于沖洪積扇前洼地，地下水在此徑流不暢，含水層中砂層顆粒變細且厚度逐漸變薄，黏土厚度逐漸增大，水量變小，形成還原環(huán)境為主的地下水環(huán)境，隨著地下水還原性增強［26］，地下水中砷質(zhì)量濃度逐漸升高；在黃河沖積平原，黃河頻繁決口泛濫為該地區(qū)提供了大量的有機質(zhì)及較好的低氧－缺氧條件，加之受黃河側(cè)向補給影響，濃縮作用及還原環(huán)境導致該地區(qū)地下水砷質(zhì)量濃度較高［27－29］，即黃河頻繁決口地區(qū)有利于含水層中砷的釋放。

年均氣溫和年降水量對地下水砷分布也有較大影響，年降水量和年均氣溫通過調(diào)節(jié)地表徑流進而影響地下水砷質(zhì)量濃度，當年降水量增大時，河流湖泊等地表水水位上升，地表水補給地下水，對地下水中砷起稀釋作用。同時地表水中含有溶解氧，將氧氣等氧化劑帶入地下水，阻礙了砷的還原性釋放，地下水砷質(zhì)量濃度不斷降低。氣溫升高導致蒸發(fā)量增大，可利用地表水不斷減少，地下水補給量相應減少，對地下水的稀釋作用減弱，而研究區(qū)所處還原環(huán)境有利于砷的釋放和運移，因此地下水砷質(zhì)量濃度上升，這與曹永生等［30］的研究結(jié)果一致。

高程和水力梯度對高砷地下水的分布也有一定影響，高程和水力梯度主要通過地下水流速影響地下水砷質(zhì)量濃度。在低海拔和低水力梯度的平原地區(qū)，沉積物顆粒較細，地下水流速緩慢，增加水巖作用時間，O2、NO－3等氧化劑缺乏使得沉積物中氧化物被還原，吸附在沉積物上面的砷被釋放，水體中砷質(zhì)量濃度相應提高。在高海拔、高水力梯度的山前補給區(qū)，沉積物顆粒較粗，地下水流速快，地下水補給也快，將O2等氧化劑帶入含水層，不利于砷的富集。

3 結(jié)論與建議

通過構(gòu)建Stacking模型，預測了黃河下游沖積平原豫北區(qū)高砷地下水空間分布，分析了環(huán)境變量對高砷地下水分布的影響。研究區(qū)地下水砷質(zhì)量濃度為0.01～190 μg／L，超標率（砷質(zhì)量濃度＞10 μg／L）為16.74%，表明研究區(qū)地下水存在威脅人體健康的風險。高砷地下水主要分布在研究區(qū)太行山?jīng)_洪積扇前洼地及黃河決口扇地區(qū)，具體分布在新鄉(xiāng)市延津縣、原陽縣和封丘縣北部，安陽市滑縣和內(nèi)黃縣南部，濮陽市范縣和濮陽縣等地。構(gòu)建的Stacking模型相較于XG?Boost、RF、SVM獨立機器學習模型，具有最大的AUC、Accuracy、Specificity和Recall，Stacking模型性能優(yōu)于獨立機器學習算法模型的性能，可以提高當前機器學習模型對地下水砷分布預測的精確度。黃河決口情況、年均氣溫、年降水量、高程、水力梯度是影響高砷地下水分布的重要環(huán)境變量，沉積環(huán)境（黏砂比、黃河決口情況、第四紀地貌類型、高程）與地下水中砷富集顯著相關(guān)。

基于本文研究結(jié)果，針對黃河下游豫北區(qū)高砷地下水提出以下防治建議：地下水和土壤是相互作用的，在地下水受砷污染嚴重的地區(qū)，通過注射井向該地區(qū)土壤中的黏土層注入改性物質(zhì)及表面活性劑，經(jīng)改性后的黏土可以吸附地下水中的砷污染物；在條件允許的地區(qū)構(gòu)建集雨系統(tǒng)，利用雨季降水補給地下水，進而稀釋地下水中砷污染物；集中開發(fā)無法飲用的高砷水作為工業(yè)用水；在高砷地下水區(qū)域，打井抽取地下水造就地下水漏斗，使地下水漏斗區(qū)的水位明顯低于周圍地區(qū)，加快漏斗區(qū)地下水的循環(huán)和更新速度，達到改善水質(zhì)的目的；在地下水砷污染較為嚴重的地區(qū)，加快地下水動態(tài)監(jiān)測網(wǎng)絡(luò)建設(shè)和污染防治技術(shù)攻關(guān)，對地下水水質(zhì)進行監(jiān)測，以便提前采取地下水污染防治措施。