劉毅,陳吉寧,杜鵬飛(清華大學(xué)環(huán)境科學(xué)與工程系環(huán)境模擬與污染控制國(guó)家重點(diǎn)聯(lián)合實(shí)驗(yàn)室,北京　
100084 ,E-mail :jchen1 @mail. tsinghua. edu. cn)

摘要:在對(duì)水文模型實(shí)例的參數(shù)不確定性分析基礎(chǔ)上,分別采用傳統(tǒng)靈敏度分析方法、HSY算法、線性回歸等方法對(duì)模型參數(shù)特性進(jìn)行了識(shí)別與比較研究. 結(jié)果表明參數(shù)優(yōu)化算法與傳統(tǒng)靈敏度分析方法不能解釋模型結(jié)構(gòu)復(fù)雜性特征,采用不確定性分析方法對(duì)環(huán)境模型參數(shù)進(jìn)行識(shí)別提供了深入分析與理解模型系統(tǒng)的有效途徑.
關(guān)鍵詞:參數(shù)優(yōu)化;不確定性;模型結(jié)構(gòu);參數(shù)識(shí)別;靈敏度;HSY算法
中圖分類號(hào):X11 　文獻(xiàn)標(biāo)識(shí)碼:A 　文章編號(hào):025023301 (2002) 0620520006
基金項(xiàng)目:高等學(xué)校優(yōu)秀青年教師教學(xué)科研獎(jiǎng)勵(lì)計(jì)劃資助項(xiàng)目
作者簡(jiǎn)介:劉毅(1975～) , 男, 博士研究生, 主要從事環(huán)境系統(tǒng)分析方向的研究.
收稿日期:2001210212 ;修訂日期:2002205208

　　隨著環(huán)境模型結(jié)構(gòu)復(fù)雜性的急劇增長(zhǎng),模型參數(shù)在高維空間表現(xiàn)出了復(fù)雜的相關(guān)性結(jié)構(gòu)并直接導(dǎo)致了優(yōu)化后驗(yàn)參數(shù)的識(shí)別問題[1 , 2 ] .因此僅僅局限于參數(shù)優(yōu)化算法效率和精度等方面的研究已經(jīng)不能滿足理論與實(shí)踐的需要. 相對(duì)于觀測(cè)數(shù)據(jù)和模型參數(shù)而言,基于現(xiàn)有科學(xué)認(rèn)知體系構(gòu)建的模型結(jié)構(gòu)是建模過程中不確定性的根本來源[3 ] . 然而由于缺乏深入研究結(jié)構(gòu)不確定性的理論基礎(chǔ)和有效技術(shù)手段,模型結(jié)構(gòu)往往只能通過識(shí)別參數(shù)后驗(yàn)統(tǒng)計(jì)分布規(guī)律間接地得到驗(yàn)證.
　　參數(shù)不確定性依賴于模型結(jié)構(gòu),并直接導(dǎo)致了靈敏度問題[4 ] . 模型參數(shù)不確定性包括參數(shù)可識(shí)別性和參數(shù)靈敏度2 個(gè)基本方面. 參數(shù)不確定性分析提供了后驗(yàn)地識(shí)別模型結(jié)構(gòu)的可行途徑[5 , 6 ] . 本文以一個(gè)經(jīng)典的水文箱式模型為實(shí)例,從優(yōu)化參數(shù)的不確定性和參數(shù)識(shí)別問題出發(fā),對(duì)復(fù)雜模型的參數(shù)靈敏度分析方法進(jìn)行了比較研究.
1 　環(huán)境模型參數(shù)不確定性與識(shí)別
　　由于參數(shù)不確定性普遍存在,根據(jù)經(jīng)驗(yàn)估計(jì)或者觀測(cè)值優(yōu)化得到的參數(shù)并不能保證模型應(yīng)用的精度和預(yù)測(cè)結(jié)果的可靠性. 在追求更加高效和穩(wěn)定優(yōu)化算法的同時(shí),所得到優(yōu)化參數(shù)的后驗(yàn)分布具有本質(zhì)上的不確定性;對(duì)于復(fù)雜模型來說,甚至無法判斷優(yōu)化結(jié)果是否達(dá)到了全局最優(yōu),也無法預(yù)測(cè)“最優(yōu)”參數(shù)對(duì)于模型預(yù)測(cè)的影響[7 ] . 引入?yún)?shù)不確定性分析可以更為深刻地理解和認(rèn)識(shí)到現(xiàn)實(shí)世界和模型系統(tǒng)特征之間的本質(zhì)區(qū)別.
　　直接參數(shù)優(yōu)化算法可以理解為在一定控制性準(zhǔn)則約束下的空間隨機(jī)搜索過程. 算法根據(jù)新點(diǎn)產(chǎn)生器和判定準(zhǔn)則在高維參數(shù)空間持續(xù)地進(jìn)行“產(chǎn)生新點(diǎn)2判斷2接受/ 舍棄”的迭代過程,直至滿足終止準(zhǔn)則. 為了說明參數(shù)不確定性及其識(shí)別問題,本文采用4 種直接參數(shù)優(yōu)化算法進(jìn)行比較分析[8 ] . 控制隨機(jī)搜索算法(CRS) 引入了幾何學(xué)中“重心”的概念,即考慮了新點(diǎn)產(chǎn)生的隨機(jī)性,又在一定程度上保證了搜索的整體性. 復(fù)合形混合演化算法( SCE- UA) 是將生物自然演化過程引入到數(shù)值計(jì)算中,模擬了生物進(jìn)化的過程,提高了計(jì)算效率和全局搜索整體最優(yōu)的能力. 模擬退火算法(SA) 則假設(shè)優(yōu)化問題的解及其目標(biāo)函數(shù)分別與固體物質(zhì)的微觀狀態(tài)及其能量所對(duì)應(yīng),采用隨機(jī)方法模擬固體穩(wěn)定“退火”的過程. 退火單純形算法(AS) 綜合了下山單純形方法和模擬退火法2 種優(yōu)化算法,更加充分地利用了單純形的形變信息,從而提高了計(jì)算效率和算法穩(wěn)定性.
　　隨著環(huán)境模型的不斷開發(fā)和廣泛應(yīng)用,環(huán)境模型的種類和數(shù)量日益豐富,模型本身所表現(xiàn)出的結(jié)構(gòu)特征也日趨復(fù)雜. 本研究?jī)H以一個(gè)兩箱式水文模型為例[8 ] ,表1 給出了模型中11個(gè)參數(shù)的先驗(yàn)取值范圍及其物理意義.
　　　　　　　　　　　　　

