段 揚(yáng)，胡伏生，王 旭，雷曉輝

(1.中國(guó)地質(zhì)大學(xué)(北京)水資源與環(huán)境學(xué)院，北京 100083； 2.中國(guó)水利水電科學(xué)研究院，北京 100038)

0 引 言

大型湖泊及水庫(kù)水質(zhì)模擬是研究其生態(tài)演變，污染物運(yùn)移，富營(yíng)養(yǎng)化等問(wèn)題的基礎(chǔ)，其模擬的關(guān)鍵是參數(shù)的不確定性、模擬的準(zhǔn)確性以及參數(shù)之間的影響程度等等一系列問(wèn)題，使水動(dòng)力模型的率定帶來(lái)了很大的難題，那么如何簡(jiǎn)單、快速地對(duì)模型各參數(shù)的敏感性進(jìn)行評(píng)價(jià)也就成為重中之重。伴隨著諸如光合作用、基礎(chǔ)代謝作用、水解作用、光解作用等一系列生化反應(yīng)的進(jìn)行，其模擬所需的參數(shù)眾多，前人在這些方面做了大量的工作[1]，尤其是在多參數(shù)體系下的反演算法。截至目前，參數(shù)之間敏感性度量的主要方法為多次單因素試驗(yàn)法以及正交試驗(yàn)法。前者方法簡(jiǎn)單，操作性強(qiáng)，但其具有較強(qiáng)局限性，尤其是在試驗(yàn)因素及因素水平較多的情況下常常會(huì)得出錯(cuò)誤的結(jié)論，使結(jié)果失真；后者是在全面試驗(yàn)基礎(chǔ)上挑選一些代表性點(diǎn)進(jìn)行試驗(yàn)，使其可以合理反映試驗(yàn)范圍內(nèi)各因素與指標(biāo)間的對(duì)應(yīng)關(guān)系，但其局限性在于當(dāng)因素水平較高時(shí)試驗(yàn)次數(shù)顯著增加，可操作性降低。而我國(guó)科學(xué)家方開(kāi)泰及王元于1978年所提出[2,3]的均勻設(shè)計(jì)法，可以較好地實(shí)現(xiàn)在多因素、多水平的試驗(yàn)中通過(guò)較少試驗(yàn)次數(shù)達(dá)到最有效的選擇試驗(yàn)點(diǎn)并得到關(guān)于系統(tǒng)盡可能充分信息的目標(biāo)。

基于此，通過(guò)大量文獻(xiàn)考究和參數(shù)權(quán)重分析，擬采用均勻設(shè)計(jì)方法進(jìn)行試驗(yàn)設(shè)計(jì)，以丹江口水庫(kù)作為研究對(duì)象，利用全回歸分析手段對(duì)參數(shù)進(jìn)行敏感性分析，確定了各參數(shù)對(duì)模型不確定性的影響大小，并為模型率定提供了基礎(chǔ)，將有利于后期水質(zhì)模擬分析。

1 研究方法

本文選取了美國(guó)William and Mary大學(xué)的John Hamrick所開(kāi)發(fā)的綜合性水質(zhì)模型EFDC建立丹江口水庫(kù)的三維水動(dòng)力學(xué)及水質(zhì)模型，其可以在河流、水庫(kù)、海洋等多種水體進(jìn)行一二三維的水動(dòng)力學(xué)水質(zhì)模擬，在國(guó)內(nèi)外獲得了廣泛應(yīng)用[1,4-9]。其在水動(dòng)力學(xué)方程中滿足Boussinesq近似、靜水壓近似以及準(zhǔn)3D近似等基本假設(shè)條件；在水質(zhì)模塊中，選用質(zhì)量守恒方程作為其控制方程模擬水力輸運(yùn)作用、吸附解析作用、藻類吸收及化學(xué)反應(yīng)作等影響水體水質(zhì)變化過(guò)程。公式(1)為水質(zhì)模塊控制方程：

