易 琦，王瑞芳，趙筱青**，李嘉旺，竇小東

(1.云南大學 地球科學學院，云南 昆明 650500；2.云南省氣象服務中心，云南 昆明 650034)

重金屬污染物毒性強、不易降解，其濃度達到一定閾值，不僅會危害水體及水生生態(tài)系統(tǒng)，還能經(jīng)過食物鏈傳遞，危害人體健康[1-3].河流重金屬元素的轉(zhuǎn)化主要為水體中溶解態(tài)、懸浮態(tài)以及底泥中沉積態(tài)三者之間的相互轉(zhuǎn)化[4-6].目前，對于大型河流水體中重金屬遷移模擬多采用二、三維模型，中小型河流多采用一、二維模型.黃本生等[7]建立河流重金屬隨水?懸浮物?底泥遷移轉(zhuǎn)化耦合模型.武萬國[8]建立了湘江長株潭河段二維水動力?水質(zhì)模型.無論哪一種模型，確定水體重金屬污染物的遷移轉(zhuǎn)化系數(shù)，均是確定水體重金屬污染物遷移規(guī)律的關(guān)鍵因素.沘江為中小河流，是瀾滄江上游的一級支流，也是沿岸居民生產(chǎn)生活的主要水源，但由于近年來重金屬對河水的污染，河水使用率急劇降低，導致沿岸的蘭坪、云龍兩縣出現(xiàn)不同程度的水質(zhì)性缺水[9].因此，摸清沘江水體重金屬污染物的遷移轉(zhuǎn)化規(guī)律，對沘江水體重金屬污染治理具有積極意義.本文在枯水期對沘江水體采樣分析基礎上，對河水中重金屬污染物的綜合衰減系數(shù)進行率定，從而建立沘江水體重金屬污染物的遷移轉(zhuǎn)化方程，研究結(jié)果可為沘江流域環(huán)境污染防治及生態(tài)恢復等提供科學依據(jù).

1 研究區(qū)概況

沘江屬瀾滄江水系，流域面積約2 400 km2，干流長約150 km，總體流向由北向南，途經(jīng)蘭坪、云龍兩縣后匯入瀾滄江[10-11].金頂鉛鋅礦位于沘江上游河段的東部，其儲量居中國已探明鉛鋅礦儲量的首位.礦區(qū)三面環(huán)山，大部分可以露天開采，雨季時非常容易匯水，導致大量的礦渣、碎石等被沖入河道，沘江水體已魚蝦絕跡，污染嚴重，無法飲用和灌溉，沘江水質(zhì)性缺水問題嚴重[12-14].

2 研究方法

2.1 樣品采集與分析2010 年1 月底，按照礦區(qū)對沘江影響大小，以距離礦區(qū)遠近為原則，距離越遠樣點越稀疏，反之樣點越密集，共取得沘江水樣15 個（圖1）.其中，礦區(qū)上游采樣點3 個，基本不受礦業(yè)活動影響，取其均值以反映礦業(yè)開采前的沘江水質(zhì)狀況，作為水體的背景值；礦區(qū)下游采樣點12 個，其中，1～11 號樣點分別位于蘭坪縣和云龍縣境內(nèi)，12 號樣點為瀾滄江與沘江交匯處水樣.水樣于岸邊采集并貯存于2.5 L 的高密度聚乙烯瓶中.水樣的重金屬質(zhì)量濃度由昆明礦產(chǎn)資源監(jiān)督檢測中心檢測分析，檢出限分別為：ρ（Zn）<1.0 μg/L，ρ（Pb）<0.01 μg/L，ρ（As）<0.01 μg/L.

圖1 研究區(qū)域及采樣點分布Fig.1 Research area and distribution of sampling points

2.2 污染物遷移擴散模擬方法