　　優(yōu)化算法的基本思想認(rèn)為在特定模型結(jié)構(gòu)下只有唯一一組最佳參數(shù)與之對(duì)應(yīng). 然而由于模型結(jié)構(gòu)復(fù)雜性與輸入數(shù)據(jù)的不確定性,優(yōu)化算法通常不能尋優(yōu)到環(huán)境模型的唯一真值,并且也無法判斷算法是否達(dá)到全局最優(yōu). 這是導(dǎo)致優(yōu)化參數(shù)不確定性的根本原因. 根據(jù)4 種直接優(yōu)化算法得到的優(yōu)化參數(shù)結(jié)果可知,在目標(biāo)函數(shù)值無顯著性差異的條件下( max δ <118 %) , “最佳”估計(jì)參數(shù)之間具有較大差異(見表2) . 由于本文采用的水文模型輸入數(shù)據(jù)序列非常完整,可以忽略其不確定性,并且?guī)追N優(yōu)化算法也是當(dāng)前最為穩(wěn)定和可靠的全局搜索算法[8 ] ,因此上述優(yōu)化結(jié)果直接證實(shí)了優(yōu)化算法不能為深入研究復(fù)雜環(huán)境模型提供有效途徑.
　　進(jìn)一步研究表明,不同算法尋優(yōu)進(jìn)程中的參數(shù)收斂軌跡也具有較大差異[8 ] . 以CRS 算法運(yùn)用于水文模型為例,圖1 是500 個(gè)最佳優(yōu)化參數(shù)樣本點(diǎn)分布圖,其中菱形符號(hào)代表最佳的10 個(gè)參數(shù)值分布(采樣總數(shù)n = 20000) ,虛線對(duì)應(yīng)于最佳目標(biāo)函數(shù)值.
　　　　　　　　　　　　　　　　
　　可以由圖直觀地看到最佳10 組參數(shù)對(duì)應(yīng)的目標(biāo)函數(shù)值沒有顯著區(qū)別,但這10 組參數(shù)值并沒有非常一致地收斂到“最佳”參數(shù)估計(jì)值.這種優(yōu)化結(jié)果通常稱為優(yōu)化參數(shù)的等效性或可置換性[7 , 9 ] . 優(yōu)化參數(shù)可置換性的產(chǎn)生是由于模型參數(shù)在高維空間具有的復(fù)雜相關(guān)性,是模型結(jié)構(gòu)復(fù)雜性和參數(shù)不確定性的集中體現(xiàn). 這種與優(yōu)化算法結(jié)構(gòu)設(shè)計(jì)、新點(diǎn)產(chǎn)生器、空間搜索方法和接受/ 舍棄判定準(zhǔn)則等高度相關(guān)的后驗(yàn)參數(shù)空間分布上的差異性,是導(dǎo)致優(yōu)化算法不能解釋模型結(jié)構(gòu)及其所產(chǎn)生的參數(shù)不確定性的主要原因.
2 　環(huán)境模型的參數(shù)靈敏度
　　如前所述,優(yōu)化方法不能解釋模型結(jié)構(gòu)復(fù)雜性與參數(shù)不確定性,由此產(chǎn)生了對(duì)于模型參數(shù)的識(shí)別問題. 參數(shù)識(shí)別很重要的一個(gè)方面就是研究參數(shù)變化所引起的模型響應(yīng),即參數(shù)靈敏度問題[5 , 9 ] . 研究參數(shù)靈敏性有助于深入理解并改進(jìn)模型結(jié)構(gòu)的穩(wěn)定性.
　　傳統(tǒng)參數(shù)靈敏度分析方法是在某個(gè)參數(shù)最佳估計(jì)值附近給定一個(gè)人工干擾,并計(jì)算參數(shù)在這一很小范圍內(nèi)產(chǎn)生波動(dòng)所導(dǎo)致模型輸出的變化率,即擾動(dòng)分析方法. 從表2 的計(jì)算結(jié)果中可以得出基于4 種直接優(yōu)化算法最佳參數(shù)估計(jì)值的靈敏度數(shù)值(δ= ±10 %) 具有較大差異,其結(jié)果強(qiáng)烈依賴于優(yōu)化算法的選擇.
　　　　　
　　如前所述,由于優(yōu)化算法自身結(jié)構(gòu)設(shè)計(jì)決定了最佳參數(shù)估計(jì)值并得到具有特定收斂特征的參數(shù)樣本,因此基于“最佳”估計(jì)參數(shù)值的靈敏度分析不能完整地描述模型參數(shù)的空間分布形態(tài);另一方面更為本質(zhì)的是,由模型結(jié)構(gòu)復(fù)雜性導(dǎo)致的參數(shù)相關(guān)性要求在靈敏度分析過程中必須考慮參數(shù)之間的相互影響,而不是僅靠變動(dòng)某一個(gè)參數(shù)得到模型響應(yīng). 考慮到模型參數(shù)之間的高度相關(guān)性,現(xiàn)代環(huán)境系統(tǒng)研究在不確定性分析思想框架下提出了更為有效的參數(shù)靈敏度分析方法, 即區(qū)域靈敏度分析方法( Regional Sensitivity) [1 , 3 ] . 與傳統(tǒng)分析方法不同,區(qū)域靈敏度分析拋棄了“尋優(yōu)”思想,承認(rèn)參數(shù)空間分布的復(fù)雜性與相關(guān)性,是一個(gè)對(duì)模型結(jié)構(gòu)在一定準(zhǔn)則下對(duì)隨機(jī)參數(shù)大樣本發(fā)生響應(yīng)的統(tǒng)計(jì)分析過程. 本文采用線性回歸和HSY 2 種區(qū)域靈敏度方法.
　　線性回歸方法是通過隨機(jī)采樣產(chǎn)生參數(shù)樣本序列,計(jì)算每個(gè)樣本對(duì)應(yīng)的模型響應(yīng),然后進(jìn)行線性回歸(LR) , 建立如下形式的線性模 型[6 ] :
　　　　　　　　
　　(1) 式中,βi 為對(duì)應(yīng)于參數(shù)pi 的線性回歸模型系數(shù),表征了參數(shù)pi 對(duì)模型響應(yīng)y 的貢獻(xiàn)率(權(quán)重) ,即為基于線性回歸方法得到的參數(shù)pi的絕對(duì)靈敏度. 由于模型參數(shù)量綱不同,通常采用(2) 式計(jì)算參數(shù)相對(duì)靈敏度βi( s) , 其中σpi和σy 分別表示參數(shù)樣本和模型輸出的方差:
　　　　　　　　　　
　　HSY算法用于區(qū)域靈敏度分析是基于模型某種形式的判定準(zhǔn)則,通過隨機(jī)采樣過程產(chǎn)生2組參數(shù),并根據(jù)統(tǒng)計(jì)檢驗(yàn)方法計(jì)算參數(shù)靈敏度[3 ] . 這里的某一準(zhǔn)則是根據(jù)已有經(jīng)驗(yàn)、數(shù)據(jù)或是在一定置信區(qū)間下的模擬準(zhǔn)則,用來判斷隨機(jī)采樣的參數(shù)樣本是否被接受,體現(xiàn)了不確定性分析中不承認(rèn)“最優(yōu)”的基本思想. 可接受與不可接受2 組參數(shù)樣本共同表征了模型本質(zhì)特征[1 , 5 ] .　　本文中采用(3) 式定義可接受的系統(tǒng)行為:Yobv[1 - αβ( n) ] ≤ Ysim ≤ Yobv[1 +αβ( n) ](3)
　　式中, n 是模型優(yōu)化結(jié)果與觀測(cè)值之間的誤差向量,δ( n) 是n 維向量的樣本方差,α表示可接受系統(tǒng)偏差的倍數(shù). Yobv和Ysim分別表示觀測(cè)值和模型輸出值. 這樣就形成了一個(gè)系統(tǒng)輸出的控制“廊道”,凡是超出這個(gè)廊道邊界的樣本,都認(rèn)為是不可接受的采樣(見圖2) .
　　　　　　　　　　　　　　
　　兩樣本非參數(shù)假設(shè)檢驗(yàn)的方法有χ22擬合優(yōu)度檢驗(yàn), K2S 檢驗(yàn),秩檢驗(yàn),游程檢驗(yàn)等,不同方法對(duì)靈敏度計(jì)算結(jié)果會(huì)產(chǎn)生一定影響[10 ] . 本文采用K2S 檢驗(yàn)方法[1 , 5 , 10 ] . K2S 檢驗(yàn)數(shù)值越大,說明可接受參數(shù)的概率分布與不可接受參數(shù)的概率分布之間的差別越大,從而參數(shù)靈敏性越大,具有明確物理意義的參數(shù)也越容易被識(shí)別. 當(dāng)涉及到大樣本采樣過程時(shí),可以采用下面的準(zhǔn)則來判斷參數(shù)靈敏度分析算法是否達(dá)到收斂[1 ] : ①重復(fù)采樣過程,樣本均值和方差不發(fā)生變化; ②重復(fù)采用過程,模型參數(shù)的靈敏度排序不發(fā)生變化.
　　　　　
　　表３給出了采用LR 與HSY算法得到水文模型實(shí)例中的11 個(gè)參數(shù)的區(qū)域靈敏度數(shù)值及排序,作為比較同時(shí)列出了參數(shù)的局部靈敏度排序. 由于HSY算法考慮了模型結(jié)構(gòu)導(dǎo)致的參數(shù)相關(guān)性,并且采用完全隨機(jī)采樣方法,因此其區(qū)域靈敏度分析結(jié)果具有可靠性. 而LR 算法是在整個(gè)參數(shù)可能空間上采樣,參數(shù)靈敏度中不加區(qū)分地包含了可接受與不可接受參數(shù)的特征信息;另外,盡管算法也采用了在整個(gè)參數(shù)空間上的隨機(jī)過程,但是從建立參數(shù)與模型響應(yīng)之間的線性模型中得到的參數(shù)靈敏度又不可避免地抵消了模型結(jié)構(gòu)自身復(fù)雜性的影響,因此其靈敏度結(jié)果是不可靠的.
　　盡管局部靈敏度數(shù)值之間存在較大差異,但某些參數(shù)的局部靈敏度排序卻出現(xiàn)一定程度上的一致性,例如局部最靈敏參數(shù)均為h1 . 如前所述,由于僅考慮了單一參數(shù)在“最優(yōu)點(diǎn)”附近沿特定方向變化所引起的模型輸出響應(yīng),參數(shù)局部靈敏度不能反映出模型結(jié)構(gòu)特征與參數(shù)相關(guān)性. 而HSY算法提出將參數(shù)劃分為系統(tǒng)可接受與不可接受樣本的思想,從根本上避免了傳統(tǒng)“尋優(yōu)”思想帶來的模型結(jié)構(gòu)和參數(shù)不確定性分析中的困惑,并通過計(jì)算2 參數(shù)樣本之間的統(tǒng)計(jì)距離,得到了包含參數(shù)空間分布特征整體的區(qū)域靈敏度,為優(yōu)化模型結(jié)構(gòu)和改進(jìn)模型穩(wěn)定性提供了有力的分析基礎(chǔ)和技術(shù)支持.
　　以部分參數(shù)為例,圖3 進(jìn)一步給出了HSY算法可接受參數(shù)的最佳500 個(gè)估計(jì)值與4 種優(yōu)化算法的最佳估計(jì)參數(shù)的空間分布比較. 取代唯一“最優(yōu)”的參數(shù)估計(jì)值,HSY 算法以降低計(jì)算效率為代價(jià)保證了參數(shù)空間搜索的整體性,避免了優(yōu)化算法定向搜索對(duì)參數(shù)特性的片面理解;并以降低計(jì)算精度為代價(jià)得到的系統(tǒng)可接受的參數(shù)空間,對(duì)應(yīng)著可接受的模型輸出,其意義在于定量地給出參數(shù)的可識(shí)別性和模型預(yù)測(cè)結(jié)果的可靠性. 特別是在模型應(yīng)用于情景分析中,傳統(tǒng)方法只能根據(jù)參數(shù)局部靈敏度絕對(duì)數(shù)值大小來計(jì)算未來各種可能情景下的模型響應(yīng),由于實(shí)際環(huán)境系統(tǒng)復(fù)雜性和參數(shù)之間的實(shí)際相關(guān)性,優(yōu)化算法并不能保證模型預(yù)測(cè)的可靠性. 而采用區(qū)域靈敏度方法得到的可接受參數(shù)樣本,定量地給出了可接受模型輸出的參數(shù)空間分布,從而保證了模型整體上的預(yù)測(cè)可靠性. 正如有學(xué)者指出,與其面對(duì)沒有可靠性的優(yōu)化算法,不如接受一個(gè)更為合理的不確定性[3 ] .
　　　　　　　　　　　　　　　　
　　　　　　　　　　　(菱形、正方形、三角形與圓形分別表示CRS、SCE - UA、SA 和AS 算
　　　　　　　　　　　法的最佳參數(shù)估計(jì)值,空間點(diǎn)為HSY算法可接受參數(shù)的最佳500 個(gè)估計(jì)值)
　　進(jìn)一步地,采用K2S 統(tǒng)計(jì)檢驗(yàn)方法所得到的區(qū)域靈敏度數(shù)值表征了參數(shù)可識(shí)別性,其為今后模型結(jié)構(gòu)不確定性研究以及模型結(jié)構(gòu)的簡(jiǎn)化和改進(jìn)提供了科學(xué)分析依據(jù)[5 ] .
　　圖4 是基于各種區(qū)域靈敏度分析方法的參數(shù)靈敏度排序分布圖. 基于前文分析,認(rèn)為只有HSY算法得到的參數(shù)靈敏度是無偏的. 無論是參數(shù)優(yōu)化還是參數(shù)靈敏度分析,優(yōu)化算法都不能很好地解釋環(huán)境模型復(fù)雜性和參數(shù)識(shí)別問題. 基于不確定性分析思想的HSY算法提供了識(shí)別復(fù)雜模型參數(shù)和深入理解模型結(jié)構(gòu)的有效方法.
　　　　　　　　　　　　　
參考文獻(xiàn):
1 　Chen J . The Modeling and Control of activated sludeg process : Towards a systematic framework. Ph. D. thesis. Imperial College , London , Engl. , U. K. 1993. 54～94.
2 　Chen J , Wheater H S. Identification and uncertainty analysis　of soil water retention models using lysimeter data. Water　Resources Research , 1999 , 35 (8) : 2401～2414.
3 　Beck M B. Water quality modeling : A review of the analysis　of uncertainty. Water Resources Research , 1987 , 23 ( 8) :1393～1442.
4 　Spear R C , Grieb T M , Shang N. Parameter uncertainty and　interaction in complex environmental models. Water Resources Research , 1994 , 30 (11) : 3159～3169.
5 　Chen J , Beck M B. Quanlity Assurance of Multi2Media　Model for Predictive Screening Tasks. USEPA , Report .1998.
6 　Klepper O. Multivariate aspects of model uncertainty analysis : Subjective calibration of incorrect models. Report . National Institute of Public Health and Environmental Protection , The Netherlands. 1999.
7 　Aalderink R H , Jovin J . Identification of the parameters describing primary production from continuous oxygen signals.　Water Science and Technology , 1997 , 36 (5) : 43～51.
8 　劉毅, 陳吉寧, 杜鵬飛. 環(huán)境模型參數(shù)優(yōu)化方法的比較研究. 環(huán)境科學(xué), 2002 , 23 (2) : 1～6.
9 　Beven K. Prophecy , reality and uncertainty in distributed　hydrological modeling. Advances in Water Resources , 1993 ,16 : 41～51.
10 　Reckhow K H , Clements J T , Dodd R C. Statistical evaluation of mechanistic water2quality models. Journal of Environmental. Engineering , 1990 , 116 (2) : 250～268.