段 揚，胡伏生，王 旭，雷曉輝

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

0 引 言

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

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

1 研究方法

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

(1)

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

2 模型率定

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

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

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

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

圖2 涼水河站水溫模擬實測值對比圖Fig.2 Comparison between simulated and measured water temperature values in Liangshuihe Point

圖3 涼水河站氨氮模擬實測值對比圖Fig.3 Comparison between simulated and measured ammonia nitrogen values in Liangshuihe Point

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

3.1 理論基礎(chǔ)

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

(2)

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

3.2 評價參數(shù)體系

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

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

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

3.3 參數(shù)變換

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

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

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

4 結(jié)果與討論

4.1 不確定性影響

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

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

表3 各觀測站點不同時段氨氮濃度方差匯總表Tab.3 Variance of ammonia nitrogen in each monitoring point

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

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

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

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

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

經(jīng)過回歸分析可以得到各個參數(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檢驗絕對值進行排序可得出參數(shù)敏感性大小順序為：硝化作用半飽和常數(shù)(KHNITDO)>藻類生長吸收氮半飽和常數(shù)(KHN)>藻類生長速率(PM)>溶解態(tài)有機氮水解速率(KDN)>最大硝化速率(NITM)>藻類捕食速率(PRR)>藻類基礎(chǔ)新陳代謝速率(BMR)。

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

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

5 結(jié) 語

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

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

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

(4)經(jīng)過對7個監(jiān)測站點4個時段結(jié)果的比較，發(fā)現(xiàn)硝化作用半飽和常數(shù)敏感度最大，其次是溶解態(tài)有機氮水解速率、最大硝化速率、藻類生長吸收氮半飽和常數(shù)，而藻類生長速率、藻類基礎(chǔ)新陳代謝速率、藻類捕食速率的敏感度較低。故今后對于其他大型水庫水質(zhì)模型率定時可對敏感度較低的參數(shù)采用經(jīng)驗參數(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] 方開泰.均勻設(shè)計與均勻設(shè)計表[M].北京:科學(xué)出版社,1994.

[3] 燕 喬,吳長彬,張 巖. 基于均勻設(shè)計的鄧肯E-B模型參數(shù)敏感性分析[J]. 中國農(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] 李云生,劉偉江,吳悅穎,等. 美國水質(zhì)模型研究進展綜述[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é)報，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)及面源伴生過程模擬研究與應(yīng)用[D]. 上海：東華大學(xué)，2012.

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

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

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