文章編號:025322468(2003) 0620781205 　　　中圖分類號:X32 　　　文獻(xiàn)標(biāo)識碼:A
杜鵬飛,鄧義祥1 ,陳吉寧　(清華大學(xué)環(huán)境科學(xué)與工程系,北京　100084)

摘要:采用SCE-UA 方法和RSA 方法,比較了在不同的觀測誤差條件下,優(yōu)化方法和不確定性分析方法對于非線性模型參數(shù)識別的影響. 分析表明,RSA 方法是在具有觀測誤差的數(shù)據(jù)條件下進(jìn)行參數(shù)識別的一種有效手段. 同時(shí),通過比較分析,發(fā)現(xiàn)RSA 方法預(yù)測的濃度平均值和概率最大值與“真實(shí)值”并不完全一致. 因此,在使用RSA 方法時(shí),應(yīng)該充分考慮預(yù)測區(qū)間,以降低決策風(fēng)險(xiǎn).
關(guān)鍵詞:觀測誤差;不確定性分析;SCE-UA;RSA

Effects of observation errors on parameter identif ication of a nonlinear environmental model
DU Pengfei ,DENG Yixiang ,CHEN Jining(Department of Environmental Sciences and Engineering Tsinghua University ,Beijing　100084)

Abstract :Parameter identification is a crucial procedure for model application. Observation errors are widely present in environmental data. Thispaper focused on the effects of observation errors on parameter identification of nonlinear environmental models. It specifically compared the identification results by the optimization and uncertainty analysis approaches , taking the SCE-UA and RSA methods as representatives respectively. It was found that the RSA was an effective approach to reduce the impacts of observation errors. It was also showed that the idntified average values or the statistically most likely values of the RSA were often not in consistence with the“true values”. It is there fore suggested that the intervals be considered as identified parameter uncertainty so as to reduce the decision2making risks.
Keywords :observation errors ; uncertainty analysis ; environmental model ; SCE-UA; RSA

作者簡介:杜鵬飛(1970 —) ,男,副教授E2mail :dupf @tsinghua. edu. cn

