流域土壤重金屬分布式運(yùn)移模型研究

王圣偉 冀 豪 王 苗 婁天瀧 張 暢

(西北師范大學(xué)計(jì)算機(jī)科學(xué)與工程學(xué)院， 蘭州 730070)

0 引言

我國(guó)各流域尤其是中小流域土壤都面臨著程度不一的水資源環(huán)境污染問(wèn)題，對(duì)農(nóng)業(yè)及經(jīng)濟(jì)發(fā)展造成了深遠(yuǎn)的影響。借助GIS技術(shù)合理估算流域土壤重金屬運(yùn)移量，能夠?qū)r(nóng)業(yè)生產(chǎn)生態(tài)環(huán)境實(shí)現(xiàn)總體目標(biāo)控制。國(guó)內(nèi)外學(xué)者對(duì)于重金屬污染運(yùn)移的預(yù)測(cè)模型研究較為關(guān)注。RODE等[1]針對(duì)水環(huán)境的未來(lái)發(fā)展趨勢(shì)，提出了一種綜合水質(zhì)模型，該模型的模擬結(jié)果表明，在河流較寬的部分，重金屬積聚在通道邊緣附近，并且由于侵蝕作用可能在未來(lái)會(huì)重新遷移。Al-MUR等[2]分析了吉達(dá)附近紅海4個(gè)沉積物巖芯中重金屬的分布，結(jié)果表明，在紅海中重金屬M(fèi)n、Cu和Pb具有較高的含量，核心沉積物中的重金屬濃度在吉達(dá)中部附近增加，并且近年來(lái)已經(jīng)逐步升高。CHO等[3]基于Clark單位線圖及其空間分解方法的組合概念，并結(jié)合精細(xì)的空間變量流動(dòng)力學(xué)，實(shí)現(xiàn)了空間分布的水文模擬，開(kāi)發(fā)了一種集總概念和分布式特征混合水文模型。MEDIERO等[4]提出了一種基于分布式水文模型參數(shù)概率表示的標(biāo)定方法，該方法用于曼薩納雷斯河流域的案例研究，以校正后的模型作為實(shí)時(shí)洪水預(yù)報(bào)的決策支持工具。UNDUCHE等[5]開(kāi)發(fā)了4種水文模型(WATFLOOD、HBV-EC、HSPF和HEC-HMS)，從復(fù)雜的完全分布模型到相對(duì)簡(jiǎn)單的集總模型，并對(duì)其有效性和準(zhǔn)確性進(jìn)行了校準(zhǔn)和測(cè)試，以作為洪水預(yù)報(bào)的實(shí)用工具。這4種模型均能準(zhǔn)確地模擬流量，可用于洪水預(yù)報(bào)，但是WATFLOOD和HEC-HMS模型在模擬高流量時(shí)表現(xiàn)不佳，HSPF和HBV-EC模型在大多數(shù)情況下對(duì)低流量預(yù)測(cè)過(guò)高，且這4種模型僅應(yīng)用于草原流域的流量預(yù)測(cè)，對(duì)其他非草原流域的預(yù)測(cè)結(jié)果有待考證。

本研究以大夏河流域及其子流域的土壤和水體中重金屬為切入點(diǎn)，對(duì)重金屬的含量和遷移特征進(jìn)行統(tǒng)計(jì)和分析，根據(jù)大夏河流域重金屬的總體分布情況，選取5種易造成污染的重金屬As、Cd、Cr、Cu和Pb進(jìn)行分析，采用GIS圖像處理技術(shù)和數(shù)據(jù)分析方法對(duì)流域遙感和水文數(shù)據(jù)進(jìn)行處理，結(jié)合大夏河流域地貌形態(tài)特征和流域DEM數(shù)據(jù)進(jìn)行分形計(jì)算，提出流域土壤重金屬分布式運(yùn)移模型，對(duì)大夏河流域土壤中重金屬的運(yùn)移量進(jìn)行估算。

1 材料與方法

1.1 研究區(qū)概況