2.2.1 水體重金屬遷移擴散模型的選取 天然中小河流多為淺寬河流，重金屬污染物進入水體后很快在垂向均勻混合[6]，因此在考慮河流重金屬污染物質(zhì)量濃度變化時只考慮橫向（y）與縱向（x）方向即二維狀態(tài)即可.而對于天然寬度較窄的小河流而言，可以近似地認為重金屬污染物進入水體后迅速在橫向（y）與垂向（z）方向均勻混合，其水體中重金屬污染物遷移轉(zhuǎn)化僅考慮縱向（x）方向即一維狀態(tài)即可，沘江即屬此類河流.

綜合考慮水體中重金屬遷移轉(zhuǎn)化的水動力學特性，及重金屬在水體?懸浮物?底泥間的吸附、解吸、沉降、再懸浮過程[15]，同時，結(jié)合重金屬自身的化學轉(zhuǎn)化因素的影響，可知沘江水體重金屬污染物遷移的一維模擬方程可用下式表示：

式中，D為分子擴散系數(shù)（單位：m2/s）；ρ為溶質(zhì)的質(zhì)量濃度（單位：mg/L）；X為河流縱向上距排污口附近距離（單位：m）；u為河水的平均流速（單位：m/s）；K為綜合衰減系數(shù)（單位：s?1）；ρ0為礦區(qū)下游排污口附近重金屬質(zhì)量濃度值（單位：mg/L）；ρb為沘江河流水體重金屬背景值（單位：mg/L）.

若不考慮水體重金屬背景值的影響，式（1）的解析解為：

若考慮到背景值的影響則公式（2）、（3）可分別表示為：

本次沘江實地采樣，1 號采樣點是礦區(qū)下游水體中重金屬與沘江水體重金屬背景值的總含量，因此ρ0的值已經(jīng)包含ρb的影響.因此，對于沘江水體重金屬遷移轉(zhuǎn)化的的解析解應采用公式（2）、（3）模擬，鑒于公式（3）中縱向擴散系數(shù)D的確定需要考慮河水深度及比降的影響，而本次采樣并未測得沘江深度及比降的具體數(shù)據(jù)，同時河流的比降及深度對重金屬遷移的影響也主要是通過流速實現(xiàn).因此，本次沘江水體重金屬遷移轉(zhuǎn)化的解析解采用公式（2）進行.

2.2.2 模型參數(shù)的率定及檢驗 模型中邊界值的選取方法是取礦區(qū)下游的1 號樣點的Zn、Pb、As 的質(zhì)量濃度作為連續(xù)面源污染的初始值ρ0.

模型中綜合衰減系數(shù)K的率定一般包括經(jīng)驗公式法、實驗法、試錯法[16].經(jīng)驗公式法是根據(jù)前人總結(jié)的一些經(jīng)驗關(guān)系式或近似計算方法確定某一參數(shù)的取值.所謂實驗法是通過相關(guān)實驗確定某一參數(shù)的具體值，而試錯法則是通過查閱相關(guān)文獻，確定該參數(shù)的范圍或該參數(shù)的取值大概在哪個數(shù)值附近，進而通過嘗試不同的賦值即試錯，并對賦值后的結(jié)果檢驗，若達到預期要求或通過檢驗，則所賦值即為該參數(shù)值，若達不到預期要求或通不過檢驗，則要繼續(xù)賦予其它的值，直到取得滿意效果或通過檢驗為止.鑒于目前并無成熟的綜合衰減系數(shù)K的經(jīng)驗公式，且由于自身實驗條件的限制，本文采取試錯法來確定K的取值，并運用PASW18.0軟件不斷對K的不同試錯值所生成的各樣點模擬值與實測值進行配對樣本T檢驗，直到通過檢驗進而確定K值的范圍，并將配對T檢驗效果最理想時所用的K值作為K的首選值.

配對T檢驗是對相同總體在兩種不同的條件下進行相同研究，對比兩種條件下總體均值是否存在顯著差異的方法.因其計算簡便、檢驗性能高效，因而得到廣泛應用[17].配對樣本T檢驗的零假設H0認為，兩總體均值之間沒有顯著差異.檢驗思路為：首先計算每對觀察值的差值，得到差值序列；再計算差值的均值；最后檢驗均值是否和0 存在顯著差異.若存在顯著差異，則認為兩總體均值間有顯著差異，反之，沒有顯著差異.