(1)

式中：C代表水質(zhì)狀態(tài)變量濃度；u，v，w代表x，y，z方向速度分量；Kx，Ky，Kz代表x，y，z方向的湍流擴(kuò)散系數(shù)；Sc代表每單位體積內(nèi)部與外部的源和匯。

2 模型率定

進(jìn)行參數(shù)敏感性分析的前提需建立起準(zhǔn)確的水動(dòng)力學(xué)及水質(zhì)模型。本文經(jīng)過(guò)不斷試算最終確定選用500 m×500 m矩形網(wǎng)格，分為142行，79列，共計(jì)2155個(gè)。垂向上分為5層。綜合考慮丹江口水庫(kù)的入出流情況，選取較大的河流作為出入庫(kù)邊界，將附近的小河并入其相鄰的主河道中，最終包含有3個(gè)入流邊界、1個(gè)出流邊界及1個(gè)水位邊界。同時(shí)設(shè)置7個(gè)監(jiān)測(cè)點(diǎn)來(lái)進(jìn)行實(shí)測(cè)與模擬數(shù)據(jù)對(duì)比。

模擬所需氣象數(shù)據(jù)選取水庫(kù)周圍3個(gè)氣象站點(diǎn)加權(quán)平均值。模擬時(shí)間從2012年3月27日-2012年12月31日，時(shí)間步長(zhǎng)為20 s，初始流速場(chǎng)為0，初始水位場(chǎng)為145 m，水質(zhì)初始場(chǎng)由分布在庫(kù)區(qū)內(nèi)的7個(gè)監(jiān)測(cè)站點(diǎn)數(shù)據(jù)經(jīng)插值得到。

模型率定考慮了丹江口水庫(kù)的水位、水溫及水質(zhì)的實(shí)測(cè)結(jié)果，模擬了其在水動(dòng)力模型基礎(chǔ)上的動(dòng)態(tài)結(jié)果(圖1、圖2及圖3)，發(fā)現(xiàn)模擬結(jié)果與實(shí)測(cè)值擬合度較好，其中水位誤差在0.1%左右，7個(gè)監(jiān)測(cè)站點(diǎn)水溫平均誤差為5.65%，氨氮平均誤差為16.44%，可以較準(zhǔn)確地反映丹江口水庫(kù)的水位、溫度及水質(zhì)在一年中的動(dòng)態(tài)變化趨勢(shì)，可以作為后續(xù)參數(shù)敏感性分析的基礎(chǔ)模型。

圖1 丹江口大壩2012年模擬實(shí)測(cè)日水位過(guò)程圖Fig.1 Water level process map of simulated and Measuredvaluesin front of the Danjiangkou Dam

圖2 涼水河站水溫模擬實(shí)測(cè)值對(duì)比圖Fig.2 Comparison between simulated and measured water temperature values in Liangshuihe Point

圖3 涼水河站氨氮模擬實(shí)測(cè)值對(duì)比圖Fig.3 Comparison between simulated and measured ammonia nitrogen values in Liangshuihe Point

3 參數(shù)敏感性分析

3.1 理論基礎(chǔ)

選取氨氮作為參數(shù)敏感性分析的目標(biāo)，其主要源項(xiàng)包括：①藻類基礎(chǔ)代謝、捕食；②溶解性有機(jī)氮的礦化作用；③沉積床和水柱界面的交換。主要匯項(xiàng)包括：①藻類吸收；②硝化作用。氨氮含量的計(jì)算方程見(jiàn)下式(2)：

(2)

式中：NH4為水中氨氮的含量；FNIx為藻類群x(包含藍(lán)藻、綠藻、硅藻，本次模擬僅考慮藍(lán)藻，后同)代謝的氮中所生成的無(wú)機(jī)氮部分；BMx為藻類基礎(chǔ)新陳代謝速率；FNIP表示被捕食的氮中所生成的溶解性有機(jī)氮部分；PRx為藻類捕食速率；PNx為藻類銨吸收偏好；Px為藻類生產(chǎn)速率；ANCx表示藻類氮對(duì)碳比例；Bx為藻類藻類生物量；KDON表示溶解性有機(jī)氮礦化速率；Nit為硝化速率；BFNH4為沉積物-水的銨交換通量；WNH4表示銨的外部負(fù)荷。