大夏河屬于黃河上游的一級(jí)支流，地處青藏高原與黃土高原過(guò)渡帶(東經(jīng)102°2′～103°23′，北緯34°51′～35°48′)。大夏河經(jīng)夏河縣城東北流入，出土門(mén)關(guān)后進(jìn)入臨夏盆地，自西向東橫貫臨夏市、至折橋轉(zhuǎn)北向東鄉(xiāng)縣泄湖峽，在康家灣喇嘛川塔張村注入劉家峽水庫(kù)[6]。主河道總長(zhǎng)203 km，總流域面積7 152 km2，多年平均流量34.3 m3/s，年均徑流量11億m3，多年含沙量3.49 kg/m3，年輸沙量4.192×1010kg，林地面積3.27×106hm2，占流域總面積的21.2%[7]。

1.2 流域分級(jí)與劃分

根據(jù)大夏河流域邊界的經(jīng)緯度，從地理空間數(shù)據(jù)云上高級(jí)檢索，下載包含大夏河流域的DEM圖像數(shù)據(jù)(共4塊)，將下載的原始DEM圖像導(dǎo)入ArcMAP軟件，需要應(yīng)用ArcMAP軟件中數(shù)據(jù)管理工具—柵格—柵格數(shù)據(jù)集—鑲嵌工具，把4塊分散的柵格數(shù)據(jù)圖合并到一起，然后根據(jù)大夏河的流域邊界把大夏河流域從4塊合并的柵格數(shù)據(jù)上裁剪，如圖1所示。

圖1 大夏河高程圖Fig.1 Digital elevation model of Daxia river

本文首先對(duì)大夏河流域選取合適的閾值[8-9]進(jìn)行河網(wǎng)水系提取，然后根據(jù)Strahler和Shreve這兩種流域分級(jí)法[10]，并結(jié)合ArcMAP軟件Spatial Analyst—水文分析—河網(wǎng)分級(jí)工具，對(duì)大夏河流域進(jìn)行分級(jí)處理。而流域的劃分首先需要確定各個(gè)子流域的出水口位置，然后根據(jù)水流流向[11-12]判斷流入該出水口的所有上游區(qū)域，即為子流域，在ArcGIS軟件中把大夏河流域劃分成7個(gè)子流域，如圖2所示。

圖2 大夏河子流域區(qū)域劃分Fig.2 Regional division of sub-basins of Daxia river basin

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

2016—2018年，使用手持式GPS記錄位置，采集大夏河流域周邊土樣和水樣，以全年存在穩(wěn)定徑流為依據(jù)，共26個(gè)采樣點(diǎn)，采樣位置具有代表性，見(jiàn)圖3。在野外將樣品放入自封袋中，編號(hào)后帶回實(shí)驗(yàn)室，對(duì)土壤樣品進(jìn)行四分法取樣后風(fēng)干，研磨過(guò)200目尼龍篩，放入對(duì)應(yīng)編號(hào)的樣品容器內(nèi)，對(duì)流域水樣搖勻，也裝入對(duì)應(yīng)編號(hào)的樣品容器內(nèi)，大約40 mL[13]。樣品以及標(biāo)準(zhǔn)樣送往中國(guó)科學(xué)院蘭州化學(xué)物理研究所進(jìn)行檢測(cè)，運(yùn)用電感耦合等離子體發(fā)射光譜儀ICP-OES測(cè)定重金屬的含量[14]，As、Cd、Cr、Cu、Pb 5種重金屬的檢出限分別為0.42、0.012、0.018、0.014、0.42 mg/L。

圖3 大夏河流域采樣點(diǎn)分布圖Fig.3 Distribution map of sampling points in Daxia river basin

1.4 重金屬運(yùn)移模型

1.4.1SCS-CN運(yùn)移模型