　　數(shù)學(xué)模型是環(huán)境規(guī)劃與管理的重要手段,參數(shù)是數(shù)學(xué)模型必不可少的組成部分. 一般情況下,數(shù)學(xué)模型的參數(shù)必須通過識別加以確定. 因此,參數(shù)識別是數(shù)學(xué)模型應(yīng)用的重要步驟. 然而,由于環(huán)境系統(tǒng)的復(fù)雜性和環(huán)境監(jiān)測條件的相對滯后,人們可以獲得的數(shù)據(jù)往往十分有限,并且不可避免地存在著觀測誤差,給數(shù)學(xué)模型的參數(shù)識別帶來了嚴(yán)重的影響.
　　優(yōu)化方法是傳統(tǒng)的參數(shù)識別技術(shù),在環(huán)境領(lǐng)域以及其它領(lǐng)域都有比較廣泛的應(yīng)用. 優(yōu)化技術(shù)的研究成果不斷涌現(xiàn),由Duan 等人提出的SCE-UA 方法是比較公認(rèn)的具有較高效率的優(yōu)化方法之一. RSA 方法是由Hornberger 和Spear 等人提出,它是基于隨機(jī)采樣的不確定性分析方法. 本文將采用上述2 種方法研究觀測誤差對非線性環(huán)境模型參數(shù)識別的影響.
1 　模型描述
　　考慮一個(gè)簡單的子過程方程,即污染物衰減模型:
　　　　　d C/d t = - k′C (1)
　　式中, C 為污染物的濃度(mg/L) , k′為衰減系數(shù)(1/d) , t 為時(shí)間(d) .
　　這里假設(shè)k′不是常數(shù),而與污染物濃度有關(guān):
　　　　　k′= kC/( C + M) (2)
　　式中, k 是系數(shù)(1/d) ,M 是限制因子(mg/L) .
　　根據(jù)式(2) ,當(dāng)C m M 時(shí), k′近似為常數(shù);當(dāng)C n M 時(shí),則k′近似與C 成正比.
　　聯(lián)立(1) (2) 兩式: d C/d t = - kC2/( C + M) (3)
　　積分可得: - ln C + M/C = kt - ln C0 + M/C0 (4)
　　其中, C0 為初始時(shí)刻污染物的濃度(mg/L) .
　　注意到式(4) 中C 是t 的單調(diào)函數(shù),因此在k 、M 和C0 已知的情況下, C 是可唯一確定的.因此記C 為t 的函數(shù):
　　　　　C = C( C0 , k ,M , t) (5)
　　選擇該模型作為研究對象,主要考慮到以下幾個(gè)方面的因素:首先,它是許多環(huán)境模型的組成部分,在實(shí)際中有廣泛的應(yīng)用;其次,模型是動態(tài)和非線性的,參數(shù)之間具有較強(qiáng)的相互作用;最后,模型的結(jié)構(gòu)比較簡單,參數(shù)之間的關(guān)系容易把握.
2 　數(shù)據(jù)序列的產(chǎn)生
　　本文采用合成的“觀測”數(shù)據(jù)進(jìn)行研究,即在已知參數(shù)值的情況下采用式(4) ,產(chǎn)生“真實(shí)值”時(shí)間序列,然后加上隨機(jī)撓動以模擬觀測誤差,得到“觀測值”時(shí)間序列,然后利用這些數(shù)據(jù)進(jìn)行模型的識別和驗(yàn)證. 這樣做的目的,是為了使系統(tǒng)的真實(shí)情況在掌握之中,同時(shí)模型沒有結(jié)構(gòu)上的誤差,參數(shù)估計(jì)的所有的誤差僅來源于初值和觀測系列. 這使得有可能在排除結(jié)構(gòu)誤差的情況下,單獨(dú)研究觀測誤差對于參數(shù)識別的影響. 合成數(shù)據(jù)方法在許多模型分析中有廣泛的使用^{[1 ]} . 為簡潔起見,在不引起混淆的情況下,不再對“觀測值”和“真實(shí)值”加上引號.
　　令C0 = 20 mg/L , k = 0．2／d ,M = 10 mg／L ,由于方程(4) 是一個(gè)非線性方程,因此對于給定的時(shí)刻t ,采用牛頓方法解非線性方程求出C 值^{[2 ]} ,再加上均值為0、標(biāo)準(zhǔn)方差為1 m/g 的隨機(jī)撓動,即撓動項(xiàng)為: ξ～ N (0 ,1) (6)
　　觀測值C′為: C′(t) = C( t) + ξ (7)
　　顯然,觀測誤差滿足獨(dú)立、正態(tài)和同方差的性質(zhì).
3 　目標(biāo)函數(shù)
　　由于具有多個(gè)時(shí)刻的采樣值,因此必須將多目標(biāo)轉(zhuǎn)化為單目標(biāo). 考慮到上述觀測誤差是正態(tài)的獨(dú)立同分布,此時(shí)最大似然方法與最小二乘法是等價(jià)的,即目標(biāo)函數(shù)為模擬值與實(shí)測值差的平方和. 大量的經(jīng)驗(yàn)表明,如果觀測誤差滿足上述關(guān)系,則最小二乘法目標(biāo)函數(shù)是最優(yōu)的目標(biāo)函數(shù). 由于模型沒有結(jié)構(gòu)誤差,本文不再對方法進(jìn)行相關(guān)性分析.
　　由于C0 的特殊性,這里考慮將C0 作為參數(shù)和不作為參數(shù)2 種方案. 不識別C0 時(shí),參數(shù)取值空間為Ω = { ( k ,M) | 0 ≤k ≤2 ,0 ≤M ≤100} ;識別C0 時(shí),參數(shù)取值空間為Ω = { ( k , M ,C0 ) | 0 ≤k ≤2 ,0 ≤M ≤100 ,15 ≤C0 ≤30}.
4 　優(yōu)化方法的識別與分析
4．1 　SCE－UA 方法
　　優(yōu)化方法是傳統(tǒng)的參數(shù)識別技術(shù)之一,它的形式很多,收斂速度各有不同,但本質(zhì)是一樣的^{[3 ,4 ]} . 本文采用SCE－UA 優(yōu)化方法.
　　SCE方法是Duan 等人于1992 年在吸收了控制隨機(jī)搜索方法和遺傳算法的基礎(chǔ)上提出 的^{[5 ]} ,1993 年對其進(jìn)行了修正,稱為SCE-UA 方法^{[6 ,7 ]} . SCE-UA 方法在水文學(xué)模型中得到了廣泛應(yīng)用,是到目前為止對于非線性復(fù)雜模型采用隨機(jī)搜索方法尋找最優(yōu)值最為成功的方法之一^{[5 —8 ]} . SCE 方法的主要步驟為:
　　步驟0 　初始化. 選擇p ≥1 和m ≥n + 1 ,這里p 是復(fù)雜形的個(gè)數(shù), m 是每一個(gè)復(fù)雜形的點(diǎn)數(shù),計(jì)算樣本數(shù)s = p ×m.
　　步驟1 　產(chǎn)生樣本. 在可行空間Ω< Rn 采取s 個(gè)點(diǎn)x1 、x2 、.、xs . 計(jì)算在點(diǎn)xi 處的函數(shù)值f i . 在缺少前驗(yàn)信息的條件下,采用均勻采樣分布.
　　步驟2 　對點(diǎn)進(jìn)行排序. 將s 個(gè)點(diǎn)以升序進(jìn)行排列,將它們存貯到數(shù)組中: D = { xi , f i , i =1 , ., s} ,因此i = 1 代表了目標(biāo)最小的函數(shù)點(diǎn).
　　步驟3 　復(fù)雜形劃分. 將D 劃分為p 個(gè)復(fù)雜形A1 , .,Ap ,每一個(gè)包括m 個(gè)點(diǎn),以使:
　　
　　步驟4 　根據(jù)競爭復(fù)雜形演化算法對每一個(gè)復(fù)雜形進(jìn)行演化.
　　步驟5 　混合復(fù)雜形. 將A1 , .,Ap 替代到D 中,以使D = { Ak , k = 1 , ., p} ,對D 按目標(biāo)函數(shù)的升序進(jìn)行排列.
　　步驟6 　檢查收斂性. 如果滿足收斂的標(biāo)準(zhǔn),則結(jié)束搜索,否則回到步驟3.
　　有關(guān)SCE－UE 方法更為詳細(xì)的介紹見參考文獻(xiàn).
4．2 　初值不作為參數(shù)
　　采用真實(shí)序列進(jìn)行參數(shù)識別,識別的結(jié)果為k = 0120 ,M = 9198 ,基本上等于模型的真實(shí)參數(shù),說明在沒有觀測誤差的條件下,SCE-UA 方法能夠有效的識別系統(tǒng)的真值. 但如果采用觀測值序列,則識別的結(jié)果分別為k = 1 , M = 8116 ,盡管模擬的結(jié)果尚可,最大相對誤差不超過10 %(見圖1 (a) ) ,但參數(shù)識別結(jié)果與真實(shí)值的誤差高達(dá)400 %和716 %.

