吳?；郏n月琪，王成林，黃 娟

(1.金陵科技學(xué)院基礎(chǔ)部，江蘇 南京 211169;2.解放軍理工大學(xué)氣象學(xué)院，江蘇 南京 211101;3.江蘇省環(huán)境監(jiān)測(cè)中心，江蘇 南京 210036)

·解析評(píng)價(jià)·

污染模型中參數(shù)訂正的集合Kalman濾波擴(kuò)展?fàn)顟B(tài)變量法

吳?；?，韓月琪2，王成林2，黃 娟3

(1.金陵科技學(xué)院基礎(chǔ)部，江蘇 南京 211169;2.解放軍理工大學(xué)氣象學(xué)院，江蘇 南京 211101;3.江蘇省環(huán)境監(jiān)測(cè)中心，江蘇 南京 210036)

污染模型中不確定參數(shù)的精確訂正對(duì)于提高模型的精度有著重要的意義。在集合Kalman濾波(EnKF)同化方法的基礎(chǔ)上，提出了對(duì)模型中不確定參數(shù)進(jìn)行訂正的EnKF擴(kuò)展?fàn)顟B(tài)變量法，將不確定參數(shù)看成和模型狀態(tài)變量一樣的量，根據(jù)觀測(cè)資料對(duì)不確定變量進(jìn)行訂正，以達(dá)到訂正參數(shù)的目的。采用一個(gè)簡(jiǎn)化的空氣質(zhì)量方程，對(duì)模型參數(shù)訂正方案進(jìn)行檢驗(yàn)，結(jié)果證明提出的方案可行和有效。同時(shí)發(fā)現(xiàn)，隨著觀測(cè)資料精度的提高，無(wú)論是參數(shù)還是模型的狀態(tài)變量，估計(jì)分析值的精度也得到相應(yīng)的提高。

污染模型;參數(shù)訂正;集合Kalman濾波;擴(kuò)展?fàn)顟B(tài)變量法

20世紀(jì)70年代以來(lái)，國(guó)內(nèi)外污染的治理實(shí)踐表明，污染的負(fù)荷定量化研究是控制、評(píng)價(jià)和管理點(diǎn)源、非點(diǎn)源污染的基礎(chǔ)，而模型模擬已成為定量研究點(diǎn)源、非點(diǎn)源污染負(fù)荷最直接和有效的方法［1-3］。國(guó)內(nèi)外學(xué)者針對(duì)污染估算開(kāi)發(fā)了大量的數(shù)學(xué)模型，如 ANSWERS，CREAMS，WEPP，AGNPS，SWRRB/SWAT，COHERENS 模型等［4-11］，國(guó)內(nèi)外利用模型進(jìn)行污染以及水文水資源的模擬研究成為熱點(diǎn)，并取得了一定的成果［12-19］。雖然這些模型在試驗(yàn)中都取得了較好的模擬結(jié)果，但地域特征明顯，對(duì)不同的區(qū)域條件需要重新確定模型參數(shù)。這些參數(shù)的確定往往決定了模型在各地域使用的成敗［20，21］。

在這些數(shù)值模型中，都有一些經(jīng)驗(yàn)系數(shù)或常用參數(shù)，針對(duì)不同的地域特點(diǎn)或下邊界條件，這些參數(shù)會(huì)發(fā)生改變。如果還是沿用原來(lái)的參數(shù)值，就會(huì)導(dǎo)致模型的模擬結(jié)果和實(shí)際情況產(chǎn)生較大的偏差。因此，如何根據(jù)實(shí)際觀測(cè)資料，對(duì)模型中的不確定參數(shù)進(jìn)行合理的估計(jì)訂正，對(duì)提高污染模型的模擬結(jié)果精度具有重要的意義。

