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

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

0 引言

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

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

1 材料與方法

1.1 研究區(qū)概況

大夏河屬于黃河上游的一級支流，地處青藏高原與黃土高原過渡帶(東經(jīng)102°2′～103°23′，北緯34°51′～35°48′)。大夏河經(jīng)夏河縣城東北流入，出土門關(guān)后進入臨夏盆地，自西向東橫貫臨夏市、至折橋轉(zhuǎn)北向東鄉(xiāng)縣泄湖峽，在康家灣喇嘛川塔張村注入劉家峽水庫[6]。主河道總長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īng)緯度，從地理空間數(shù)據(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

本文首先對大夏河流域選取合適的閾值[8-9]進行河網(wǎng)水系提取，然后根據(jù)Strahler和Shreve這兩種流域分級法[10]，并結(jié)合ArcMAP軟件Spatial Analyst—水文分析—河網(wǎng)分級工具，對大夏河流域進行分級處理。而流域的劃分首先需要確定各個子流域的出水口位置，然后根據(jù)水流流向[11-12]判斷流入該出水口的所有上游區(qū)域，即為子流域，在ArcGIS軟件中把大夏河流域劃分成7個子流域，如圖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個采樣點，采樣位置具有代表性，見圖3。在野外將樣品放入自封袋中，編號后帶回實驗室，對土壤樣品進行四分法取樣后風(fēng)干，研磨過200目尼龍篩，放入對應(yīng)編號的樣品容器內(nèi)，對流域水樣搖勻，也裝入對應(yīng)編號的樣品容器內(nèi)，大約40 mL[13]。樣品以及標準樣送往中國科學(xué)院蘭州化學(xué)物理研究所進行檢測，運用電感耦合等離子體發(fā)射光譜儀ICP-OES測定重金屬的含量[14]，As、Cd、Cr、Cu、Pb 5種重金屬的檢出限分別為0.42、0.012、0.018、0.014、0.42 mg/L。

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

1.4 重金屬運移模型

1.4.1SCS-CN運移模型

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

流域重金屬地表徑流運移量公式為

(1)

其中

S=25 400/CN-254

式中γi——第i類重金屬的運移量，t/a

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

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

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

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

S——徑流開始前的潛在最大滯留量，mm/a

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

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

1.4.2分布式運移模型

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

(2)

式中R——徑流深，mm

Q——年徑流量，m3

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

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

(3)

式中Sj——第j個子流域的面積，km2

考慮到徑流深的特點，將重金屬運移模型公式修改為

(4)

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

1.5 模型因子

1.5.1土地利用類型分類

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

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

1.5.2徑流曲線數(shù)

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

1.5.3網(wǎng)格法

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

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

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

lnN(a)=-Dlna+C

(5)

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

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

C——比例常數(shù)

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

1.5.4地勢起伏比的計算

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

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

(6)

式中Rm——地勢起伏比

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

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

1.5.5相似流域的識別模型

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

先分流域生成圖像的像素分布函數(shù)，再通過分布函數(shù)的相似性，得到兩幅圖像的相似度。若用1表示圖像中被點出的點，用0表示圖像中未被點出的點，那么，該圖像其實就是一個由1和0組成的0-1矩陣。該圖像轉(zhuǎn)換成0-1矩陣后可應(yīng)用Matlab進行操作。因此，各個圖像的結(jié)構(gòu)相似問題就轉(zhuǎn)換為比較矩陣的不同區(qū)域之間的相似性問題[27]。大夏河流域以及與大夏河流域比較相似的9個流域的圖像如圖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點的值。因為每幅圖像的寬度為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，它們對應(yīng)的0-1 矩陣為M1、M2，對應(yīng)的特征向量分別為T1、T2，由此定義

(10)

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

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

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

1.5.5.2特征相似性

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

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

2 結(jié)果與分析

2.1 水溶出率

采用BCR連續(xù)分級提取法[29]，將水溶態(tài)和可交換態(tài)合并的酸溶態(tài)含量(質(zhì)量比)作為水溶出率計算因子，在提取過程中可保持穩(wěn)定，準確地描述污染物的運移能力。以5組土壤重金屬樣本作為測試提取對象，得到平均提取量和水溶出率如表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]，并對數(shù)據(jù)進行描述性統(tǒng)計分析，結(jié)果如表5所示。

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

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

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

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

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

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

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

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

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

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

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

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

2.4 流域相似度的計算

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

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

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

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

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

2.5 兩種模型對比

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

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

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

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

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

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

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

3 結(jié)論

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

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

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