李 平，馬 放

(哈爾濱工業(yè)大學(xué)城市水資源與水環(huán)境國家重點實驗室，黑龍江哈爾濱150090)

河流水質(zhì)模擬在進行流域污染控制、水環(huán)境規(guī)劃、水質(zhì)預(yù)測及評價中發(fā)揮著重要的作用，為污染物總量控制和水功能區(qū)水質(zhì)目標環(huán)境管理提供了有效的理論支撐［1-2］.水質(zhì)模擬過程是利用水體歷史數(shù)據(jù)、污染物的物理化學(xué)演變規(guī)律，建立相應(yīng)的數(shù)學(xué)、物理、化學(xué)及試驗?zāi)Ｐ兔枋鏊w中污染物質(zhì)與河流水質(zhì)之間的響應(yīng)關(guān)系，對河流水環(huán)境系統(tǒng)行為進行模擬、預(yù)測和評價，是水環(huán)境監(jiān)督管理的重要工具［3-4］.河流水質(zhì)模擬中模型各項參數(shù)(生化耗氧速率、大氣復(fù)氧速率、生化耗氧量的沉淀和懸浮系數(shù)、氮化合物的硝化速率等)的測定和估計是建立水質(zhì)模型的重要工作前提和基礎(chǔ)，模型參數(shù)是否精確直接關(guān)系水質(zhì)模擬能否真實反應(yīng)污染實際過程［5-6］.

為了進一步準確模擬、預(yù)測和評價松花江流域河流水環(huán)境系統(tǒng)行為，本研究利用O'Connor水質(zhì)模型建立松花江河段BOD-DO之間的響應(yīng)關(guān)系，為進一步開展流域水環(huán)境水質(zhì)模擬研究和松花江流域水環(huán)境保護及污染物總量優(yōu)化分配奠定基礎(chǔ).

1 材料與方法

1.1 河段水文概況

松花江佳木斯段全長201 km，流經(jīng)佳木斯、樺川、綏濱、富錦、同江等市縣，并與黑龍江匯合于同江市.松花江下游地處三江平原，河道寬，分叉較多.歷史上徑流量年變化幅度較大，存在明顯的豐、枯周期.枯水年與豐水年徑流量比為差異較大，年內(nèi)徑流量極不平衡，一年四季存在豐、平、枯3個水期.豐、枯期流量相差百倍以上.試驗江段為松花江下游佳木斯-富錦江段，長167 km，河道彎曲系數(shù)1.17.

1.2 羅丹明B示蹤及水團追蹤試驗

羅丹明B示蹤試驗:淡水、河水、海水中彌散僅從理論上來分析比較困難，采用現(xiàn)場示蹤試驗，選用易溶于水、稀溶液呈強熒光的工業(yè)染料堿性羅丹明B為示蹤劑，采用熒光分光光度法進行測定，確定縱向彌散系數(shù)，研究干流擴散規(guī)律.示蹤劑羅丹明B(分析純)瞬時投放，在河流的下游布置斷面采樣，樣品送回化驗室檢測，波長578 nm定量，水質(zhì)指標采用國家標準方法［7］測定.

水團追蹤試驗:利用人工投放的方式在河道中造成一個含有機物的污水團，進行水團跟蹤監(jiān)測，利用水團質(zhì)量濃度沿程衰變情況，來推定河流對污染物質(zhì)的降解規(guī)律.BOD(20)為現(xiàn)場測定，冰封期、明水期在佳木斯江段前后共進行15次水團追蹤作業(yè).冰封期水溫為0℃，以初始投加為第1天計，共監(jiān)測14 d，各斷面監(jiān)測時刻如下:佳木斯下1∶00，樺川6∶00，新城14∶30，富錦18∶00;明水期為14 ℃，以初始投加為第1天計，各斷面監(jiān)測時刻如下:佳木斯下8∶00，樺川19∶00，新城8∶00，富錦400.表1列出了15次實測數(shù)據(jù).

表1 松花江佳木斯江段水團追蹤數(shù)據(jù)表 mg·L-1

2 結(jié)果與分析

2.1 BOD-DO模型選擇及改進

