李 平，馬 放

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

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

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

1 材料與方法

1.1 河段水文概況

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

1.2 羅丹明B示蹤及水團(tuán)追蹤試驗(yàn)

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

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

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

2 結(jié)果與分析

2.1 BOD-DO模型選擇及改進(jìn)

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

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

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

式中:ρ(CBOD)為河水質(zhì)量濃度，mg·L-1;ρ(NBOD)為河水 NBOD 質(zhì)量濃度，mg·L-1;ρ(O)為河水DO質(zhì)量濃度，mg·L-1;ρ(Os)為河水t℃時(shí)飽和溶解氧質(zhì)量濃度，mg·L-1;x為距排污口河水流動(dòng)距離，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的物理沉浮作用表達(dá)清楚，綜合Thomas，引入K3表征懸浮作用的CBOD物理變化和O'Connor的CBOD總衰減率的改進(jìn)方法，改進(jìn)O'Connor模型，以對(duì)松花江佳木斯江段BOD-DO進(jìn)行精確模擬:

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

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

CBOD方程解為

NBOD方程解為

則假定為一維無(wú)離散，溶解氧方程變?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'為任意斷面時(shí)間.由式(10)-(13)構(gòu)成改進(jìn)O'Connor模型解析解.松花江佳木斯段溶解氧模型由以上方程組建立.

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

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

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

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

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

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

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

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

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

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

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

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

式中:θ為溫度修正系數(shù)，一般介于1.047～1.140之間;K(t)為任意溫度t時(shí)的參數(shù);K(20)為20℃時(shí)的參數(shù).根據(jù)表3的參數(shù)識(shí)別結(jié)果，確定相應(yīng)溫度修正系數(shù)，進(jìn)行溫度修正，可得松花江佳木斯段溶解氧模型參數(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ù)的驗(yàn)證

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

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

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

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

對(duì)改進(jìn)O'Connor模型溶解氧有關(guān)的各項(xiàng)參數(shù)進(jìn)行靈敏度分析.檢驗(yàn)的方式是保持其他參數(shù)不變，即只改變被檢驗(yàn)參數(shù)，令模型中Kd，Kn，Ka，Kr，Km等各項(xiàng)參數(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兩個(gè)斷面時(shí)刻DO變化幅度考察其靈敏性，結(jié)果如圖2所示.

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

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

3 結(jié)論

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

2)對(duì)松花江佳木斯段BOD-DO模型參數(shù)進(jìn)行了任意溫度下的參數(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)對(duì)模型預(yù)測(cè)結(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模型計(jì)算河流水功能區(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ì)模型在流域水資源保護(hù)管理中的應(yīng)用［J］.水利學(xué)報(bào)，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］ 國(guó)家環(huán)境保護(hù)總局.水和廢水監(jiān)測(cè)分析方法［M］.4版.北京:中國(guó)環(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ù)反問(wèn)題偶合模型的研究［J］.水科學(xué)進(jìn)展，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］ 杜彥良，褚君達(dá).河流連續(xù)點(diǎn)源的垂向混合規(guī)律研究及其應(yīng)用［J］.環(huán)境科學(xué)學(xué)報(bào)，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)濟(jì)，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)