應用PASW18.0 軟件分別進行Zn、Pb、As 的實測值與不同K值下其模擬值的配對樣本T檢驗，并根據(jù)T分布表給出T值所對應的概率值P.依據(jù)為若P<顯著水平α，則拒絕零假設，認為兩總體均值有顯著差異；若P>顯著水平α，則不能拒絕零假設，認為兩總體均值不存在顯著差異，判定K值范圍[18].

2.2.3 基于MATLAB 對沘江水體重金屬遷移全程模擬 MATLAB 軟件即矩陣實驗室，作為科技應用軟件，其集程序設計語言、圖形處理以及數(shù)學于一體[19-20]，把科學計算、結(jié)果的可視化和編程集中在一個非常方便的環(huán)境中，該軟件已經(jīng)在許多領域得到了廣泛的應用.宋新山等[21]對應用該軟件構(gòu)建地表水環(huán)境系統(tǒng)的模擬模型做了詳細的實例演示，效果較好.因此，本文選用該軟件對沘江水體重金屬遷移進行全程模擬.

3 結(jié)果與分析

3.1 沘江水體中Zn 的模擬結(jié)果與實測數(shù)據(jù)的配對樣本T 檢驗模型中各樣點距1 號樣點的距離由采樣時GPS 所測得的各樣點的地理坐標經(jīng)Arcgis 軟件求得.由于天然河流的水文條件及周邊環(huán)境各不相同，因此目前并無普遍適用的K值范圍，但可以確認在有連續(xù)穩(wěn)定污染源排放的天然河流中K>0.本文嘗試性選取步長0.001 s?1逐次遞增K值，運用公式（2）計算2～12 號采樣點Zn 的模擬值.PASW18.0 軟件對K值取0.001～0.04 s?1時，2～12 號采樣點Zn 的模擬值與實測值的配對T檢驗結(jié)果顯示（圖2（a））：當K值介于0.013～0.017 s?1時，P>α=0.05，模擬值與實測值不存在顯著差異，即當沘江水體中Zn 的綜合衰減系數(shù)K值介于0.013～0.017 s?1時，公式（2）可以較好地反映沘江水體中Zn 的遷移轉(zhuǎn)化規(guī)律；當K=0.015 s?1時，P=0.762，公式（2）的模擬值最接近真實值.因此，沘江水體中Zn 的綜合衰減系數(shù)K的最優(yōu)值取0.015 s?1.

3.2 沘江水體中Pb 的模擬結(jié)果與實測數(shù)據(jù)的配對樣本T 檢驗鑒于10～12 號采樣點間Pb 的遷移轉(zhuǎn)化受上游攔河壩的影響較大的事實，沘江水體中Pb 的綜合衰減系數(shù)K值采取分段率定，即分別率定2～9 號采樣點與10～12 號采樣點的綜合衰減系數(shù)K值.仍然選取步長0.001 逐次遞增K值，并運用公式（2）分別計算2～9 號及10～12 號采樣點Pb 的模擬值.

PASW18.0 軟件對K值取0.001～0.04 s?1時，2～9 號采樣點Pb 的模擬值與實測值的配對T檢驗結(jié)果顯示（圖2（b））：當K值介于0.003～0.009 s?1時，P>α=0.05，模擬值與實測值不存在顯著差異，即當沘江水體2～9 號采樣點間Pb 的綜合衰減系數(shù)K值介于0.003～0.009 s?1時，公式（2）可以較好地反映沘江水體中2～9 號采樣點間Pb 的遷移轉(zhuǎn)化規(guī)律；當K=0.006 s?1時，P=0.692，公式（2）的模擬值最接近真實值.因此，沘江水體中2～9 號采樣點間Pb 的綜合衰減系數(shù)K的最優(yōu)值取0.006 s?1.