污染物在河流水動力的作用下發(fā)生稀釋、擴散、沉降、降解等過程.根據(jù)河流的水文特點、污染物性質(zhì)建立BOD-DO模型，來描述水體中污染物隨時間和空間遷移轉(zhuǎn)化規(guī)律的數(shù)學(xué)方程［8］，對這一過程進行預(yù)測.溶解氧質(zhì)量濃度是衡量現(xiàn)有水質(zhì)與水體能否保證水生物系統(tǒng)良性循環(huán)的重要標志.因此建立溶解氧模型描述河流水質(zhì)狀態(tài)，進而建立排污量與水質(zhì)狀態(tài)的定量關(guān)系，則可將河流納污能力與污染排放，用模型的模擬解來表示，構(gòu)成水環(huán)境容量的定量計算方法［9-10］.取6月松花江佳木斯某段水樣，檢測其20 d的BOD質(zhì)量濃度變化過程，如圖1所示.

圖1 松花江佳木斯段20 d的BOD變化過程

由圖1可知:該段BOD質(zhì)量濃度存在2個明顯階段變化，即CBOD降解和NBOD降解，顯示松花江該段存在著硝化過程，且硝化過程較為超前，故而選擇同時考慮CBOD和NBOD的衰減與耗氧作用的O'Connor模型進行松花江佳木斯江段BOD-DO模擬.基本方程如下:

式中:ρ(CBOD)為河水質(zhì)量濃度，mg·L-1;ρ(NBOD)為河水 NBOD 質(zhì)量濃度，mg·L-1;ρ(O)為河水DO質(zhì)量濃度，mg·L-1;ρ(Os)為河水t℃時飽和溶解氧質(zhì)量濃度，mg·L-1;x為距排污口河水流動距離，m;u為河水流速，m·s-1;Kd為CBOD衰減系數(shù)，d-1;Ka為河水復(fù)氧系數(shù)，d-1;K3(d)為CBOD沉浮系數(shù)，d-1;Kn為NBOD衰減系數(shù)，d-1.

但O'Connor基本模型并未將CBOD，NBOD的物理沉浮作用表達清楚，綜合Thomas，引入K3表征懸浮作用的CBOD物理變化和O'Connor的CBOD總衰減率的改進方法，改進O'Connor模型，以對松花江佳木斯江段BOD-DO進行精確模擬:

則式(1)變?yōu)?/p>

式中:Kr為CBOD總衰減速率，d-1;Km為NBOD總衰減速率，d-1;K3(n)為NBOD物理沉浮系數(shù)，d-1.

CBOD方程解為

NBOD方程解為

則假定為一維無離散，溶解氧方程變?yōu)槭?12).為方便式(12)求解，令，根據(jù)常數(shù)變易法，可求出式(12)解析解，即溶解氧方程為(13).公式如下所示:

式中:ρo(CBOD)為初始斷面河水CBOD質(zhì)量濃度，mg·L-1;ρo(NBOD)為初始斷面河水NBOD質(zhì)量濃度，mg·L-1;ρo(OO)為初始斷面溶解氧質(zhì)量濃度，mg·L-1;t'為任意斷面時間.由式(10)-(13)構(gòu)成改進O'Connor模型解析解.松花江佳木斯段溶解氧模型由以上方程組建立.

2.2 BOD-DO模型的參數(shù)識別分析

2.2.1 模型參數(shù)的識別

復(fù)氧系數(shù)Ka對于河流溶解氧質(zhì)量濃度變化影響較大，國內(nèi)一般均以采用經(jīng)驗公式估值最普遍，但由于時空限制，模擬精度較差，不同公式估值相差較大，另外利用示蹤劑實測Ka，因價格昂貴，適用范圍窄，很難廣泛采用.CBOD的衰減系數(shù)Kd、NBOD衰減系數(shù)Kn識別是通過利用BOD的20 d過程線進行試驗室模擬，完成相應(yīng)參數(shù)識別.根據(jù)松花江佳木斯下游斷面BOD的20 d實測過程數(shù)據(jù)，采用最小二乘法既可確定目標Kd和Kn，其中Kn一般在0.1～0.3 d-1.結(jié)合試驗前期資料獲得了耗氧系數(shù)常數(shù)Kd(20℃)=0.250 d-1，Kn(20℃)=0.172 d-1，簡化了松花江佳木斯段水質(zhì)模型參數(shù)Ka，Kr和Km識別過程，可以實現(xiàn)利用水團追蹤試驗中溶解氧實測數(shù)據(jù)結(jié)合多維參數(shù)的最優(yōu)化估值法識別Ka，Kr和Km.