流域污染物的運(yùn)移研究直接影響環(huán)境科學(xué)的進(jìn)展，通常采用數(shù)學(xué)模型預(yù)測(cè)污染物的運(yùn)移規(guī)律[15]。掌握流域重金屬運(yùn)移規(guī)律可以有效地控制與治理環(huán)境污染[16]。國(guó)外對(duì)于城市小流域進(jìn)行地表徑流污染特征的提取應(yīng)用較為廣泛[17]，國(guó)內(nèi)相關(guān)資料較少。流域土壤的質(zhì)地對(duì)重金屬的遷移有一定的影響[18]。SCS-CN模型由美國(guó)農(nóng)業(yè)部土壤保持局(USDA SCS)提出，由于對(duì)輸入數(shù)據(jù)量要求不高、模擬精度高而被許多國(guó)家和地區(qū)廣泛應(yīng)用[19]。本文選取大夏河流域及其支流，以As、Cd、Cr、Cu、Pb 5種重金屬為估算對(duì)象，結(jié)合不同土地利用類型，統(tǒng)計(jì)流域降雨徑流的SCS-CN模型因子，計(jì)算不同土地利用類型下的地表徑流量，估算大夏河流域的重金屬運(yùn)移能力[20]。

流域重金屬地表徑流運(yùn)移量公式為

(1)

其中

S=25 400/CN-254

式中γi——第i類重金屬的運(yùn)移量，t/a

αi——第i類重金屬的水溶出率

Aj——第j個(gè)水文響應(yīng)單元(HRU)的面積，km2

P——研究區(qū)域年降雨量，mm/a

CN——徑流曲線數(shù)

S——徑流開(kāi)始前的潛在最大滯留量，mm/a

根據(jù)大夏河流域多年降雨量，把大夏河流域按照子流域進(jìn)行劃分，統(tǒng)計(jì)各個(gè)子流域的年降雨量，并根據(jù)子流域劃分圖，確定子流域的流域面積，如表1所示。

表1 子流域統(tǒng)計(jì)結(jié)果Tab.1 Sub-basin statistical results

1.4.2分布式運(yùn)移模型

定義河網(wǎng)分形維數(shù)與地表徑流量的函數(shù)關(guān)系，可以得到大夏河各個(gè)子流域用分維值預(yù)測(cè)的徑流量[21]?？紤]到重金屬運(yùn)移量模型的單位換算，應(yīng)用徑流深這一概念消除單位不統(tǒng)一的影響，徑流深指在某一時(shí)間段內(nèi)通過(guò)流域指定斷面的徑流量除以該斷面以上的流域面積的值[22]，公式為

(2)

式中R——徑流深，mm

Q——年徑流量，m3

F——該斷面以上的流域面積，km2

根據(jù)分維值求得的年徑流量，將流域重金屬運(yùn)移量公式定義為

(3)

式中Sj——第j個(gè)子流域的面積，km2

考慮到徑流深的特點(diǎn)，將重金屬運(yùn)移模型公式修改為

(4)

式中Qj——第j個(gè)流域分維值模擬的徑流量，m3

1.5 模型因子

1.5.1土地利用類型分類

在ArcGIS 10.4下，大夏河流域的土地利用類型經(jīng)過(guò)遙感監(jiān)督分類解析，主要選取了旱地、林地、草地、水域、居民用地和沼澤地，如圖4所示。

圖4 土地利用類型圖Fig.4 Land use type map

1.5.2徑流曲線數(shù)

徑流曲線數(shù)CN取決于研究區(qū)域的土地利用類型和土壤類型，本文參照《美國(guó)國(guó)家工程手冊(cè)》中的CN值[23]，結(jié)合流域多年降雨量和土壤濕潤(rùn)程度，確定CN值，如表2所示。參照美國(guó)水土保持局手冊(cè)，選定CN為A類土壤組。

1.5.3網(wǎng)格法

應(yīng)用網(wǎng)格法計(jì)算大夏河及子流域的分維數(shù)，網(wǎng)格法基本原理為使用不同邊長(zhǎng)的正方形網(wǎng)格覆蓋被

表2 不同土地利用類型下的CNTab.2 Value of CN under different land use types

測(cè)流域等線體，當(dāng)改變正方形的邊長(zhǎng)時(shí)，被測(cè)線體的網(wǎng)格數(shù)必然會(huì)出現(xiàn)相應(yīng)的變化[24-25]，因此，擬合方程式為

lnN(a)=-Dlna+C

(5)

圖5 相似流域圖像Fig.5 Similar watershed image

式中a——邊長(zhǎng)N(a)——網(wǎng)格數(shù)

C——比例常數(shù)

D——流域的分維數(shù)

1.5.4地勢(shì)起伏比的計(jì)算