PASW18.0 軟件對K值取0.001～0.02 s?1時，10～12 號采樣點Pb 的模擬值與實測值的配對T檢驗結(jié)果顯示（圖2（c））：當K值介于0.002～0.008 s?1時，P>α=0.05，模擬值與實測值不存在顯著差異，即當沘江水體10～12 號采樣點間Pb 的綜合衰減系數(shù)K值介于0.002～0.008 s?1時，公式（2）可以較好地反映沘江水體中10～12 號采樣點間Pb 的遷移轉(zhuǎn)化規(guī)律；當K=0.004 s?1時，P=0.667，公式（2）的模擬值最接近真實值.因此，沘江水體中10～12 號采樣點間Pb 的綜合衰減系數(shù)K的最優(yōu)值取0.004 s?1.

圖2 不同K 值時2～12 號采樣點模擬值與實測值配對T 檢驗結(jié)果Fig.2 PairedT test results of simulated and measured values at samples 2-12 at differentK values

3.3 沘江水體中As 的模擬結(jié)果與實測數(shù)據(jù)的配對樣本T 檢驗10、11、12 號采樣點間As 的遷移轉(zhuǎn)化亦受上游攔河壩的影響較大，沘江中As 的綜合衰減系數(shù)亦采取分段率定，即分別率定2～9 號采樣點與10～12 號采樣點的綜合衰減系數(shù)K.同樣選取步長0.001 逐次遞增K值，并運用公式（2）分別計算2～9 號及10～12 號采樣點As 的模擬值.

PASW18.0 軟件對K值取0.001～0.04 s?1時，2～9 號采樣點As 的模擬值與實測值的配對T檢驗結(jié)果顯示（圖2（d））：當K值介于0.013～0.027 s?1時，P>α=0.05，模擬值與實測值不存在顯著差異，即當沘江水體2～9 號采樣點間As 的綜合衰減系數(shù)K值介于0.013～0.027 s?1時，公式（2）可以較好地反映沘江水體中2～9 號采樣點間As 的遷移轉(zhuǎn)化規(guī)律；當K=0.018 s?1時，P=0.918，公式（2）的模擬值最接近真實值.因此，沘江水體中2～9 號采樣點間Pb 的綜合衰減系數(shù)K的最優(yōu)值取0.018 s?1.

PASW18.0 軟件對K值取0.001～0.04 s?1時，10～12 號采樣點As 的模擬值與實測值的配對T檢驗結(jié)果顯示（圖2（e））：當K值介于0.001～0.04 s?1時，P>α=0.05，模擬值與實測值不存在顯著差異，實際上，當K=0.023 s?1時，10～12 號采樣點As 模擬值的數(shù)量級已達10?5，已非常接近0，即K>0.023 s?1時，K值的增加對T檢驗結(jié)果已不會形成影響.因此，可將K值限定在0.001～0.023 s?1之間，即當沘江水體10～12 號采樣點中As 的綜合衰減系數(shù)介于0.001～0.023 s?1時，公式（2）可以較好地反映沘江水體中10～12 號采樣點間As 的遷移轉(zhuǎn)化規(guī)律；當K=0.002 s?1時，P=0.859，公式（2）的模擬值最接近真實值.因此，沘江水體中10～12 號采樣點間As 的綜合衰減系數(shù)K的最優(yōu)值取0.002 s?1.

綜上所述，礦區(qū)下游沘江水體中Zn 的綜合衰減系數(shù)K值介于0.013～0.017 s?1時，其最優(yōu)值為0.015 s?1.礦區(qū)下游至果朗村（9 號采樣點）附近江段水體Pb 的綜合衰減系數(shù)K值介于0.003～0.009 s?1時，其最優(yōu)值為0.006 s?1，As 的綜合衰減系數(shù)K值介于0.013～0.027 s?1時，其最優(yōu)值為0.018 s?1.此外，果朗村附近至沘江與瀾滄江交匯江段水體Pb的綜合衰減系數(shù)K值介于0.002～0.008 s?1時，其最優(yōu)值為0.004 s?1，As 的綜合衰減系數(shù)K值介于0.001～0.023 s?1時，其最優(yōu)值為0.002 s?1.