3.2 評(píng)價(jià)參數(shù)體系

通過(guò)文獻(xiàn)調(diào)研[10-14]，選擇藻類基礎(chǔ)新陳代謝速率(BMR)、藻類捕食速率(PRR)、藻類生長(zhǎng)吸收氮半飽和常數(shù)(KHN)、藻類生長(zhǎng)速率(PM)、溶解態(tài)有機(jī)氮水解速率(KDN)、硝化作用半飽和常數(shù)(KHNITDO)、最大硝化速率(NITM)等7個(gè)參數(shù)作為敏感性分析的輸入?yún)?shù)。假定參數(shù)間相互獨(dú)立，根據(jù)已有文獻(xiàn)[11-15]確定參數(shù)最大值與最小值范圍(見(jiàn)表1)。

表1 參數(shù)敏感性分析模型輸入?yún)?shù)的范圍Tab.1 Ranges of parameters used for the parameter sensitivity analysis model

本次采用的均勻設(shè)計(jì)試驗(yàn)方法，其優(yōu)點(diǎn)在于每個(gè)均勻設(shè)計(jì)表都會(huì)有一個(gè)使用表，可指導(dǎo)我們?nèi)绾芜x擇均勻設(shè)計(jì)表中的某些列來(lái)進(jìn)行試驗(yàn)設(shè)計(jì)以及由這些列所組成試驗(yàn)方案的均勻度，而均勻度的好壞一般用均勻度的偏差D來(lái)衡量，其值越小，表明均勻度越好。

3.3 參數(shù)變換

通過(guò)查閱均勻設(shè)計(jì)表[12]篩選出符合條件的所有表，對(duì)比當(dāng)S=7時(shí)的均勻度偏差D，選出當(dāng)D最小時(shí)的方案，最終確定為U*28(288)，D=0.155 0。考慮到所有參數(shù)取值范圍中中最大值與最小值相差達(dá)340倍，直接做回歸分析效果欠佳，需首先對(duì)各參數(shù)進(jìn)行對(duì)數(shù)變換[13]，變換后具體試驗(yàn)參數(shù)列表見(jiàn)表2。

表2 均勻變換參數(shù)表Tab.2 Parameters for uniform design

將表2中參數(shù)分別帶入EFDC模型中，運(yùn)算28次，得到28組特定目標(biāo)的結(jié)果。將7個(gè)監(jiān)測(cè)點(diǎn)的氨氮濃度序列輸出，并選取4月1日，7月15日，10月1日，12月1日4個(gè)時(shí)間點(diǎn)來(lái)反映氨氮濃度一年內(nèi)變化情況。之后將7個(gè)監(jiān)測(cè)點(diǎn)在4個(gè)具體時(shí)間的數(shù)值與其對(duì)應(yīng)參數(shù)耦合后分別進(jìn)行回歸分析，研究這7個(gè)參數(shù)對(duì)于模型模擬結(jié)果的不確定性影響。

4 結(jié)果與討論

4.1 不確定性影響

以7個(gè)監(jiān)測(cè)點(diǎn)分別在4個(gè)具體時(shí)段的氨氮含量作為輸出值，對(duì)28種參數(shù)變化環(huán)境下所得結(jié)果進(jìn)行不確定性分析。以何家灣站點(diǎn)4月1日計(jì)算結(jié)果為例，在28次試驗(yàn)過(guò)程中其氨氮含量最大值為0.122 mg/L，最小值僅為0.005 8 mg/L，二者相差21倍，表明由于參數(shù)不同對(duì)于氨氮含量產(chǎn)生較大變化，證明參數(shù)不確定性將直接導(dǎo)致模型結(jié)果不確定性。通過(guò)各點(diǎn)氨氮濃度變化方差可以更明確地反映其不確定性大小，詳見(jiàn)表3。