多維參數(shù)最優(yōu)基本方法是從給定的初始參數(shù)值Ki出發(fā)，每次增減一定的量，逐步改善目標函數(shù)J(K)，直到其滿足目標值收斂的誤差要求.研究采用梯度搜索法，在起點的目標函數(shù)值J(Ki)下降速率最快的方向(一階負梯度方向)，按定步長進行搜索，每次改進目標函數(shù)值，得到新的起點，反復(fù)迭代計算，直到滿足要求.目標函數(shù)取溶解氧的模型計算值與實測值之差的平方和為

用梯度搜索法尋出最優(yōu)解.

式中:ρij(O)為模型計算的第i次監(jiān)測第j個斷面的河水溶解氧質(zhì)量濃度，mg·L-1;ρij(Om)為實測得到的第i次監(jiān)測第j斷面的河水溶解氧質(zhì)量濃度，mg·L-1.

改進的O'Connor模型參數(shù)估值采用最優(yōu)搜索計算框圖［11-12］，參數(shù)識別分冰封期和明水期分別進行，獲得15組實測數(shù)據(jù)，10組用于識別，5組用于參數(shù)驗證.1-10號為冰封期實測數(shù)據(jù)選擇1-7號數(shù)進行參數(shù)識別，8-10號數(shù)據(jù)進行參數(shù)驗證;11-15號為明水期實測數(shù)據(jù)選擇11，13和14號數(shù)據(jù)進行參數(shù)識別，12和15號數(shù)據(jù)進行參數(shù)驗證.利用表1水團追蹤數(shù)據(jù)以及Kd=0.250 d-1，Kn=0.172 d-1輸入值推求冰封期和明水期的Ka，Kr，Km值，相應(yīng)參數(shù)識別結(jié)果如表2所示.

表2 參數(shù)識別結(jié)果 d-1

分別將冰封期、明水期參數(shù)識別值取平均值，作為冰封期和明水期參數(shù)優(yōu)選值，如表3所示.

表3 冰封期和明水期參數(shù)值 d-1

2.2.2 模型參數(shù)的溫度修正

模型參數(shù)修正方程基本型式為

式中:θ為溫度修正系數(shù)，一般介于1.047～1.140之間;K(t)為任意溫度t時的參數(shù);K(20)為20℃時的參數(shù).根據(jù)表3的參數(shù)識別結(jié)果，確定相應(yīng)溫度修正系數(shù)，進行溫度修正，可得松花江佳木斯段溶解氧模型參數(shù)修正方程如下:Kd=0.25×1.047t-20;Kn=0.172×1.088t-20;Kr=0.833×1.011t-20;Km=0.651×1.025t-20;Ka=0.87×1.024t-20.

2.2.3 模型參數(shù)的驗證

分別利用冰封期8-10號數(shù)和明水期12，15號數(shù)據(jù)進行參數(shù)驗證.驗證結(jié)果如表4所示.

表4驗證結(jié)果可以表明，分不同水期識別并驗證得到的模型參數(shù)值，在精度方面可以滿足水質(zhì)模擬的要求，利用改進的模型推求得模擬計算值與實測值的相對誤差，除兩點由于實際溶解氧ρ(DO)較低而偏大外，分別為14.79%，50.69%，其余均在10%之內(nèi).15個模擬計算值與實測值的相對誤差平均在10%以下，所確定的模型參數(shù)應(yīng)用在研究松花江佳木斯段水質(zhì)模擬中較為合適.

表4 松花江溶解氧模型驗證結(jié)果

2.3 參數(shù)靈敏度分析

對改進O'Connor模型溶解氧有關(guān)的各項參數(shù)進行靈敏度分析.檢驗的方式是保持其他參數(shù)不變，即只改變被檢驗參數(shù)，令模型中Kd，Kn，Ka，Kr，Km等各項參數(shù)分別增減幅度10%，ρo(CBOD)=6.50 mg·L-1，ρo(NBOD)=11.55 mg·L-1，ρo(O)=8.00 mg·L-1，由模型輸出12∶00和24∶00兩個斷面時刻DO變化幅度考察其靈敏性，結(jié)果如圖2所示.

圖2 模型參數(shù)靈敏度分析

由圖2模型參數(shù)變化分析可以判定模型對應(yīng)靈敏程度，松花江佳木斯段改進的O'Connor溶解氧模型中，大氣復(fù)氧系數(shù)Ka和Kn變化所引起的DO變化幅度最大，CBOD耗氧系數(shù)Kd次之，BOD總衰減速率Km和Kr靈敏度最低.