對(duì)于污染模型中的不確定參數(shù)，目前主要是根據(jù)文獻(xiàn)資料、典型區(qū)調(diào)查和現(xiàn)場(chǎng)實(shí)驗(yàn)的結(jié)果來(lái)決定。這些方法基本上都是定性的，只能給出一個(gè)大概的結(jié)果。為了根據(jù)觀測(cè)資料對(duì)污染模型中的不確定參數(shù)進(jìn)行定量的訂正，使參數(shù)值更接近實(shí)際情況，引入集合 Kalman濾波(EnKF)數(shù)據(jù)同化方法［22］。EnKF同化方法具有隨流型演變的背景誤差協(xié)方差陣，不需要求解模型的切線性及其伴隨方程的優(yōu)點(diǎn)，受到許多人的檢驗(yàn)、改進(jìn)和應(yīng)用［23］，成為近10年來(lái)數(shù)據(jù)同化方向研究的熱點(diǎn)。但EnKF同化方法主要是根據(jù)觀測(cè)資料對(duì)模型狀態(tài)變量進(jìn)行估計(jì)，還不能直接用來(lái)訂正估計(jì)模式中的不確定參數(shù)，因?yàn)檫@些不確定參數(shù)并不是模型狀態(tài)變量。筆者對(duì)這種方法進(jìn)行了改進(jìn)，提出可以根據(jù)觀測(cè)資料對(duì)模式中不確定參數(shù)進(jìn)行訂正的擴(kuò)充狀態(tài)變量法，即把這些不確定參數(shù)當(dāng)成和模型狀態(tài)變量一樣，然后利用觀測(cè)資料進(jìn)行同化訂正。

1 EnKF 的原理方法

EnKF是Evensen根據(jù)Epstein的隨機(jī)動(dòng)力預(yù)報(bào)理論［22，24］，提出的用 Monte Carlo 方法和 Kalman估計(jì)方法相結(jié)合的集合資料同化方法［25］。EnKF的基本思想是:根據(jù)背景場(chǎng)和觀測(cè)場(chǎng)誤差統(tǒng)計(jì)性質(zhì)，隨機(jī)產(chǎn)生有限樣本的集合擾動(dòng)，然后將這些擾動(dòng)分別加到背景場(chǎng)和觀測(cè)場(chǎng)上，從而產(chǎn)生背景場(chǎng)和觀測(cè)場(chǎng)集合，并用背景場(chǎng)的預(yù)報(bào)集合樣本來(lái)統(tǒng)計(jì)計(jì)算預(yù)報(bào)誤差協(xié)方差矩陣，最后利用Kalman最小方差估計(jì)理論對(duì)新時(shí)刻的觀測(cè)資料進(jìn)行同化，得到一組分析場(chǎng)集合，用此分析場(chǎng)集合作預(yù)報(bào)，到下一個(gè)觀測(cè)時(shí)刻再同化，如此循環(huán)。公式如下:

用一組模式預(yù)報(bào)集合擾動(dòng)估計(jì)的背景場(chǎng)誤差協(xié)方差矩陣為:

2 用于不確定參數(shù)訂正的擴(kuò)展?fàn)顟B(tài)變量EnKF的方法

前面介紹EnKF分析過(guò)程，主要是針對(duì)模型狀態(tài)變量估計(jì)的。如果要根據(jù)觀測(cè)資料對(duì)模型中的不確定參數(shù)和模型狀態(tài)變量一并進(jìn)行估計(jì)，就需要采用擴(kuò)展?fàn)顟B(tài)變量的方法。即如何把模型中的不確定參數(shù)當(dāng)成和模型狀態(tài)變量一樣，然后再利用EnKF方法對(duì)這些擴(kuò)展后的“模型狀態(tài)變量”進(jìn)行估計(jì)分析。具體過(guò)程如下:

針對(duì)包含不確定參數(shù)的數(shù)值模型，公式為:

式中:x——原模型狀態(tài)變量;a——模型中的參數(shù);w——模型誤差。可以把(7)式改寫(xiě)為如下公式:

再利用EnKF方法，就可以根據(jù)觀測(cè)資料對(duì)包括不確定參數(shù)在內(nèi)的“狀態(tài)變量”進(jìn)行訂正。