徑流是流域地貌形成的外營(yíng)力之一，并且參與地殼中的地球化學(xué)過(guò)程，不僅影響植物的生長(zhǎng)和湖泊、沼澤的形成，還影響土壤的發(fā)育程度。流域分維值又能夠很好地反映流域地貌形態(tài)[26]，應(yīng)用流域分維值擬合徑流存在可行性。應(yīng)用ArcGIS計(jì)算和大夏河流域形態(tài)相似的流域水文站的分維值，并計(jì)算該水文站的徑流量。

由于研究區(qū)域所處的地勢(shì)起伏度比較大，對(duì)地表徑流量的影響也比較大，所以需要從ArcGIS中獲取該研究區(qū)域最大與最小高程，然后計(jì)算地勢(shì)起伏比，從而將分維值乘以地勢(shì)起伏比進(jìn)行歸一化，公式為

(6)

式中Rm——地勢(shì)起伏比

Hmax——研究區(qū)域的最大高程

Hmin——研究區(qū)域的最小高程

1.5.5相似流域的識(shí)別模型

1.5.5.1結(jié)構(gòu)相似性

先分流域生成圖像的像素分布函數(shù)，再通過(guò)分布函數(shù)的相似性，得到兩幅圖像的相似度。若用1表示圖像中被點(diǎn)出的點(diǎn)，用0表示圖像中未被點(diǎn)出的點(diǎn)，那么，該圖像其實(shí)就是一個(gè)由1和0組成的0-1矩陣。該圖像轉(zhuǎn)換成0-1矩陣后可應(yīng)用Matlab進(jìn)行操作。因此，各個(gè)圖像的結(jié)構(gòu)相似問(wèn)題就轉(zhuǎn)換為比較矩陣的不同區(qū)域之間的相似性問(wèn)題[27]。大夏河流域以及與大夏河流域比較相似的9個(gè)流域的圖像如圖5所示。

為了比較兩幅圖像的相似程度，可以定義任一圖像的像素分布函數(shù)。設(shè)一幅圖像F的 0-1矩陣為M。定義向量

(7)

式中vk——圖像的豎直像素分布向量

mjk——矩陣M中j行k列的元素

由此，定義函數(shù)

(8)

式中k——矩陣的維數(shù)

f(k/n)——圖像F的像素分布函數(shù)

根據(jù)式(8)的定義，可得到函數(shù)f(x)在x=k/n點(diǎn)的值。因?yàn)槊糠鶊D像的寬度為230像素，高度為230像素，所以1≤k≤230。

T=(u1,u2,…,u230)

(9)

式中T——向量F的特征向量

uk——函數(shù)f(x)在區(qū)間(0.01(k-1),0.01k)的平均值

假設(shè)兩幅圖像為F1和F2，它們對(duì)應(yīng)的0-1 矩陣為M1、M2，對(duì)應(yīng)的特征向量分別為T(mén)1、T2，由此定義

(10)

式中 〈·,·〉——兩個(gè)向量的內(nèi)積

α——兩幅圖像的結(jié)構(gòu)相似度

如果α大于0.99，說(shuō)明兩幅圖像結(jié)構(gòu)相似。

1.5.5.2特征相似性

相似流域的識(shí)別不僅在結(jié)構(gòu)上比較相似，在流域面積、河網(wǎng)長(zhǎng)度、年平均氣溫、年平均降水量等某些特征指標(biāo)方面也具有相似性[28]。表3對(duì)比了大夏河流域及與大夏河相似的9個(gè)流域的一些特征參數(shù)。

表3 相似流域的特征對(duì)比Tab.3 Characteristic comparison of similar basins

2 結(jié)果與分析

2.1 水溶出率

采用BCR連續(xù)分級(jí)提取法[29]，將水溶態(tài)和可交換態(tài)合并的酸溶態(tài)含量(質(zhì)量比)作為水溶出率計(jì)算因子，在提取過(guò)程中可保持穩(wěn)定，準(zhǔn)確地描述污染物的運(yùn)移能力。以5組土壤重金屬樣本作為測(cè)試提取對(duì)象，得到平均提取量和水溶出率如表4所示。

表4 BCR提取水溶態(tài)含量和水溶出率結(jié)果Tab.4 BCR extraction results of water soluble content state and water dissolution rate