3.4 沘江全程水體中重金屬污染物遷移的沿程模擬在MATLAB 中調(diào)用Function 函數(shù)，并對ρ(x)=各相關(guān)參數(shù)進行定義后，即可得到?jīng)a江全程水體中重金屬污染物沿江遷移的模擬模型，也即礦區(qū)下游沘江水體中Zn 的綜合衰減系數(shù)K=0.015 s?1；礦區(qū)下游至果朗村(9 號采樣點)江段水體Pb 的K=0.006 s?1，As 的K=0.018 s?1；果朗村9 號采樣點至沘江與瀾滄江交匯江段Pb 的K=0.004 s?1，As 的K=0.002 s?1.依此模型及擬定的相應K值，只需在礦區(qū)下游附近河流水體中某處分別測得該點Zn、Pb、As 的質(zhì)量濃度值ρ0，并在模型中輸入該值，模型即可自動模擬該測點下游水體中任一點的Zn、Pb、As 質(zhì)量濃度值，并繪制該點下游河流沿程水體中Zn、Pb、As 質(zhì)量濃度變化曲線.

以本次采樣1 號樣點作為ρ0，則礦區(qū)下游ρ0點以下沿程Zn、Pb、As 的模擬值與本次采樣點實測值的對比結(jié)果見圖3.Zn 的沿程模擬值與本次采樣點的實測值極為相近；而Pb 的沿程模擬值與本次采樣點實測值的變化趨勢亦較為接近.雖然個別樣點的實測值與模擬值略有偏差，但下游波動幅度較大的10 號采樣點與沿程模擬值（ρ0點下游88 km處）尚且較為吻合，可見Pb 的沿程模擬值基本能客觀反映水體中Pb 的質(zhì)量濃度變化，模擬效果較好；除個別采樣點外，礦區(qū)下游ρ0點以下As 的沿程模擬值與對應的實測值也基本吻合，即便在下游出現(xiàn)實測值顯著波動的10～12 號采樣點也能較好反映其變化趨勢，模擬計算能夠較好地擬合水體中As的質(zhì)量濃度沿程變化.由此，經(jīng)本文率定綜合衰減系數(shù)K值后的沘江水體重金屬遷移模型，能夠較好地模擬沘江枯水期水體重金屬Zn、Pb、As 的質(zhì)量濃度變化，該模型可用于對沘江枯水期水體重金屬Zn、Pb、As 質(zhì)量濃度值的模擬.

圖3 沿程模擬值與金頂?shù)V區(qū)下游采樣點實測值對比Fig.3 Comparison of simulated values along the way with measured values at downstream sampling points in the Jinding Mining area

4 討論

沘江受蘭坪金頂鉛鋅礦（鉛鋅儲量居世界第二位、亞洲第一位）礦區(qū)影響，水體重金屬Zn、Pb、As 污染嚴重[10]，而沘江流域的生態(tài)健康是維系瀾滄江上游區(qū)域生態(tài)安全的重要保障.因此，摸清沘江水體重金屬Zn、Pb、As 沿程遷移擴散規(guī)律，建立適合沘江水體Zn、Pb、As 遷移轉(zhuǎn)化的模擬模型對流域生態(tài)安全等的研究至關(guān)重要，建模后準確率定模型中各重金屬的綜合衰減系數(shù)，則是模型能否準確模擬水體中重金屬遷移擴散的關(guān)鍵所在.