之所以能夠用擴(kuò)展?fàn)顟B(tài)向量法對(duì)數(shù)值模型中的不確定參數(shù)進(jìn)行訂正，是因?yàn)椴淮_定參數(shù)的值會(huì)對(duì)每一時(shí)刻狀態(tài)變量的值產(chǎn)生影響，而它們之間的影響關(guān)系，可以通過(guò)Monte Carlo集合預(yù)報(bào)的方法統(tǒng)計(jì)出來(lái)，具體表現(xiàn)形式就是擴(kuò)展后狀態(tài)變量的誤差協(xié)方差矩陣B∧b。有了這個(gè)關(guān)系，采用EnKF同化方法就可以根據(jù)觀測(cè)資料值對(duì)模型中不確定的參數(shù)進(jìn)行訂正估計(jì)。下面通過(guò)一個(gè)空氣質(zhì)量變化的數(shù)值試驗(yàn)對(duì)提出的方案進(jìn)行檢驗(yàn)。

3 一個(gè)空氣質(zhì)量模型的數(shù)值試驗(yàn)

下面是描述空氣質(zhì)量變化的一個(gè)簡(jiǎn)單的、無(wú)量綱概念模型［26］:

式中:Q——隨時(shí)間變化的污染源;C——污染物的濃度;a——降解系數(shù)。試驗(yàn)中假設(shè)真實(shí)的污染源為Qtrue(t)=sin(t)，降解系數(shù)真值為0.2。將真實(shí)初始濃度Ctrue(0)=1.0在時(shí)間段t∈［0，5］的積分值視為真值，積分的時(shí)間步長(zhǎng)取0.1，試驗(yàn)中不考慮模型誤差。

觀測(cè)值是在真實(shí)濃度上加白噪音，觀測(cè)誤差均方差為0.1。集合成員個(gè)數(shù)為100。這里將濃度的初猜值設(shè)為Cf(0)=0.2，誤差均方差設(shè)為0.3;實(shí)際模式中降解系數(shù)為0，誤差均方差為0.1。在下面的試驗(yàn)中都沒(méi)有考慮污染源排放先驗(yàn)誤差的時(shí)間相關(guān)。

為了對(duì)前面提出的不確定參數(shù)訂正的擴(kuò)展變量EnKF方案進(jìn)行檢驗(yàn)，首先進(jìn)行第一組試驗(yàn)，共包含4個(gè)試驗(yàn)。試驗(yàn)1為真值;試驗(yàn)2為不作En-KF資料同化;試驗(yàn)3為只對(duì)模型狀態(tài)變量(即空氣污染物濃度)進(jìn)行訂正;試驗(yàn)4為對(duì)模型的狀態(tài)變量和不確定參數(shù)(即降解系數(shù))同時(shí)進(jìn)行訂正。

圖1為第一組各個(gè)試驗(yàn)中污染物濃度隨時(shí)間的變化。由圖中可以發(fā)現(xiàn)，如果不考慮模型中不確定參數(shù)的誤差而只考慮狀態(tài)變量的誤差，即使進(jìn)行了觀測(cè)資料同化，數(shù)值模擬結(jié)果仍然有較大的誤差。如果根據(jù)觀測(cè)資料對(duì)狀態(tài)變量和不確定參數(shù)同時(shí)進(jìn)行訂正，則數(shù)值模擬結(jié)果具有較大的改善。這說(shuō)明在參數(shù)存在誤差而不訂正的情況下，在同化或數(shù)值模擬過(guò)程中會(huì)對(duì)污染物狀態(tài)變量的計(jì)算造成很大的誤差。因此根據(jù)觀測(cè)資料對(duì)模式中的不確定參數(shù)進(jìn)行訂正對(duì)于準(zhǔn)確模擬污染物的擴(kuò)散和改進(jìn)模式的性能具有重要的意義。另外，數(shù)值結(jié)果表明對(duì)不確定參數(shù)訂正的擴(kuò)展變量EnKF方案不但可以根據(jù)觀測(cè)資料改善模式的狀態(tài)變量模擬效果，而且可以根據(jù)觀測(cè)資料對(duì)模式中的不確定參數(shù)同時(shí)進(jìn)行訂正。

圖1 第一組試驗(yàn)中污染物濃度隨時(shí)間的變化