2.2 重金屬含量特征

由于流域位于甘肅省境內(nèi)，因此采用按行政區(qū)域劃分的甘肅省土壤背景值的幾何平均值作為參考[30]，并對(duì)數(shù)據(jù)進(jìn)行描述性統(tǒng)計(jì)分析，結(jié)果如表5所示。

表5 重金屬含量描述性統(tǒng)計(jì)分析Tab.5 Descriptive statistical analysis of heavy metal content

從表5可以看出，5種重金屬元素的平均含量從大到小依次為：Cr、As、Cu、Pb、Cd。

2.3 分維值模擬徑流量結(jié)果

閾值對(duì)應(yīng)的是區(qū)域的匯流量，也就是匯入該區(qū)域的柵格流量數(shù)，本文在ArcGIS 10.4中對(duì)該水系的矢量圖取5 000的閾值，并對(duì)該矢量圖進(jìn)行要素轉(zhuǎn)柵格、網(wǎng)格分析，在輸入像元大小(正方形網(wǎng)格邊長(zhǎng))里填入10～500不同的10組數(shù)值，得到不同邊長(zhǎng)所對(duì)應(yīng)的柵格圖，統(tǒng)計(jì)柵格圖的網(wǎng)格數(shù)，對(duì)求得的數(shù)值取對(duì)數(shù)進(jìn)行線性擬合，生成擬合方程為lnN(a)=-1.016 3lna+14.331，決定系數(shù)R2=0.999 8，擬合方程系數(shù)的斜率即為流域的分維值，故大夏河流域的分維值為1.016 3，小于1.6，判定大夏河流域地貌處于幼年期，此階段流域水系發(fā)育不充分，地面較為完整。同理，分別計(jì)算7個(gè)子流域的分維值。

與大夏河流域比較相似的9個(gè)流域，其分維值基本與大夏河各子流域分維值相近，且某些特征指標(biāo)也較為相近，所以求得相似流域分維值與徑流量如表6所示。

對(duì)這9個(gè)水文站的歸一化分維值和徑流量做線性擬合，得到擬合曲線如圖6所示。

由圖6可知，相似流域分維值和徑流量擬合方程為：y=-8.3×109x+6.1×109，決定系數(shù)R2=0.761 9，說(shuō)明流域分維值和徑流量相關(guān)性比較顯著，該研究方法具有可行性。

表6 相似流域分維值和徑流量Tab.6 Similar basin fractal dimension and runoff

圖6 歸一化的分維值和徑流量擬合曲線Fig.6 Normalized fractal dimension and runoff fitting curve

根據(jù)擬合曲線求大夏河各個(gè)子流域的徑流量，如表7所示。

從全國(guó)水雨情網(wǎng)站獲取水文站的流量數(shù)據(jù)，以及從甘肅省水利廳官網(wǎng)下載《甘肅省水資源公報(bào)》，查詢到折橋和夏河水文站的平均年徑流量分別為7.265億m3和3.817億m3(2016年和2017年)，模擬的年徑流量與實(shí)測(cè)數(shù)據(jù)基本吻合，說(shuō)明用分維值模擬的徑流數(shù)據(jù)可以真實(shí)反映流域徑流量。

表7 大夏河流域模擬徑流量Tab.7 Simulating runoff in Daxia river basin

2.4 流域相似度的計(jì)算

先對(duì)大夏河流域以及與大夏河流域比較相似的9個(gè)流域的圖像進(jìn)行旋轉(zhuǎn)與裁剪，使每幅流域圖像的寬度為230像素，高度為230像素。進(jìn)而在Matlab R2014a中，把各個(gè)流域的圖像導(dǎo)入并轉(zhuǎn)換成0-1矩陣，再把數(shù)據(jù)導(dǎo)出，最后通過(guò)相似流域的識(shí)別模型進(jìn)行計(jì)算。

根據(jù)相似流域的識(shí)別模型可以求得結(jié)構(gòu)相似度α，如表8所示。

表8 流域相似度Tab.8 Watershed similarity

