流域土壤重金屬分布式運移模型研究

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

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

0 引言

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

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

1 材料與方法

1.1 研究區概況

大夏河屬于黃河上游的一級支流，地處青藏高原與黃土高原過渡帶(東經102°2′～103°23′，北緯34°51′～35°48′)。大夏河經夏河縣城東北流入，出土門關后進入臨夏盆地，自西向東橫貫臨夏市、至折橋轉北向東鄉縣泄湖峽，在康家灣喇嘛川塔張村注入劉家峽水庫[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 流域分級與劃分

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

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

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

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

1.3 數據處理與分析

2016—2018年，使用手持式GPS記錄位置，采集大夏河流域周邊土樣和水樣，以全年存在穩定徑流為依據，共26個采樣點，采樣位置具有代表性，見圖3。在野外將樣品放入自封袋中，編號后帶回實驗室，對土壤樣品進行四分法取樣后風干，研磨過200目尼龍篩，放入對應編號的樣品容器內，對流域水樣搖勻，也裝入對應編號的樣品容器內，大約40 mL[13]。樣品以及標準樣送往中國科學院蘭州化學物理研究所進行檢測，運用電感耦合等離子體發射光譜儀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運移模型

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

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

(1)

其中

S=25 400/CN-254

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

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

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

P——研究區域年降雨量，mm/a

CN——徑流曲線數

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

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

表1 子流域統計結果Tab.1 Sub-basin statistical results

1.4.2分布式運移模型

定義河網分形維數與地表徑流量的函數關系，可以得到大夏河各個子流域用分維值預測的徑流量[21]。考慮到重金屬運移量模型的單位換算，應用徑流深這一概念消除單位不統一的影響，徑流深指在某一時間段內通過流域指定斷面的徑流量除以該斷面以上的流域面積的值[22]，公式為

(2)

式中R——徑流深，mm

Q——年徑流量，m3

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

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

(3)

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

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

(4)

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

1.5 模型因子

1.5.1土地利用類型分類

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

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

1.5.2徑流曲線數

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

1.5.3網格法

應用網格法計算大夏河及子流域的分維數，網格法基本原理為使用不同邊長的正方形網格覆蓋被

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

測流域等線體，當改變正方形的邊長時，被測線體的網格數必然會出現相應的變化[24-25]，因此，擬合方程式為

lnN(a)=-Dlna+C

(5)

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

式中a——邊長N(a)——網格數

C——比例常數

D——流域的分維數

1.5.4地勢起伏比的計算

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

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

(6)

式中Rm——地勢起伏比

Hmax——研究區域的最大高程

Hmin——研究區域的最小高程

1.5.5相似流域的識別模型

1.5.5.1結構相似性

先分流域生成圖像的像素分布函數，再通過分布函數的相似性，得到兩幅圖像的相似度。若用1表示圖像中被點出的點，用0表示圖像中未被點出的點，那么，該圖像其實就是一個由1和0組成的0-1矩陣。該圖像轉換成0-1矩陣后可應用Matlab進行操作。因此，各個圖像的結構相似問題就轉換為比較矩陣的不同區域之間的相似性問題[27]。大夏河流域以及與大夏河流域比較相似的9個流域的圖像如圖5所示。

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

(7)

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

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

由此，定義函數

(8)

式中k——矩陣的維數

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

根據式(8)的定義，可得到函數f(x)在x=k/n點的值。因為每幅圖像的寬度為230像素，高度為230像素，所以1≤k≤230。

T=(u1,u2,…,u230)

(9)

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

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

假設兩幅圖像為F1和F2，它們對應的0-1 矩陣為M1、M2，對應的特征向量分別為T1、T2，由此定義

(10)

式中 〈·,·〉——兩個向量的內積

α——兩幅圖像的結構相似度

如果α大于0.99，說明兩幅圖像結構相似。

1.5.5.2特征相似性

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

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

2 結果與分析

2.1 水溶出率

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

表4 BCR提取水溶態含量和水溶出率結果Tab.4 BCR extraction results of water soluble content state and water dissolution rate

2.2 重金屬含量特征

由于流域位于甘肅省境內，因此采用按行政區域劃分的甘肅省土壤背景值的幾何平均值作為參考[30]，并對數據進行描述性統計分析，結果如表5所示。

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

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

2.3 分維值模擬徑流量結果

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

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

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

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

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

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

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

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

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

2.4 流域相似度的計算

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

根據相似流域的識別模型可以求得結構相似度α，如表8所示。

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

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

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

2.5 兩種模型對比

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

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

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

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

分布式運移模型優點在于從機理上去除了重金屬河道滯留量和土壤入滲保留量，可以結合少量采樣數據對遠程遙感進行估算，缺點在于需要根據實際區域的土壤顆粒類型和水土保持因子來做進一步的剔除。從表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可以看出，應用分維值估算的徑流量分布式模型求取的重金屬運移量和應用SCS-CN模型求取的重金屬運移量結果基本一致，說明分布式重金屬運移模型相對合理。

3 結論

(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)采用網格法對大夏河流域及其子流域進行分析，其水系分維值為1.016 3。根據相似流域水文站的分維值和徑流數據，計算分維值與徑流量之間的線性關系，擬合曲線求出大夏河子流域的徑流數據，對大夏河流域及與大夏河流域比較相似的9個流域的相似性進行了建模驗證。

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