圖2為試驗(yàn)4模式中不確定參數(shù)降解系數(shù)a的估計(jì)值隨時(shí)間的變化。從圖中可以看出，隨著估計(jì)次數(shù)的增多，降解系數(shù)逐漸向真值0.2靠近，最后穩(wěn)定在0.18左右。

圖2 降解系數(shù)估計(jì)值隨時(shí)間的變化

由以上的試驗(yàn)結(jié)果可以看出，對(duì)參數(shù)進(jìn)行估計(jì)的擴(kuò)展?fàn)顟B(tài)變量法是有效的。通過(guò)觀測(cè)資料，在對(duì)狀態(tài)變量進(jìn)行估計(jì)的同時(shí)，此種方法可以對(duì)數(shù)值模型中的不確定參數(shù)進(jìn)行訂正。當(dāng)然從前面的試驗(yàn)還可以看出，無(wú)論是狀態(tài)變量(即污染物濃度)還是降解系數(shù)，它們的估計(jì)值最后都還和真值有一定的偏差。為此我們把觀測(cè)資料的方差定為0.001(相比第一組試驗(yàn)提高了觀測(cè)資料的精度)，進(jìn)行了第二組試驗(yàn)。試驗(yàn)還是4個(gè)試驗(yàn)，設(shè)計(jì)同第一組。

圖3為第二組試驗(yàn)中各個(gè)試驗(yàn)中污染物濃度隨時(shí)間的變化，圖4為不同試驗(yàn)情況下，狀態(tài)變量(即污染物濃度)的誤差函數(shù)隨時(shí)間的變化。其中誤差函數(shù)定義為:

式中:X——污染物濃度，是第二組中4個(gè)不同試驗(yàn)情況下的計(jì)算值;t——真解，由真實(shí)的初始場(chǎng)和降解系數(shù)時(shí)間積分得到。

從圖3、4中可以得出和第一組試驗(yàn)相同的結(jié)論，另外與第一組試驗(yàn)中的圖1相比，還可以看出在對(duì)模型的狀態(tài)變量和不確定參數(shù)同時(shí)進(jìn)行訂正的情況下，因?yàn)橛^測(cè)資料精度的提高，狀態(tài)變量污染物濃度的估計(jì)值與真值的偏差減小，估計(jì)值的精度得到了提高。這說(shuō)明提高觀測(cè)資料的精度，可以有效改善估計(jì)變量(狀態(tài)變量和參數(shù))的同化效果。

圖3 第二組試驗(yàn)中污染物濃度隨時(shí)間的變化

圖4 誤差函數(shù)隨時(shí)間的變化

圖5為觀測(cè)資料方差為0.001的情況下作為不確定參數(shù)的降解系數(shù)估計(jì)值隨時(shí)間的變化。與圖2相比可以發(fā)現(xiàn)，隨著觀測(cè)資料精度的提高降解系數(shù)估計(jì)值的精度也得到相應(yīng)的提高。

圖5 降解系數(shù)估計(jì)值隨時(shí)間的變化

4 結(jié)論與討論

污染數(shù)值模型中不確定參數(shù)的訂正對(duì)于提高數(shù)值模型的模擬精度有著重要的意義。該文在EnKF同化方法的基礎(chǔ)上，提出了EnKF擴(kuò)展?fàn)顟B(tài)變量法來(lái)對(duì)模型中不確定的參數(shù)進(jìn)行訂正，即將不確定參數(shù)看成和狀態(tài)變量一樣的量，然后根據(jù)觀測(cè)資料對(duì)這些變量進(jìn)行訂正。采用一個(gè)簡(jiǎn)化的空氣質(zhì)量方程對(duì)所提出的參數(shù)訂正方案進(jìn)行檢驗(yàn)，結(jié)果證明提出的方案是可行且有效的。并且發(fā)現(xiàn)，隨著觀測(cè)資料精度的提高，無(wú)論是不確定參數(shù)還是狀態(tài)變量，訂正分析值的精度也得到相應(yīng)的提高。