由于方差越大表征參數(shù)對(duì)于模型計(jì)算不確定性影響越大[15]，在春季各點(diǎn)方差值大致相同，均保持在0.03左右，僅在江北大橋處略?。欢募居捎诟黝惿磻?yīng)進(jìn)行較旺盛，由于參數(shù)不確定所造成的模型結(jié)果差異程度被放大，方差范圍為0.004～0.062，且?guī)靺^(qū)北部由于受到漢江及丹江入流影響其不確定性增大；秋季同夏季情況類似，但由于入流量降低且主要入流量來(lái)自于漢江使得西北部肖川站點(diǎn)不確定性達(dá)到最大；而冬季隨著溫度及日照減少，生化反應(yīng)速率降低，模型趨于穩(wěn)定，各點(diǎn)方差相差不大，僅有涼水河相較其他站點(diǎn)略大。

表3 各觀測(cè)站點(diǎn)不同時(shí)段氨氮濃度方差匯總表Tab.3 Variance of ammonia nitrogen in each monitoring point

4.2 參數(shù)敏感性分析

將7個(gè)樣本點(diǎn)的計(jì)算進(jìn)行回歸分析，并利用方差檢驗(yàn)來(lái)判斷試驗(yàn)結(jié)果的可靠性。首先，將之前設(shè)置的7項(xiàng)敏感性參數(shù)候選項(xiàng)作為自變量，將氨氮濃度作為因變量帶入回歸分析軟件中進(jìn)行自動(dòng)建模，可以初步對(duì)于自變量相對(duì)重要性進(jìn)行判斷。以陶岔站7月15日計(jì)算結(jié)果為例，經(jīng)過(guò)自動(dòng)線性建模得到其準(zhǔn)確度為90.7%，之后得出各變量相關(guān)重要性圖，即各個(gè)參數(shù)對(duì)于模型不確定性的貢獻(xiàn)率，見(jiàn)圖4。

圖4 7個(gè)參數(shù)對(duì)于氨氮含量不確定性貢獻(xiàn)圖Fig.4 Level of contributions to the uncertainty for the 7 parameters in ammonia nitrogen

可以看出：硝化作用半飽和常數(shù)對(duì)于氨氮含量不確定性最大，達(dá)到39%；其次是藻類生長(zhǎng)吸收氮半飽和常數(shù)以及藻類生長(zhǎng)速率，分別為21%和19%；前3項(xiàng)占到總貢獻(xiàn)度的80%，后面四項(xiàng)參數(shù)貢獻(xiàn)度較低。

由于自動(dòng)建模無(wú)法求得模型方程及每項(xiàng)參數(shù)的回歸系數(shù)，所以之后將進(jìn)行全回歸分析。由于篇幅限制，仍以陶岔站7月15日結(jié)果為例，取顯著性水平α=0.05，樣本數(shù)量為28，檢驗(yàn)值F=23.245，通過(guò)查表查得F7,20(0.05)=2.514，可以看出F>F7,20(0.05)，表明氨氮的變化確實(shí)由于上述7個(gè)參數(shù)的變化所引起，回歸方程顯著。另外通過(guò)查驗(yàn)7個(gè)參數(shù)的相關(guān)系數(shù)發(fā)現(xiàn)，其中任意2個(gè)參數(shù)相關(guān)系數(shù)最大值為0.3，相關(guān)性較低，驗(yàn)證了之前關(guān)于7個(gè)參數(shù)相互獨(dú)立的假設(shè)。

經(jīng)過(guò)回歸分析可以得到各個(gè)參數(shù)的標(biāo)準(zhǔn)回歸系數(shù)并建立回歸方程：

NH4=0.057+0.125×BMR-0.131×PRR+