3 結(jié)論

1)應(yīng)用改進的O'Connor溶解模型研究松花江佳木斯段BOD-DO變化關(guān)系，水團追蹤試驗驗證了模型耗氧系數(shù)、復(fù)氧系數(shù)及總衰減系數(shù)等重要參數(shù)取值的合理性，其精密可以滿足水質(zhì)模擬的要求.

2)對松花江佳木斯段BOD-DO模型參數(shù)進行了任意溫度下的參數(shù)修正，Kd=0.25×1.047t-20，Kn=0.172 ×1.088t-20，Kr=0.833 ×1.011t-20，Km=0.651×1.025t-20，Ka=0.87×1.024t-20.

3)對模型預(yù)測結(jié)果影響較大的參數(shù)為大氣復(fù)氧系數(shù)Ka和Kn，CBOD耗氧系數(shù)Kd次之，BOD總衰減速率Km和Kr靈敏度最低.

References)

［1］ Sullivan A B，Jager H I，Myers R.Modeling white sturgeon movement in a reservoir:the effect of water quality and sturgeon density ［J］.Ecological Modelling，2003(167):97-114.

［2］ Madsen H.Parameter estimation in distributed hydrological catchment modelling using automatic calibration with multiple objectives［J］.Advances in Water Resources，2003，26:205-216.

［3］ 劉 偉，劉洪超，徐海巖.基于MIKE11模型計算河流水功能區(qū)納污能力方法［J］.東北水利水電，2009(8):69-72.Liu Wei，Liu Hongchao，Xu Haiyan.Calculation method of water environment capacity for water function area based on MIKE 11model［J］.Water Resources＆Hydropower of Northeast China，2009(8):69-72.(in Chinese)

［4］ 金春久，王 超，范曉娜，等.松花江干流水質(zhì)模型在流域水資源保護管理中的應(yīng)用［J］.水利學(xué)報，2010，41(1):86-92.Jin Chunjiu，Wang Chao，F(xiàn)an Xiaona，et al.Application of water quality model to the river basin water resources protection and management for the mainstream of the Songhua River［J］.Journal of Hydraulic Engineering，2010，41(1):86-92.(in Chinese)

［5］ 張永祥，蔣 源，施同平，等.基于WASP 7.2河流水質(zhì)模型的應(yīng)用研究［J］.北京水務(wù)，2010(1):22-24.Zhang Yongxiang，Jiang Yuan，Shi Tongping，et al.Study on the application of model about river water quality［J］.Beijing Water，2010(1):22-24.(in Chinese)

［6］ 宋新山，鄧 偉.環(huán)境數(shù)學(xué)模型［M］.北京:科學(xué)出版社，2004.

［7］ 國家環(huán)境保護總局.水和廢水監(jiān)測分析方法［M］.4版.北京:中國環(huán)境科學(xué)出版社，2002.

［8］ McAvoy D C，Masscheleyn P，Peng C，et al.Risk assessment approach for untreated wastewater using the QUAL2E water quality model［J］.Chemosphere，2003(52):55-66.

［9］ 李 蘭，李志永，劉金才.BOD5-DO參數(shù)反問題偶合模型的研究［J］.水科學(xué)進展，2000，11(3):255-259.Li Lan，Li Zhiyong，Liu Jincai.A BOD5-DO parameter inverse coincidence model［J］.Advances in Water Science，2000，11(3):255-259.(in Chinese)

［10］ 杜彥良，褚君達.河流連續(xù)點源的垂向混合規(guī)律研究及其應(yīng)用［J］.環(huán)境科學(xué)學(xué)報，2000，20(6):699-703.Du Yanliang，Chu Junda.Study on verticalmixing law of continuous point source in nature river and its application ［J］.ACTA Scientiae Circumstantiae，2000，20(6):699-703.(in Chinese)

［11］ 芮孝芳.河流水文學(xué)［M］.南京:河海大學(xué)出版社，2003:143-156.

［12］ 馬常仁，徐得潛，周 慧，等.基于最小二乘法與最速下降法的綜合衰減系數(shù)的率定［J］.水利科技與經(jīng)濟，2012，18(1):20-24.Ma Changren，Xu Deqian，Zhou Hui，et al.Rating at synthetic attenuation coefficient based on least-square method and the steepest descent method ［J］.Water Conservancy Science and Technology and Economy，2012，18(1):20-24.(in Chinese)