當(dāng)然這種參數(shù)訂正估計(jì)的方法可以用在統(tǒng)計(jì)回歸模型系數(shù)的確定上。文中只是采用一個(gè)簡(jiǎn)單的空氣質(zhì)量模型對(duì)所提出的參數(shù)訂正方案進(jìn)行檢驗(yàn)。對(duì)于更復(fù)雜的污染數(shù)值模型，因?yàn)榭紤]的參數(shù)、試用的條件均不同，該文提出的訂正方案是否有效以及如何將這種方案應(yīng)用到更復(fù)雜的模型中還有待進(jìn)一步研究。

［1］ 郝芳華，楊勝天，程紅光，等.大尺度區(qū)域非點(diǎn)源污染負(fù)荷計(jì)算方法［J］.環(huán)境科學(xué)學(xué)報(bào)，2006，26(3):375-383.

［2］ 鄭東海，王凌，曾紅娟，等.松濤水庫(kù)流域非點(diǎn)源污染負(fù)荷模擬模型［J］.環(huán)境科學(xué)學(xué)報(bào)，2009，29(6):1311-1320.

［3］ 李強(qiáng)坤，李懷恩，胡亞偉，等.黃河干流潼關(guān)斷面非點(diǎn)源污染負(fù)荷估算［J］.水科學(xué)進(jìn)展，2008，19(4):460-466.

［4］ BEASLEY D B，HUGGIN L F.ANSWERS user′s manual［M］.West Layette:Purdue University，1982:33-40.

［5］ KNISEL W G，CREAM S.A field scale model for chemicals，runoff and erosion from agriculture management system［R］.Washington，D.C.:USDA，1983.

［6］ FLANAGAR D C，NEARING G R，LANE L J.USDA-WEPP:hillslope profile and watershed model documentation［R］.Washington，D.C.:USDA-ARS，1995.

［7］ YOUNG R A.AGNPS:A nonpoint source pollution model for evaluating agricultural watershed ［J］.Journal of Soil and Water Conservation，1989，44(2):168-173.

［8］ ECKHARDT K，HAVERKAMP S，F(xiàn)OHRER N，et al.SWAT-G，a version of SWAT99.2 modified for app lication to low mountain range catchments［J］.Physics and Chemistry of the Earth:Parts A/B/C，2002，27:641-644.

［9］ HUA Z L，F(xiàn)AN X P，JiING Z Y，et al.Three-Dimensional Numerical Simulation of Flow Field and Pollutant Transport in Haizhou Bay［C］.Proceeding of 15thCongress of the Asia and Pacific Division of the International of Hydraulic Engineering and Research International Symposium on Maritime Hydraulics，Allied Publishers Pvt.Ltd，2006:1437-1448.

［10］ 李艷蕓，李紹武.風(fēng)暴潮預(yù)報(bào)模式在渤海海域中的應(yīng)用研究［J］.海洋技術(shù)，2006，25(1):101-106.

［11］ 華祖林，顧莉，查玉含，等.基于COHERENS模型的污染物質(zhì)輸運(yùn)數(shù)值模擬［J］.環(huán)境科學(xué)與技術(shù)，2009，32(4):14-18.

［12］ IOANNIS P，MARIA M，MARIA K.Estimation of nitrogen and phosphorus losses to surface water and groundwater through the implementation of the SWAT model for Norwegian soils［J］.Journal of Soils and Sediments，2007，7(4):223-231.

［13］ DEBELE B，SRINIVASAN R，PARLANGE J Y.Coupling up land watershed and downstream waterbody hydrodynamic and water quality models(SWAT and CE-QUAL-W2)for better water resources management in complex river basins［J］.Environmental Modeling and Assessment，2008，13(1):135-153.

［14］ 王中根，劉昌明，黃友波.SWAT模型的原理結(jié)構(gòu)及應(yīng)用研究［J］.地理科學(xué)進(jìn)展，2003，22(1):79-87.

［15］ 楊玨，錢(qián)新，張玉超，等.兩種新型流域非點(diǎn)源污染負(fù)荷估算模型的比較［J］.中國(guó)環(huán)境科學(xué)，2009，29(7):762-766.