0.483×KHN+0.347×PM+0.319KDN+

0.645KHNITDO-0.232×NITM

(3)

其中，藻類捕食速率(PRR)及最大硝化速率(NITM)與氨氮濃度均呈負(fù)相關(guān)關(guān)系。

將各參數(shù)回歸系數(shù)及t檢驗(yàn)絕對(duì)值進(jìn)行排序可得出參數(shù)敏感性大小順序?yàn)椋合趸饔冒腼柡统?shù)(KHNITDO)>藻類生長(zhǎng)吸收氮半飽和常數(shù)(KHN)>藻類生長(zhǎng)速率(PM)>溶解態(tài)有機(jī)氮水解速率(KDN)>最大硝化速率(NITM)>藻類捕食速率(PRR)>藻類基礎(chǔ)新陳代謝速率(BMR)。

順序與圖5中對(duì)不確定性貢獻(xiàn)大小順序相同。對(duì)比同一監(jiān)測(cè)點(diǎn)4個(gè)時(shí)段的回歸系數(shù)結(jié)果可發(fā)現(xiàn)在一年的不同時(shí)段，參數(shù)敏感程度略有不同。以肖川站為例，表4顯示的是其在各時(shí)段內(nèi)各參數(shù)回歸系數(shù)大小?？梢钥闯觯趸饔冒腼柡统?shù)在4個(gè)時(shí)段回歸系數(shù)均位于首位且遠(yuǎn)大于后面各參數(shù)，然而之后排位與季節(jié)有較大關(guān)系，春季由于水中氨氮含量不高使得硝化作用進(jìn)行程度較低，由于經(jīng)過(guò)冬季有機(jī)殘?jiān)姆纸饧霸孱愃劳鏊尫盘幍娜芙庑杂袡C(jī)氮的礦化水解作用得到增強(qiáng)；而夏季由于氣溫與日照較強(qiáng)，藻類生長(zhǎng)茂盛，其生長(zhǎng)過(guò)程所吸收氮成為影響氨氮濃度的一個(gè)重要因素，且同時(shí)硝化作用作為好氧反應(yīng)受到夏季溶解氧含量降低影響而減弱；到了秋季，上游地區(qū)存在許多未被作物吸收的氮肥隨降雨徑流進(jìn)入庫(kù)區(qū)，且此時(shí)入流量的快速增加且其中含有大量氨氮，氨氮數(shù)量充足，使得硝化作用進(jìn)行充分，此時(shí)最大硝化速率將在很大程度上影響氨氮濃度變化。其他站點(diǎn)結(jié)果類似，僅在少數(shù)排位上存在偏差，但并沒(méi)有發(fā)生較大地區(qū)差異。

表4 肖川站各時(shí)段參數(shù)回歸系數(shù)表Tab.4 Regression coefficients in different period in XiaoChuan

5 結(jié) 語(yǔ)

(1)根據(jù)丹江口水庫(kù)的底部高程、氣象、水文、環(huán)境實(shí)測(cè)數(shù)據(jù)，基于EFDC模型，建立三維水動(dòng)力水質(zhì)模型，并將模擬結(jié)果與實(shí)測(cè)數(shù)據(jù)進(jìn)行對(duì)比取得了良好的效果，表明模型具有很強(qiáng)的精確度。

(2)在已率定模型基礎(chǔ)上，以氨氮為目標(biāo)，選擇藻類基礎(chǔ)新陳代謝速率等7個(gè)參數(shù)進(jìn)行參數(shù)敏感度分析。利用均勻設(shè)計(jì)理論設(shè)計(jì)試驗(yàn)，可以在較少試驗(yàn)次數(shù)的情況下充分利用模型信息，快速準(zhǔn)確地實(shí)現(xiàn)參數(shù)敏感度大小判斷。

(3)通過(guò)對(duì)計(jì)算結(jié)果進(jìn)行不確定性分析可知，參數(shù)改變對(duì)于氨氮含量影響較大，其最大最小值相差可達(dá)20余倍。