何紹福等[22]在2015 年研究沙溪沙縣段水環(huán)境容量時，建立了該河段的一維水質(zhì)模型，并采用試錯法率定了CODCr的耗氧速率、氨氮的硝化速率、無機磷的沉積速率等水質(zhì)因子的衰減速率，進而采用試錯法確定了沙溪沙縣段水體中CODCr、氨氮、總磷等的環(huán)境容量，取得了較好的效果.徐凌云等[23]在2017 年研究溫嶺市平原河網(wǎng)水環(huán)境容量時，對溫嶺市平原河網(wǎng)水量及水質(zhì)進行了一維建模，并采用試錯法率定了CODMn、氨氮等的衰減系數(shù)，模擬值與監(jiān)測值擬合較好.本文同樣采用試錯法率定了沘江流域水體重金屬Zn、Pb、As 沿程遷移擴散的綜合衰減系數(shù)，率定結(jié)果也較為理想.因此，在進行水體重金屬遷移擴散的綜合衰減系數(shù)率定時，若實驗條件不允許，則采用試錯法確定水體重金屬遷移擴散的綜合衰減系數(shù)也是一個較為穩(wěn)妥的方法，但也要注意試錯法所采取的步長選擇難免會對K值的準確性產(chǎn)生一定的影響.

由于沘江為天然淺窄的河流，因此本次沘江水體重金屬遷移擴散建模采用的是一維模型，即只考慮縱向（x）重金屬的微量擴散，并未考慮重金屬在橫向（y）及垂向（z）的微量擴散，這也會在一定程度上增加研究結(jié)果的不確定性.本文未考慮河流沿岸工農(nóng)業(yè)生產(chǎn)對沘江水體造成的影響，有待于今后進一步的研究.

此外，沘江豐水期和枯水期的水體流速等水文參數(shù)存在一定的差異，重金屬在水體?懸浮物?底泥間的吸附、解吸、沉降、再懸浮過程及重金屬污染物的化學轉(zhuǎn)化過程亦會不同.若能夠在豐水期對沘江水體再次采樣，并與本次采樣數(shù)據(jù)對比分析，模型將會取得更好的適用性.因此，在考慮對沘江豐水期重金屬Zn、Pb、As 遷移轉(zhuǎn)化進行模擬時，應在豐水期再次采樣，并對模型中Zn、Pb、As 的綜合衰減系數(shù)K分別在其取值范圍內(nèi)重新率定，或取個別樣點進行適用性驗證后，再應用本模型進行模擬.

5 結(jié)論

根據(jù)水體中重金屬沿程遷移的普適模型，結(jié)合沘江水體重金屬的實測數(shù)據(jù)，本文通過試錯法率定了該數(shù)學模型的綜合衰減系數(shù)K的取值范圍及其最優(yōu)值，確定了沘江水體中重金屬遷移的數(shù)學模型，并運用MATLAB 建立了適合沘江全程的水體重金屬遷移轉(zhuǎn)化的模擬模型，得到如下結(jié)論.

（1）枯水期沘江流域金頂?shù)V區(qū)下游Zn 的綜合衰減系數(shù)K值介于0.013～0.017 s?1時，其最優(yōu)值為0.015 s?1；礦區(qū)下游至果朗村（2～9 號采樣點）附近江段水體Pb 的綜合衰減系數(shù)K值介于0.003～0.009 s?1時，其最優(yōu)值為0.006 s?1，As 的綜合衰減系 數(shù)K值介于0.013～0.027 s?1時，其最優(yōu)值為0.018 s?1；果朗村附近至沘江與瀾滄江交匯江段（10～12 號采樣點）水體Pb 的綜合衰減系數(shù)K值介于0.002～0.008 s?1時，其最優(yōu)值為0.004 s?1，As的綜合衰減系數(shù)K值介于0.001～0.023 s?1時，其最優(yōu)值為0.002 s?1.即枯水期沘江流域金頂?shù)V區(qū)下游重金屬Zn、Pb、As 的綜合衰減系數(shù)K分別為Zn 全程0.015 s?1；果朗村上游Pb 為0.006 s?1，As為0.018 s?1；下游Pb 為0.004 s?1，As 為0.002 s?1.

（3）受限于實測數(shù)據(jù)采樣時段的影響，針對沘江豐水期水體重金屬Zn、Pb、As 質(zhì)量濃度沿程遷移變化的模擬，還須對綜合衰減系數(shù)K值重新進行率定或補充個別樣點采樣數(shù)據(jù)進行驗證，以進一步優(yōu)化模擬效果.