從表8可以看出，多個(gè)流域圖塊與參考圖的結(jié)構(gòu)相似度α均大于0.99(其中相似度為1即參考圖所在的區(qū)域，是參考圖與自身相比的結(jié)果)，說(shuō)明大夏河流域以及與大夏河流域比較相似的9個(gè)流域具有結(jié)構(gòu)相似性。

通過(guò)流域相似度的驗(yàn)證與分析，大夏河流域以及與大夏河流域相似的9個(gè)流域在結(jié)構(gòu)上是相似的，并且根據(jù)表3可以看出，這些流域在不同尺度下的年平均氣溫、年平均降水量、流域面積與河網(wǎng)長(zhǎng)度比值也具有相似性。

2.5 兩種模型對(duì)比

根據(jù)SCS-CN模型估算出大夏河流域的重金屬運(yùn)移量，如表9所示。

表9 基于SCS-CN模型估算的重金屬運(yùn)移量Tab.9 Heavy metal migration estimated based on SCS-CN model t/a

SCS-CN模型的優(yōu)點(diǎn)在于結(jié)構(gòu)簡(jiǎn)單、輸入?yún)?shù)少，缺點(diǎn)在于數(shù)據(jù)需要長(zhǎng)期持續(xù)性觀測(cè)。從表9可以看出，大夏河各個(gè)子流域的重金屬運(yùn)移量，As、Cd、Cr、Cu、Pb 5種重金屬的總運(yùn)移量為0.626 8、0.197 3、5.805 1、1.596 0、1.714 1 t/a。

根據(jù)分布式預(yù)測(cè)模型估算出大夏河流域的重金屬運(yùn)移量，如表10所示。

分布式運(yùn)移模型優(yōu)點(diǎn)在于從機(jī)理上去除了重金屬河道滯留量和土壤入滲保留量，可以結(jié)合少量采樣數(shù)據(jù)對(duì)遠(yuǎn)程遙感進(jìn)行估算，缺點(diǎn)在于需要根據(jù)實(shí)際區(qū)域的土壤顆粒類型和水土保持因子來(lái)做進(jìn)一步的剔除。從表10可以看出，大夏河各個(gè)子流域的重金屬運(yùn)移量，As、Cd、Cr、Cu、Pb 5種重金屬的總運(yùn)移量為0.671 0、0.209 9、6.281 6、1.746 5、1.837 7 t/a。

表10 基于分布式模型估算的重金屬運(yùn)移量Tab.10 Heavy metal migration based on distributed model estimation t/a

對(duì)比表9、10可以看出，應(yīng)用分維值估算的徑流量分布式模型求取的重金屬運(yùn)移量和應(yīng)用SCS-CN模型求取的重金屬運(yùn)移量結(jié)果基本一致，說(shuō)明分布式重金屬運(yùn)移模型相對(duì)合理。

3 結(jié)論

(1) Cu、Pb元素的實(shí)測(cè)平均含量分別為20.25、15.37 mg/kg，分別低于甘肅省Cu、Pb含量的地方背景值22.5、17.9 mg/kg，說(shuō)明這兩種元素沒(méi)有超標(biāo)。As、Cr、Cd元素的實(shí)測(cè)平均含量分別為20.43、94.48、1.43 mg/kg，分別高于甘肅省As、Cr、Cd含量的地方背景值11.7、69.3、0.110 6 mg/kg，說(shuō)明這3種元素在大夏河流域存在不同程度的歷史累積。

(2)采用網(wǎng)格法對(duì)大夏河流域及其子流域進(jìn)行分析，其水系分維值為1.016 3。根據(jù)相似流域水文站的分維值和徑流數(shù)據(jù)，計(jì)算分維值與徑流量之間的線性關(guān)系，擬合曲線求出大夏河子流域的徑流數(shù)據(jù)，對(duì)大夏河流域及與大夏河流域比較相似的9個(gè)流域的相似性進(jìn)行了建模驗(yàn)證。

(3)基于SCS-CN和分布式運(yùn)移兩種模型的估算結(jié)果基本一致，但總體上應(yīng)用分布式運(yùn)移模型估算重金屬的運(yùn)移量偏高，比例較為合理，尤其對(duì)于遠(yuǎn)程遙感估算更具有實(shí)際應(yīng)用前景。