(4)經(jīng)過(guò)對(duì)7個(gè)監(jiān)測(cè)站點(diǎn)4個(gè)時(shí)段結(jié)果的比較，發(fā)現(xiàn)硝化作用半飽和常數(shù)敏感度最大，其次是溶解態(tài)有機(jī)氮水解速率、最大硝化速率、藻類生長(zhǎng)吸收氮半飽和常數(shù)，而藻類生長(zhǎng)速率、藻類基礎(chǔ)新陳代謝速率、藻類捕食速率的敏感度較低。故今后對(duì)于其他大型水庫(kù)水質(zhì)模型率定時(shí)可對(duì)敏感度較低的參數(shù)采用經(jīng)驗(yàn)參數(shù)，提高率定工作效率。

□

[1] Cerco C F, Cole T M. Three-dimensional eutrophication model of Chesapeake Bay[J]. Journal of Environment Engineering, 1993,119:1 006-1 025.

[2] 方開(kāi)泰.均勻設(shè)計(jì)與均勻設(shè)計(jì)表[M].北京:科學(xué)出版社,1994.

[3] 燕 喬,吳長(zhǎng)彬,張 巖. 基于均勻設(shè)計(jì)的鄧肯E-B模型參數(shù)敏感性分析[J]. 中國(guó)農(nóng)村水利水電，2010,(7):82-85.

[4] Ji Z-G. Review of hydrodynamics and water quality: modeling rivers lakes[J]. Journal of Hydrologic Engineering, 2009,(8):892-893.

[5] Jin K R, Hamrick J H, Tisdale T. Application of three-dimensional hydrodynamic model for Lake Okeechobee[J]. Journal of Hydraulic Engineering-Asce, 2000,126:758-771.

[6] 李云生,劉偉江,吳悅穎,等. 美國(guó)水質(zhì)模型研究進(jìn)展綜述[J]. 水利水電技術(shù),2006,(2):68-73.

[7] Jin K R, J M Harmrick, T S Tisdale. Application of a Three-dimensional hydrodynamic model for Lake Okeechobee[J]. Journal of Hydraulic Engineering, 2000,106(11):758-772.

[8] Hamrick J M, Wm B Mills. Analysis of temperatures in connoting pond as influenced by the peach bottom atomic power plant thermal discharge[J].Environmental Science and Policy,2001(3):197-209.

[9] 任華堂,于 良,夏建新,等. 黃河內(nèi)蒙古段水污染事故應(yīng)急預(yù)警模型研究[J].應(yīng)用基礎(chǔ)與工程科學(xué)學(xué)報(bào)，2012(9).67-75.

[10] WU G, XU Z. Prediction of algal blooming using EFDC model: Case study in the Daoxiang Lake [J]. Ecol Model, 2011,222(6):1 245-52.

[11] WANG Y, JIANG Y, LIAO W, et al. 3-D hydro-environmental simulation of Miyun reservoir, Beijing[J]. Journal of Hydro-environment Research, 2014,8(4):383-95.

[12] 王宇暉. 基于分布式水文模型的流域水循環(huán)及面源伴生過(guò)程模擬研究與應(yīng)用[D]. 上海：東華大學(xué)，2012.

[13] 李 興. 內(nèi)蒙古烏梁素海水質(zhì)動(dòng)態(tài)數(shù)值模擬研究[D]. 呼和浩特：內(nèi)蒙古農(nóng)業(yè)大學(xué)，2009.

[14] 李 杰. Falls Lake水庫(kù)水質(zhì)的數(shù)值模擬研究及分析[D]. 山東青島：中國(guó)海洋大學(xué)，2010.

[15] 李一平，唐春燕，余鐘波,等. 大型淺水湖泊水動(dòng)力模型不確定性和敏感性分析[J]. 水科學(xué)進(jìn)展，2012,23(2):271-277.