黃河下游豫北區(qū)高砷地下水空間分布研究

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

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

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

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

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

1 研究方法

1.1 研究區(qū)概況

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

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

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

表1 模型預(yù)測(cè)變量及對(duì)砷分布重要性影響占比

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

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

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

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

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

2 結(jié)果與討論

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

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

2.2 地下水高砷分布特征

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

表2 不同模型性能對(duì)比 %

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

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

2.3 驅(qū)動(dòng)因素分析

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

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

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

3 結(jié)論與建議

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

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