　　采用C0 = 1 mg/L 和100 mg/L 的系統(tǒng)狀態(tài)進(jìn)行模型驗(yàn)證,則最大相對誤差分別為1016 %和5418 % ,見圖1 (b) 和圖1 (c) . 可見,當(dāng)模型進(jìn)行外推時(shí),誤差將顯著增大.
4．3 　初值作為參數(shù)
　　當(dāng)k 和M 確定以后,只要C0 確定,模擬時(shí)間序列就唯一確定了. 因此, C0 既是觀測值,也具有參數(shù)的特點(diǎn). 因此可將初值也作為未知參數(shù)進(jìn)行識別. 采用SCE-UA 方法進(jìn)行搜索,3 個(gè)參數(shù)的識別結(jié)果分別為: k = 0126 ,M = 1511 , C0 = 2011 ,誤差分別為3015 %、5017 %和0125 %. 最優(yōu)參數(shù)的擬合曲線與觀測的對比見圖2 (a) .

　　將k 和M 的識別結(jié)果代入到C0 = 1 mg/L 和100 mg/L 進(jìn)行模型驗(yàn)證,見圖2 (b) 和圖2 (c) .與不將初值作為參數(shù)進(jìn)行識別相比,模型識別精度的改進(jìn)非常明顯. 對C0 = 1 mg/L 和100 mg/L的系統(tǒng)狀態(tài)進(jìn)行檢驗(yàn),則最大相對誤差分別為3110 %和19192 %.
5 　RSA 方法的識別與分析
5．1 　RSA 方法
　　20 世紀(jì)70 年代末、80 年代初,由于認(rèn)識到模型識別的困難,Hornberger 和Spear 將過于強(qiáng)硬的優(yōu)化條件弱化,即將其轉(zhuǎn)化為一些可以用定量或定性的語言描述的條件來決定參數(shù)的取舍,從而在一定程度上克服了采用優(yōu)化方法進(jìn)行參數(shù)識別而帶來的不確定性問題,這就是RSA方法. RSA 方法是基于行為和非行為的二元劃分進(jìn)行參數(shù)識別的,換句話說,給定一組參數(shù),如果系統(tǒng)的模擬行為滿足事先設(shè)定的條件,那么這組參數(shù)就是可接受的,否則是不可接受的. 該方法的主要步驟如下:1．確定參數(shù)的初始分布F(θ) . 如果對參數(shù)的分布特性沒有比較確切的了解,可假設(shè)參數(shù)在取值區(qū)間上均勻分布. 2．根據(jù)模擬值和觀測值的對比,確定行為參數(shù)( B)和非行為參數(shù)( NB) 的劃分準(zhǔn)則. RSA 方法對行為參數(shù)的劃分準(zhǔn)則非常靈活. 這提高了RSA 方法的適用范圍. 3．從參數(shù)的初始分布隨機(jī)采樣,取得一組參數(shù),帶入模型進(jìn)行模擬. 如果模擬的結(jié)果是可以接受的,那么認(rèn)為這組參數(shù)是行為參數(shù),是可以接受的;否則是非行為參數(shù),是不可以接受的. 4．重復(fù)步驟2 和3 ,直到取得足夠的樣本數(shù)為止.
　　在參數(shù)空間隨機(jī)采樣1 萬次,取對應(yīng)目標(biāo)函數(shù)最小的5 %參數(shù)為可接受參數(shù). 為提高采樣的效率,采用Latin Hypercube 方法[9 ,10 ] 替代傳統(tǒng)的Monte Carlo 采樣.
5．2 　識別結(jié)果
　　從4．2 節(jié)和4．3 節(jié)的討論可以看出,將初值作為參數(shù)能夠更好地進(jìn)行系統(tǒng)識別,因此RSA方法只考慮將初值作為參數(shù)的情況.
RSA 方法的模型識別結(jié)果見圖3 (a) . 由于采用了可接受行為的定義,所獲得的不再是單一的參數(shù),而是由大量可行參數(shù)所組成的集合,其分布具有明顯的區(qū)域性特點(diǎn).
　　當(dāng)C0 = 1 mg/L 和100 mg/L 時(shí)的模型驗(yàn)證見圖3 (b) 和圖3 (c) . 可以看出,各時(shí)刻物質(zhì)濃度值的預(yù)測區(qū)間覆蓋了真實(shí)值,從而使模型能充分地估計(jì)預(yù)測結(jié)果的風(fēng)險(xiǎn). 但是,與優(yōu)化方法相比,RSA 方法在一定程度上犧牲了模擬精度以換取可靠性,這體現(xiàn)了在參數(shù)識別過程中估計(jì)精度和可靠性之間的矛盾,對其仍然需要進(jìn)一步的研究.
　　從圖3 (b) 和圖3 (c) 可看出,RSA 方法預(yù)測的濃度平均值與真實(shí)值存在一定的偏差. 因此,在應(yīng)用RSA 方法時(shí),使用均值或概率最大值仍具有潛在的風(fēng)險(xiǎn),應(yīng)充分考慮模型預(yù)測結(jié)果的范圍,以降低決策風(fēng)險(xiǎn).
　　　　　　　　　　　　
6 　結(jié)論
　　本文采用SCE-UA 方法和RSA 方法,研究了觀測誤差對模型參數(shù)識別和預(yù)測的影響. 為避免模型結(jié)構(gòu)誤差的影響,引入了合成數(shù)據(jù)方法. 在具有觀測誤差的情況下,采用優(yōu)化方法不能很好地識別模型的參數(shù). 在本例中優(yōu)化方法的分析表明,將初值作為參數(shù)進(jìn)行識別,可提高系統(tǒng)的可識別性. 采用RSA 方法,參數(shù)識別結(jié)果不再是單一值,而是由大量可行參數(shù)所形成的集合. RSA 在一定程度上克服了觀測誤差帶來的預(yù)測風(fēng)險(xiǎn),但模型預(yù)測的精度也因此而降低. 這體現(xiàn)了參數(shù)識別時(shí)預(yù)測精度和可靠性之間的固有矛盾.
參考文獻(xiàn):
[ 1 ] 　Ratto M,Tarantola S ,Saltelli A. Sensitivity analysis in model calibration :GSA2GLUE approach[J ] . Computer Physics Communication ,2001 ,136 :212 —224
[ 2 ] 　徐士良. Fortran 算法常用算法程序集[M] . 北京:清華大學(xué)出版社,1995. 123 —125
[ 3 ] 　弄文訓(xùn),謝金星. 現(xiàn)代優(yōu)化計(jì)算方法[M] . 北京:清華大學(xué)出版社,2001. 1 —245
[ 4 ] 　劉　毅,陳吉寧,杜鵬飛. 環(huán)境模型參數(shù)優(yōu)化方法的比較研究[J ] . 環(huán)境科學(xué),2002 ,23(2) :1 —6
[ 5 ] 　Qingyun Duan ,Soroosh Sorooshian ,Vijai Gupta. Effective and efficient global optimization for conceptual rainfall2runoff models[J ] .Wat Res Res ,1992 ,28 (4) :1015 —1031
[ 6 ] 　Soroosh Sorooshian ,Qingyun Duan ,Vijai Kumar Gupta. Calibration of rainfall2runoff models :application of global optimization to the　sacramento soil moisture accounting model [J ] .Wat Res Res ,1993 ,29 (4) :1185 —1197
[ 7 ] 　Duan D Y,Gupta V K,Sorooshian S. Shuffled comples evolution approach for effective and efficient global minimization[J ] . Journal of　Optimization Theory and Applications ,1993 ,76 (3) :501 —520
[ 8 ] 　Qingyun Duan ,Soroosh Sorooshian ,Vijai K Gupta. Optimal use of the SCE-UA global optimization method for calibration watershed　model [J ] . Journal of Hydrology ,1994 ,158 :265 —284
[ 9 ] 　Jining Chen. Modeling and control of the activated sludge process :towards a systematic framework[D] . The University of London ,　1993. 130 —165
[10 ] 　鄧義祥,陳吉寧,杜鵬飛. HSY算法在水質(zhì)模型參數(shù)識別中的應(yīng)用[J ] . 上海環(huán)境科學(xué),2002 ,21(8) :497 —500