［16］ 李強(qiáng)坤，李懷恩，胡亞偉，等.黃河干流潼關(guān)斷面非點(diǎn)源污染負(fù)荷估算［J］.水科學(xué)進(jìn)展，2008，19(4):460-466.

［17］ 秦耀民，胥彥玲，李懷恩.基于SWAT模型的黑河流域不同土地利用情景的非點(diǎn)源污染研究［J］.環(huán)境科學(xué)學(xué)報(bào)，2009，29(2):440-448.

［18］ 薛罡，劉亞男，汪永輝，等.曝氣充氧條件下受污染河道的水質(zhì)模型建立及應(yīng)用［J］.環(huán)境科學(xué)，2010，31(3):653-659.

［19］ 鄭東海，王凌，曾紅娟，等.松濤水庫(kù)流域非點(diǎn)源污染負(fù)荷模擬模型［J］.環(huán)境科學(xué)學(xué)報(bào)，2009，29(6):1311-1320.

［20］ JOHNES P J.Evaluation and management of the impact of land use change on the nitrogen and phosphorus load delivered to surface waters:the export coefficient modeling approach ［J］.Journal of Hydrology，1996，183:323-349.

［21］ 蔡明，李懷恩，莊詠濤，等.改進(jìn)的輸出系數(shù)法在流域非點(diǎn)源污染負(fù)荷估算中的應(yīng)用［J］.水利學(xué)報(bào)，2004(7):40-45.

［22］ EVENSEN G.Sequential data assimilation with a nonlinear Quasigeostrophic model using Monte Carlo methods to forecast error statistics［J］.Journal of Geophysical Research Oceans，1994，99:10143-10162.

［23］ SNYDER C，ZHANG F.Assimilation of simulated Doppler radar observations with an ensemble Kalman filter［J］.American Meteorological Society，2003，131:1663-1677.

［24］ EPSTEIN E S.Stochastic dynamic prediction［J］.Tellus，1969，21:739-759.

［25］ 林彩燕，朱江，陸春谷.集合KALMAN濾波和最優(yōu)插值方法在不同觀測(cè)分布的比較理想試驗(yàn)［J］.氣候與環(huán)境研究，2006，5:553-563.

［26］ 朱江，汪萍.集合卡爾曼平滑和集合卡爾曼濾波在污染源反演中的應(yīng)用［J］.大氣科學(xué)，2006，30(5):871-882.

The Research of Extended State Variable Method in Ensemble Kalman Filter for Estimating Uncertain Parameters in Pollution Model

WU Zhu-hui1，HAN Yue-qi2，WANG Cheng-lin2，HUANG Juan3
(1.Basis Department，Jinling Institute of Technology，Nanjing，Jiangsu 211169，China;2.Institute of Meteorology，PLA University of Science and Technology，Nanjing，Jiangsu 211101，China;3.Jiangsu Provincial Environmental Monitoring Center，Nanjing，Jiangsu 210036，China)

The exact estimation of uncertain parameters in pollution model makes great sense in enhance the precision of numerical model.The method of extended state variable based on Ensemble Kalman Filter(EnKF)is introduced to estimate the uncertain parameters.That is the uncertain parameters as the model state variables that can be corrected by observational data.A simple quality equation of air is used to test the method of correcting model parameters.The results of experiment show that the method is feasible and effective.At the same time，the precision of estimated value of parameters and state variables is improved with the elevation of the observational data′s precision.

pollution model;parameter estimation;Ensemble Kalman Filter;extended state variable method

X823

B

1674-6732(2012)-03-0036-05

10.3969/j.issn.1674-6732.2012.03.009

2010-12-01;

2010-12-18

國(guó)家自然科學(xué)基金項(xiàng)目(40805046);江蘇省自然科學(xué)基金項(xiàng)目(BK2010128);公益性行業(yè)(氣象)專項(xiàng)課題(GYHY(QX)2007-6-15，GYHY200906009)。

吳?；?1977—)，女，講師，碩士，從事數(shù)理統(tǒng)計(jì)研究。