基于模型參數(shù)不確定性的河流突發(fā)污染事故動態(tài)預(yù)測方法
【專利摘要】本發(fā)明公開了基于模型參數(shù)不確定性的河流突發(fā)污染事故動態(tài)預(yù)測方法。S1.利用參數(shù)率定、經(jīng)驗數(shù)據(jù)以及查閱文獻(xiàn)的方法生成不確定性參數(shù)；S2.利用S1中不確定性參數(shù)值和初始斷面污染物濃度計算事故點以后各斷面各時刻的污染物濃度值；S3.選取似然函數(shù)計算S2中不同參數(shù)組所對應(yīng)的似然值；S4.利用S3中的似然值可以估算出在一定置信度水平下模型預(yù)測結(jié)果的不確定范圍，擴(kuò)展到事故發(fā)生以后的整個事件維度可以獲得模型預(yù)測結(jié)果的不確定性區(qū)間；S5.利用實測數(shù)據(jù)對模擬預(yù)測結(jié)果不斷進(jìn)行更新和校正，重復(fù)上述步驟，可求得最新更新的預(yù)測結(jié)果。該方法結(jié)合不確定性方法、動態(tài)更新理論、普適似然不確定性算法，實現(xiàn)了基于模型參數(shù)不確定性的河流突發(fā)污染事故水質(zhì)預(yù)測。
【專利說明】
基于模型參數(shù)不確定性的河流突發(fā)污染事故動態(tài)預(yù)測方法
技術(shù)領(lǐng)域
[0001] 本發(fā)明涉及污染物模擬領(lǐng)域，具體提出一種基于普適似然不確定性算法和動態(tài)更 新的突發(fā)性水污染事故預(yù)警方法。
【背景技術(shù)】
[0002] 河流、湖泊以及水庫是人類重要的淡水生態(tài)資源，沿江臨湖地區(qū)也通常是人類活 動頻繁和各種生物生息繁衍的重要區(qū)域。然而，隨著經(jīng)濟(jì)的高速發(fā)展，人類的生產(chǎn)和生活活 動對水資源安全造成了嚴(yán)重影響，這引發(fā)了眾多經(jīng)濟(jì)問題和社會問題。目前，尚無有效手段 徹底阻止人類活動對于水環(huán)境的惡劣影響，在今后相當(dāng)長的時間范圍內(nèi)，水環(huán)境問題仍將 是我們不得不面對的突出問題。
[0003] 在眾多水環(huán)境問題中，河流突發(fā)污染事故是一類發(fā)生頻率高，危害性嚴(yán)重的事故。 河流突發(fā)污染事故通常具有以下特征：突發(fā)、不易預(yù)知，不確定性高以及危害巨大等，這使 得突發(fā)污染事故發(fā)生后的水質(zhì)預(yù)測工作和應(yīng)急處理工作充滿困難。為使應(yīng)急決策者更好地 應(yīng)對突發(fā)污染事故，對事故發(fā)生后污染物濃度時空變化規(guī)律的描述顯得尤為重要。
[0004] 污染事故的發(fā)生、發(fā)展以及演變具有很大的不確定性，主要表現(xiàn)為：（1)發(fā)生時間、 空間的不確定性；（2)污染物的不確定性；（3)事故流域特性的不確定性；（4)污染方式的不 確定性；（5)事故記錄數(shù)據(jù)的不確定性等。這些不確定性的客觀存在為污染物濃度的預(yù)測工 作帶來了困難和挑戰(zhàn)。
[0005] 確定性水質(zhì)模型雖然理論系統(tǒng)完善、甚至可以準(zhǔn)確細(xì)致地表現(xiàn)污染物迀移擴(kuò)散過 程，但由于河流環(huán)境本身的復(fù)雜性使得機(jī)理建模很難準(zhǔn)確，而初始條件和水文數(shù)據(jù)又很難 全面獲知，再加上模型參數(shù)的率定工作很難進(jìn)行，尤其是突發(fā)污染事故要求模型盡快給出 預(yù)測數(shù)據(jù)，使得短期建立起優(yōu)秀的確定性水質(zhì)模型變得非常困難。在這樣的背景下，隨著不 確定性理論研究的逐步成熟，不確定性水質(zhì)模型逐漸成為了相關(guān)工作者的研究熱點。
[0006] 不確定性水質(zhì)模型主要包括以統(tǒng)計分析為主要研究方法的隨機(jī)理論模型，以時間 序列分析、人工神經(jīng)網(wǎng)絡(luò)、遺傳算法等為代表的主要基于數(shù)據(jù)處理的數(shù)學(xué)模型，以及以灰色 系統(tǒng)、模糊數(shù)學(xué)為處理手段的灰色預(yù)測模型和模糊模型。
[0007] 目前關(guān)于參數(shù)或者變量方面的不確定性分析，用的最多的是蒙特卡羅方法，該方 法基于隨機(jī)抽樣進(jìn)行模擬，其結(jié)果與參數(shù)分布情況的關(guān)系較大，該方法適用于各類復(fù)雜的 非線性系統(tǒng)。

【發(fā)明內(nèi)容】

[0008] 為了解決現(xiàn)有技術(shù)的不足，本發(fā)明提供了一種基于模型參數(shù)不確定性的河流突發(fā) 污染事故動態(tài)預(yù)測方法。
[0009] -種基于模型參數(shù)不確定性的河流突發(fā)污染事故動態(tài)預(yù)測方法，包括以下步驟：
[0010] S1.不確定性參數(shù)組的生成；
[0011] S2.模擬預(yù)測：利用S1中不確定性參數(shù)值和初始斷面污染物濃度計算事故點以后 各斷面各時刻的污染物濃度值；
[0012] S3.似然值計算:選取似然函數(shù)計算S2中污染物濃度值所對應(yīng)的似然值；
[0013] S4.概率密度函數(shù)計算:利用S3中的似然值估算出在模型預(yù)測結(jié)果的不確定范圍， 擴(kuò)展到事故發(fā)生以后的整個事件維度獲得模型預(yù)測結(jié)果的不確定性區(qū)間；
[0014] S5.利用實測濃度數(shù)據(jù)比較模擬預(yù)測結(jié)果后作出判斷，如果實測濃度數(shù)據(jù)與模擬 預(yù)測結(jié)果之間具有誤差，則進(jìn)入S2，重復(fù)上述步驟，得到更新的預(yù)測結(jié)果；
[0015] 如果實測濃度數(shù)據(jù)與模擬預(yù)測結(jié)果相符，則結(jié)束。
[0016] 所述的S1步驟中，根據(jù)模型的不確定性參數(shù)，并確定其初始分布;對于選定的河流 突發(fā)污染事故模型，模型的不確定性參數(shù)為E，u，k，其中E為縱向彌散系數(shù)，u為縱向水流流 速，k為綜合衰減速率系數(shù)，E，u，k的初始分布取均勻分布，根據(jù)歷史數(shù)據(jù)和相關(guān)經(jīng)驗分別確 定E，u，k的參數(shù)取值范圍，然后按照蒙特卡羅方法隨機(jī)獲取E，u，k的組合0，0表示如下：
[0017] 9 = [RandE,Randu,Randk] (1)
[0018]式中，Rand代表對模型參數(shù)按照其分布進(jìn)行隨機(jī)取值。
[0019]所述的方法，將初始斷面各時刻的污染物濃度值和參數(shù)組0代入公式(2)，迭代算 得事故點以后各斷面各時刻的污染物濃度值，代入不同參數(shù)組Mi = l，2,…，N)，N為蒙特 卡羅模擬次數(shù)，分別獲得與其對應(yīng)的各斷面各時刻的污染物濃度值序列，再利用蒙特卡羅 方法進(jìn)行統(tǒng)計得到河流中任意位置、任意時刻的污染物濃度值分布，即概率密度函數(shù)；
[0021 ]式中：C為平均污染物濃度，mg/L; u為縱向水流流速，單位為km2/h ;E為縱向彌散系 數(shù)，單位為km2/h; x為河水流動距離，單位為km; k為綜合衰減速率系數(shù)，單位為IT1; t為時間， 單位為h;i表示第i段河段，j來表示第j個時間間隔；AX和At分別表示空間和時間間隔大 小。
[0022]所述的方法，過程如下：
[0023]設(shè)有統(tǒng)計獨立的隨機(jī)變量Xi(i = 1，2，…，k)，其對應(yīng)的概率密度函數(shù)分別為fxl， 1^2，'"，；1^，功能函數(shù)式為￥=；1；'(叉1，叉2，"，叉1〇;
[0024]首先根據(jù)各隨機(jī)變量&的分布，產(chǎn)生N組隨機(jī)數(shù)X1，X2,…，處值，然后計算功能函數(shù) 值Yj = f(xi，X2,…，Xk)( j = l，2,…，N)，接下來根據(jù)計算結(jié)果，等間隔地選取m組間隔為S的Yj 的取值區(qū)間，m組區(qū)間互相連接、不重疊，并且包含全部L取值，再分別統(tǒng)計落入每個取值區(qū) 間的Yi個數(shù)，即Quantity (Y.i)，并按式(3)統(tǒng)計變量密度，
(3)
[0026]當(dāng)N足夠大，s足夠小，m足夠大時，根據(jù)伯努利大數(shù)定理和正態(tài)隨機(jī)變量的特性以 及頻率分布直方圖與概率密度函數(shù)之間的關(guān)系有:全部Dens i ty (Yj)構(gòu)成了 Yj的概率密度函 數(shù)，而該概率密度函數(shù)即為被估計變量Y的實際取值分布。
[0027]似然函數(shù)表征某個參數(shù)組所預(yù)測出來的濃度值與實測濃度值之間的差異程度，采 用Nash-Sutcliffe確定性系數(shù)： (4)
[0029] 式中：L代表似然值;R2代表確定性系數(shù);of代表預(yù)測序列的誤差方差;W代表實測 序列的方差;計算似然值之前需要首先獲取斷面的實測污染物濃度值序列，結(jié)合預(yù)測值，代 入公式(4)即可得到不同參數(shù)組0:所對應(yīng)的似然值。
[0030] 按照似然值的大小對不同參數(shù)組0:的預(yù)測結(jié)果進(jìn)行統(tǒng)計加權(quán)，即按照式3統(tǒng)計時， Quantity(Yj)不再代表落入某個取值區(qū)間的結(jié)果數(shù)目，而是代表計算結(jié)果落入某個取值區(qū) 間的似然值之和，相應(yīng)的N代表全部似然值之和，進(jìn)而獲得各斷面各時刻的預(yù)測結(jié)果分布函 數(shù);在此基礎(chǔ)上，估算出在置信度水平下模型預(yù)測結(jié)果的不確定性范圍，擴(kuò)展到事故發(fā)生以 后的整個時間維度獲得模型預(yù)測結(jié)果的不確定性區(qū)間。
[0031] 所述的S5是利用實測數(shù)據(jù)對模擬預(yù)測結(jié)果不斷進(jìn)行更新、校正的過程，其具體計 算過程是:新的數(shù)據(jù)序列獲取以后，先按照步驟5計算參數(shù)組0:的觀測似然值，然后使用貝 葉斯公式對前一次計算得到的似然值結(jié)果，即先驗似然值，進(jìn)行更新，得到后驗似然值結(jié) 果，再按照步驟6進(jìn)行更新后的模擬預(yù)測，貝葉斯公式表示如下：
(5)
[0033]式中：L(Y | 0i)為后驗似然值;L(0i | Y)為觀測似然值;L(0i)為先驗似然值;C為歸一 化因子;Y為預(yù)測因變量(即預(yù)測結(jié)果he:為第i組參數(shù)組；
[0034] 當(dāng)再有新的實測數(shù)據(jù)獲取時，前一次更新過的后驗似然值結(jié)果即變?yōu)楫?dāng)前計算的 先驗似然值結(jié)果，重復(fù)上述步驟，即可求得更新的預(yù)測結(jié)果。
[0035] 本發(fā)明的有益效果：
[0036]給出了河流突發(fā)污染事故不確定性的河流突發(fā)污染事故水質(zhì)預(yù)測一般框架，結(jié)合 蒙特卡羅方法和GLUE算法實現(xiàn)前述框架，然后以此為基礎(chǔ)，引入動態(tài)更新方法對次水質(zhì)預(yù) 測框架進(jìn)行改進(jìn)，并可利用事故發(fā)生后的污染物濃度實時數(shù)據(jù)對模型預(yù)測結(jié)果進(jìn)行動態(tài)校 正，使得模型預(yù)測結(jié)果隨著污染事故的持續(xù)進(jìn)行仍然保持良好的準(zhǔn)確性和可靠性。
【附圖說明】
[0037] 圖1為基于模型參數(shù)不確定性的河流突發(fā)污染事故動態(tài)預(yù)測方法的流程圖；
[0038] 圖2為模型預(yù)測結(jié)果的動態(tài)更新流程圖；
[0039] 圖3為河流Q以及事故發(fā)生點、監(jiān)測點位置示意圖；
[0040] 圖4為模擬效率X隨著似然值臨界值m選取的變化曲線；
[0041] 圖5為監(jiān)測點A的實測值與不確定性預(yù)測區(qū)間；
[0042] 圖6為監(jiān)測點B的實測值與不確定性預(yù)測區(qū)間；
[0043] 圖7為監(jiān)測點C的實測值與不確定性預(yù)測區(qū)間；
[0044] 圖8為監(jiān)測點B的污染物濃度實測值與不確定性預(yù)測區(qū)間對比。 具體實施方案
[0045]下面結(jié)合附圖對本發(fā)明進(jìn)一步描述。
[0046] 基于模型參數(shù)不確定性的河流突發(fā)污染事故動態(tài)預(yù)測的基本步驟如圖1所示，主 要包含:參數(shù)組生成，模擬預(yù)測，似然函數(shù)選取與似然值計算以及概率密度函數(shù)計算。
[0047] 模擬預(yù)測雖然實現(xiàn)了不確定性預(yù)測，但是尚未利用后續(xù)測得的污染物濃度值數(shù)據(jù) 對預(yù)測結(jié)果進(jìn)行動態(tài)校正，其預(yù)測準(zhǔn)確性隨著時間的推移而無法保證。所以，只有不斷利用 后續(xù)測得的濃度值來實時更新、校正預(yù)測結(jié)果才能使模型預(yù)測結(jié)果一直保持較高的準(zhǔn)確 性。因此，在圖1的基礎(chǔ)上，對預(yù)測結(jié)果動態(tài)更新，其過程如圖2所示：當(dāng)沒有新的數(shù)據(jù)對似然 函數(shù)結(jié)果進(jìn)行動態(tài)更新時，直接利用這個似然函數(shù)值進(jìn)行模擬預(yù)測，即圖1的一次計算；當(dāng) 需要新數(shù)據(jù)對似然函數(shù)結(jié)果進(jìn)行動態(tài)更新時，以這個函數(shù)值作為先驗似然值結(jié)合觀測似然 值利用貝葉斯公式計算后驗似然值從而達(dá)到動態(tài)更新的目的。
[0048] -種基于模型參數(shù)不確定性的河流突發(fā)污染事故動態(tài)預(yù)測方法，包括以下步驟：
[0049] S1.不確定性參數(shù)組的生成；
[0050] S2.利用S1中不確定性參數(shù)值和初始斷面污染物濃度計算事故點以后各斷面各時 刻的污染物濃度值；
[0051] S3.選取似然函數(shù)計算S2中不同參數(shù)組所對應(yīng)的似然值；
[0052] S4.利用S3中的似然值可以估算出在一定置信度水平下模型預(yù)測結(jié)果的不確定范 圍，擴(kuò)展到事故發(fā)生以后的整個事件維度可以獲得模型預(yù)測結(jié)果的不確定性區(qū)間；
[0053] S5.利用實測數(shù)據(jù)對模擬預(yù)測結(jié)果不斷進(jìn)行更新和校正，重復(fù)上述步驟，可求得最 新更新的預(yù)測結(jié)果。
[0054]所述的S1步驟中，參數(shù)組生成部分需要首先明確模型有哪些不確定性參數(shù)，并確 定其初始分布。對于選定的河流突發(fā)污染事故模型，模型的不確定性參數(shù)為E，u，k。通常由 于缺少足夠的數(shù)據(jù)支持，E，u，k的初始分布可取均勻分布。所以實際的計算過程是根據(jù)歷史 數(shù)據(jù)和相關(guān)經(jīng)驗分別確定E，u，k的參數(shù)取值范圍，然后按照蒙特卡羅方法隨機(jī)獲取E，u，k的 組合 9，9可表不如下：
[0055] 9 = [ RandE, Randu, Randk ] (1)
[0056]式中，Rand代表對模型參數(shù)按照其分布進(jìn)行隨機(jī)取值。
[0057]將初始斷面各時刻的污染物濃度值和參數(shù)組0代入公式(2)，可以迭代算得事故點 以后各斷面各時刻的污染物濃度值。代入不同參數(shù)組0i(i = l，2,…，N)(N為蒙特卡羅模擬 次數(shù))可以分別獲得與其對應(yīng)的各斷面各時刻的污染物濃度值序列，再利用蒙特卡羅方法 進(jìn)行統(tǒng)計即可得到河流中任意位置、任意時刻的污染物濃度值分布(概率密度函數(shù)）。通常 突發(fā)污染事故發(fā)生后，對于整條河流而言只有部分?jǐn)嗝嬖O(shè)置了采樣點，即只有部分河流位 置有實測數(shù)據(jù)，所以模擬預(yù)測部分更加關(guān)注這些有采樣點的河道位置的污染物濃度值分 布。
[0059]所使用的蒙特卡羅方法其基本原理如下：
[0060]由概率定義知，某事件發(fā)生的概率可以用大量試驗中該事件發(fā)生的頻率來估算， 當(dāng)試驗足夠多時，可以認(rèn)為該事件發(fā)生的頻率即為其發(fā)生的概率。因此，可以先對影響某一 事件是否發(fā)生的隨機(jī)變量進(jìn)行大量隨機(jī)抽樣，然后把這些抽樣值分別代入功能函數(shù)式，并 統(tǒng)計計算結(jié)果，進(jìn)而求得某一事件的發(fā)生概率。蒙特卡羅方法正是基于此思路進(jìn)行分析的， 其計算過程如下：
[0061 ]設(shè)有統(tǒng)計獨立的隨機(jī)變量Xi(i = 1，2，…，k)，其對應(yīng)的概率密度函數(shù)分別為fxl， 1^2，'"，；1^，功能函數(shù)式為￥=；1；'(叉1，叉2，"，叉1〇。
[0062]首先根據(jù)各隨機(jī)變量t的分布，產(chǎn)生N組隨機(jī)數(shù)X1，X2,…，處值，然后計算功能函數(shù) 值Yj = f(xi，X2,…，Xk)( j = l，2,…，N)，接下來根據(jù)計算結(jié)果，等間隔地選取m組間隔為S的Yj 的取值區(qū)間(m組區(qū)間互相連接、不重疊，并且包含全部L取值），再分別統(tǒng)計落入每個取值 區(qū)間的Yj個數(shù)(Quantity (Yj))，并按式⑶統(tǒng)計變量密度，
C3)
[0064]當(dāng)N足夠大、s足夠小、m足夠大時，根據(jù)伯努利大數(shù)定理和正態(tài)隨機(jī)變量的特性以 及頻率分布直方圖與概率密度函數(shù)之間的關(guān)系可有:全部Dens i ty (Yj)構(gòu)成了 Yj的概率密度 函數(shù)，而該概率密度函數(shù)即為被估計變量Y的實際取值分布。
[0065] 蒙特卡羅方法主要用于解決兩類問題，即純數(shù)學(xué)求解問題和隨機(jī)性問題，其對于 解決具有隨機(jī)性或者不確定性的實際問題有很強(qiáng)的適用能力。本申請所研究的污染物濃度 模擬預(yù)測問題屬于隨機(jī)性問題，通過大量的隨機(jī)試驗從模型參數(shù)的先驗分布中抽取參數(shù)代 入河流突發(fā)污染事故預(yù)測模型進(jìn)行運算，然后統(tǒng)計得到污染物濃度在各時空維度下的概率 密度函數(shù)，這是對結(jié)果進(jìn)行進(jìn)一步不確定性分析的基礎(chǔ)。
[0066] 似然函數(shù)表征某個參數(shù)組所預(yù)測出來的濃度值與實測濃度值之間的差異程度，有 多種可供選擇，而如何選取帶有一定的主觀性，通常可采用Nash-Sutcliffe確定性系數(shù)：
(4)
[0068] 式中：L代表似然值;R2代表確定性系數(shù); <代表預(yù)測序列的誤差方差;<代表實測 序列的方差。
[0069] 計算似然值之前需要首先獲取某斷面的實測污染物濃度值序列，結(jié)合前面計算得 到的預(yù)測值，代入公式(4)即可算得不同參數(shù)組0:所對應(yīng)的似然值，然后按照一定標(biāo)準(zhǔn)(似 然值大于某一閾值，閾值的選取帶有一定主觀性)剔除明顯無法模擬污染物擴(kuò)散情況的參 數(shù)組，以備后續(xù)使用。
[0070] 按照似然值的大小對不同參數(shù)組0:的預(yù)測結(jié)果進(jìn)行統(tǒng)計加權(quán)（即按照式3統(tǒng)計時， Quantity(Yj)不再代表落入某個取值區(qū)間的結(jié)果數(shù)目，而是代表計算結(jié)果落入某個取值區(qū) 間的似然值之和，相應(yīng)的N代表全部似然值之和），進(jìn)而獲得各斷面各時刻的預(yù)測結(jié)果分布 函數(shù)。在此基礎(chǔ)上，可以估算出在一定置信度(如95%)水平下模型預(yù)測結(jié)果的不確定性范 圍，擴(kuò)展到事故發(fā)生以后的整個時間維度可以獲得模型預(yù)測結(jié)果的不確定性區(qū)間。
[0071]動態(tài)更新過程是利用實測數(shù)據(jù)對模擬預(yù)測結(jié)果不斷進(jìn)行更新、校正的過程，其具 體計算過程是:新的數(shù)據(jù)序列獲取以后，先按照步驟5計算參數(shù)組0:的觀測似然值，然后使 用貝葉斯公式對前一次計算得到的似然值結(jié)果(先驗似然值)進(jìn)行更新，得到后驗似然值結(jié) 果，再按照步驟6進(jìn)行更新后的模擬預(yù)測。貝葉斯公式在本應(yīng)用中可表示如下：
(5)
[0073 ]式中：L (Y | 9i)為后驗似然值;L (9 i | Y)為觀測似然值;L (9 i)為先驗似然值;C為歸一 化因子;Y為預(yù)測因變量(即預(yù)測結(jié)果he:為第i組參數(shù)組。
[0074] 當(dāng)再有新的實測數(shù)據(jù)獲取時，前一次更新過的后驗似然值結(jié)果即變?yōu)楫?dāng)前計算的 先驗似然值結(jié)果，重復(fù)上述步驟，即可求得最新更新的預(yù)測結(jié)果。
[0075] 實施例
[0076] 以河流Q上游某處發(fā)生的化學(xué)品泄漏污染事故為實例進(jìn)行模擬計算與結(jié)果分析， 各項數(shù)據(jù)源于文獻(xiàn)。事故發(fā)生點以后的河段全長約為160km，將河道均分為320段(每段500 米，即A x = 500m)進(jìn)行模擬計算（A t = 60s)。按照計算流程，首先需要明確不確定性分析中 的參數(shù)分布和范圍。根據(jù)大量歷史數(shù)據(jù)和現(xiàn)場測定結(jié)果可以確定事發(fā)河流流速u的取值范 圍為:0.1 -1.0(m/s);通過參考相關(guān)文獻(xiàn)，確定河流縱向彌散系數(shù)E為：100- 300(m2/s);綜合 降解系數(shù)k為:0.58\10_6-1.74\10_ 6(8_1)。模型參數(shù)￡，11汰的初始分布取為均勻分布。選取 污染事故發(fā)生后，監(jiān)測數(shù)據(jù)記錄相對完整的3個河流斷面(分別記為監(jiān)測點A，B，C)作為研究 對象，河道與采樣斷面示意圖如附圖3所示。
[0077] 1.有效參數(shù)組合的臨界值選取
[0078]計算完所有參數(shù)組所對應(yīng)的似然值后需要根據(jù)似然值的大小對參數(shù)組進(jìn)行取舍， 取舍標(biāo)準(zhǔn)帶有一定的主觀性。可以采用模擬效率(用符號X表示)對取舍標(biāo)準(zhǔn)進(jìn)行分析，模擬 效率表示被選擇用于進(jìn)行不確定性預(yù)測的參數(shù)組占全部由蒙特卡羅方法生成的參數(shù)組的 百分比，其計算公式如下：
(6)
[0080] 式中，X表示模擬效率;m表示似然值選取臨界值;分子表示似然值大于臨界值m的 參數(shù)組數(shù)目；分母表示GLUE算法中由蒙特卡羅方法所生成的全部參數(shù)組數(shù)目。
[0081] 以監(jiān)測點B的實測數(shù)據(jù)計算不同似然值臨界值m選取下的模擬效率X，可以得到如 附圖4所示的曲線。
[0082] 從圖中可見，模擬效率X隨著臨界值m的增大而逐漸下降。若m值選取過大，則X偏 小，即可接受的參數(shù)組數(shù)目過少，此時可理解為尋找"最優(yōu)解"，由前文分析可知這種情況會 忽視各種不確定性的影響，其預(yù)測結(jié)果反而不準(zhǔn)確。若m值選取過小，則X偏大，此時過多考 慮了各種不確定性因素以至于讓一些無法良好模擬當(dāng)前河流污染物擴(kuò)散狀況的參數(shù)組也 對最終預(yù)測結(jié)果有所貢獻(xiàn)，很明顯這種情況也是不合理的。當(dāng)m取值0.5時，模擬效率近似為 60%，此時既讓模型模擬效率保證在50%以上，又使計算結(jié)果包含了曲線切線附近區(qū)域所 對應(yīng)的參數(shù)組(該區(qū)域即使m值變化較小，X值也會發(fā)生較大變化，若不包含該區(qū)域則會丟掉 很多有價值的參數(shù)組），所以選定似然值臨界值為0.5。
[0083] 2.在一定置信度下的預(yù)測與實際結(jié)果對比
[0084]進(jìn)行5000次蒙特卡羅模擬，分別對監(jiān)測點A，B，C取事故發(fā)生后15h，40h，50h以內(nèi)的 實測數(shù)據(jù)對預(yù)測結(jié)果動態(tài)更新，再分別統(tǒng)計監(jiān)測點A，B，C的預(yù)測結(jié)果，可以獲得監(jiān)測點A，B， C在各時間維度下的模型預(yù)測濃度值分布(概率密度函數(shù)）。對于得到的概率密度函數(shù)，可以 計算出在95%置信度下的最優(yōu)置信區(qū)間，即預(yù)測結(jié)果的不確定性范圍，其與實測值的對比 圖分別如附圖5、6、7所示。
[0085]為評價不確定性預(yù)測區(qū)間的優(yōu)劣，引入覆蓋率(CR)的概念。覆蓋率是指預(yù)測區(qū)間 覆蓋實測數(shù)據(jù)的比率，它是最常用的預(yù)測區(qū)間評價指標(biāo)。CR越大，則表示預(yù)測區(qū)間越能包含 全部實測數(shù)據(jù)。
[0086] 從圖中可以看出，大部分的實測值(CRa=97% ;CRB = 92% ;CRc = 90% )都落在置信 度為95%的不確定性預(yù)測區(qū)間以內(nèi)，但也有少部分的實測值落于不確定性預(yù)測區(qū)間外部。 實測值落于預(yù)測區(qū)間外部的原因有很多:采樣測量誤差，歷史數(shù)據(jù)信息偏差，模型本身準(zhǔn)確 性偏差等。
[0087] 而若把置信度設(shè)定為90%，則有相對較少的實測值落于不確定性預(yù)測區(qū)間內(nèi)，但 由于不確定性預(yù)測區(qū)間變小，其預(yù)測不確定性反而降低，有助于模型結(jié)果的預(yù)測工作。總結(jié) 可有，置信度取值越低，則不確定性預(yù)測范圍越小，預(yù)測區(qū)間的可靠性越低，但有利于進(jìn)行 相對準(zhǔn)確的結(jié)果預(yù)測；反之，置信度取值越高，則不確定性預(yù)測范圍越大，預(yù)測區(qū)間的可靠 性越高，但較難對結(jié)果進(jìn)行準(zhǔn)確預(yù)測。所以，在具體選定預(yù)測置信度時，需要綜合考慮各種 方面的因素。
[0088] 3.實際應(yīng)用過程
[0089] 3中計算得到的不確定性預(yù)測結(jié)果是以某監(jiān)測點自身測得的濃度數(shù)據(jù)對該監(jiān)測點 的不確定性預(yù)測結(jié)果進(jìn)行動態(tài)更新，這樣可以得到相對準(zhǔn)確的預(yù)測結(jié)果（以監(jiān)測點B的不確 定性預(yù)測結(jié)果為例，其覆蓋率高達(dá)92%)。然而，河流突發(fā)污染事故發(fā)生后，倘若等到監(jiān)測點 有可用實測數(shù)據(jù)對不確定性預(yù)測結(jié)果進(jìn)行更新時往往時間過久(尤其對于河流中、下游的 監(jiān)測點而言時間更久，以監(jiān)測點B為例需要等到污染事故發(fā)生后約30h才有可用實測數(shù)據(jù)獲 ?。?，而應(yīng)急決策者卻希望更早地得知污染物在整條河流中的擴(kuò)散規(guī)律，以便從容準(zhǔn)備應(yīng)急 方案，避免污染事故造成更大的經(jīng)濟(jì)損失和社會負(fù)面影響。
[0090] 為解決上述矛盾，在實際應(yīng)用時，可采取如下改進(jìn)做法：
[0091] (1)先采用靠近源頭的監(jiān)測點實測數(shù)據(jù)計算模型參數(shù)組的似然值分布，并粗略認(rèn) 為該似然值分布可以描述河流其他斷面的污染物擴(kuò)散規(guī)律，然后即可據(jù)此分布分別計算河 流各關(guān)注斷面(通常被關(guān)注的河流斷面都會設(shè)置監(jiān)測點）的不確定性預(yù)測區(qū)間。
[0092] (2)當(dāng)有更靠近某關(guān)注斷面的監(jiān)測點獲取到實測數(shù)據(jù)時，則以該監(jiān)測點實測數(shù)據(jù) 對應(yīng)的參數(shù)組似然值分布替換之前的似然值分布，并計算不確定性預(yù)測區(qū)間。
[0093] (3)最后，當(dāng)該關(guān)注斷面自身獲取到實測數(shù)據(jù)后，則可如3-樣計算較為準(zhǔn)確的不 確定性預(yù)測結(jié)果。
[0094] 由于河流環(huán)境本身所具有的不確定性，不同河流斷面所對應(yīng)的模型參數(shù)組似然值 分布本應(yīng)不同，但是為了更早地進(jìn)行不確定性預(yù)測，可以先忽略這些差異的影響而進(jìn)行較 為粗略的不確定性預(yù)測，以便讓應(yīng)急決策者盡早得到數(shù)據(jù)支持，并進(jìn)行前期應(yīng)急決策;當(dāng)某 關(guān)注斷面自身獲取到實測數(shù)據(jù)后，即可按照本專利方法進(jìn)行更為準(zhǔn)確地不確定性預(yù)測，此 時的預(yù)測結(jié)果可以幫助決策者校正之前的相關(guān)決策，進(jìn)而更有針對性地應(yīng)對污染事故。
[0095] 以監(jiān)測點B為例，倘若要進(jìn)行準(zhǔn)確地不確定性預(yù)測，需要等到事故發(fā)生后40h才有 足夠的數(shù)據(jù)支持，所需時間明顯過久。實際應(yīng)用時，可以先用監(jiān)測點A獲取的實測濃度數(shù)據(jù) 計算模型參數(shù)組的似然值分布，然后以該分布計算監(jiān)測點B的不確定性預(yù)測區(qū)間，其不確定 性預(yù)測結(jié)果如附圖8
[0096]按照這種辦法，在事故發(fā)生后15h即可對監(jiān)測點B的污染物濃度變化規(guī)律進(jìn)行不確 定性預(yù)測，雖然預(yù)測結(jié)果準(zhǔn)確性較差(覆蓋率不到43%)，但卻足以為應(yīng)急決策者在對污染 物在河流中的擴(kuò)散規(guī)律絲毫不清楚時提供數(shù)據(jù)參考，便于其準(zhǔn)備充足的方案應(yīng)對污染事 故。事故發(fā)生后30h，監(jiān)測點B自身獲取到實測數(shù)據(jù)后，即可據(jù)此計算得到3所述的更為準(zhǔn)確 的不確定性預(yù)測結(jié)果，決策者可以據(jù)此采取更進(jìn)一步的污染事故應(yīng)對策略。
[0097]以上闡述的是本發(fā)明給出的實例，仿真結(jié)果體現(xiàn)了本發(fā)明所提出的技術(shù)方案對于 A河流突發(fā)污染事故的動態(tài)預(yù)測結(jié)果。需要指出的是，本發(fā)明不只限于上述實施例，對于其 他突發(fā)污染事故，采用本發(fā)明的技術(shù)方案也能給出很好的預(yù)警效果。
【主權(quán)項】
1. 一種基于模型參數(shù)不確定性的河流突發(fā)污染事故動態(tài)預(yù)測方法，其特征在于，包括 以下步驟:S1.不確定性參數(shù)組的生成；52. 模擬預(yù)測：利用S1中不確定性參數(shù)值和初始斷面污染物濃度計算事故點以后各斷 面各時刻的污染物濃度值；53. 似然值計算:選取似然函數(shù)計算S2中污染物濃度值所對應(yīng)的似然值；54. 概率密度函數(shù)計算:利用S3中的似然值估算出在模型預(yù)測結(jié)果的不確定范圍，擴(kuò)展 到事故發(fā)生以后的整個事件維度獲得模型預(yù)測結(jié)果的不確定性區(qū)間；55. 利用實測濃度數(shù)據(jù)比較模擬預(yù)測結(jié)果后作出判斷，如果實測濃度數(shù)據(jù)與模擬預(yù)測 結(jié)果之間具有誤差，則進(jìn)入S2，重復(fù)上述步驟，得到更新的預(yù)測結(jié)果；如果實測濃度數(shù)據(jù)與 模擬預(yù)測結(jié)果相符，則結(jié)束。2. 根據(jù)權(quán)利要求1所述的方法，其特征在于，所述的S1步驟中，根據(jù)模型的不確定性參 數(shù)，并確定其初始分布;對于選定的河流突發(fā)污染事故模型，模型的不確定性參數(shù)為E，u，k， 其中E為縱向彌散系數(shù)，u為縱向水流流速，k為綜合衰減速率系數(shù)，E，u，k的初始分布取均勻 分布，根據(jù)歷史數(shù)據(jù)和相關(guān)經(jīng)驗分別確定E，u，k的參數(shù)取值范圍，然后按照蒙特卡羅方法隨 機(jī)獲取E，u，k的組合θ，θ表示如下： Θ = [RandE,Randu,Randk] (1) 式中，Rand代表對模型參數(shù)按照其分布進(jìn)行隨機(jī)取值。3. 根據(jù)權(quán)利要求2所述的方法，其特征在于，將初始斷面各時刻的污染物濃度值和參數(shù) 組Θ代入公式(2)，迭代算得事故點以后各斷面各時刻的污染物濃度值，代入不同參數(shù)組Θ, (? = 1，2，···，Ν)，Ν為蒙特卡羅模擬次數(shù)，分別獲得與其對應(yīng)的各斷面各時刻的污染物濃度 值序列，再利用蒙特卡羅方法進(jìn)行統(tǒng)計得到河流中任意位置、任意時刻的污染物濃度值分 布，即概率密度函數(shù)；式中：C為平均污染物濃度，mg/L; u為縱向水流流速，單位為km2/h; Ε為縱向彌散系數(shù)，單 位為km2/h; X為河水流動距離，單位為km; k為綜合衰減速率系數(shù)，單位為IT1; t為時間，單位 為h;i表示第i段河段，j來表示第j個時間間隔；Δχ和At分別表示空間和時間間隔大小。4. 根據(jù)權(quán)利要求3所述的方法，其特征在于，過程如下:設(shè)有統(tǒng)計獨立的隨機(jī)變量t (i = 1，2,…，k)，其對應(yīng)的概率密度函數(shù)分別為fxl，fx2,…，fxk，功能函數(shù)式為Y = f(xi，x2，···， xk);首先根據(jù)各隨機(jī)變量Xi的分布，產(chǎn)生N組隨機(jī)數(shù)XI，X2，…，xk值，然后計算功能函數(shù)值Yj = f(xi，X2，···，xk)( j = l，2,…，N)，接下來根據(jù)計算結(jié)果，等間隔地選取m組間隔為S的Yj的取 值區(qū)間，m組區(qū)間互相連接、不重疊，并且包含全部乃取值，再分別統(tǒng)計落入每個取值區(qū)間的 Yj個數(shù)，即Quantity (Yj)，并按式(3)統(tǒng)計變量密度，當(dāng)N足夠大，s足夠小，m足夠大時，根據(jù)伯努利大數(shù)定理和正態(tài)隨機(jī)變量的特性以及頻 率分布直方圖與概率密度函數(shù)之間的關(guān)系有:全部Dens i ty (Yj)構(gòu)成了 Yj的概率密度函數(shù)， 而該概率密度函數(shù)即為被估計變量Y的實際取值分布。5. 根據(jù)權(quán)利要求4所述的方法，其特征在于，似然函數(shù)表征某個參數(shù)組所預(yù)測出來的濃 度值與實測濃度值之間的差異程度，采用Nash-Sutcliffe確定性系數(shù)：式中：L代表似然值;R2代表確定性系數(shù);代表預(yù)測序列的誤差方差;erf代表實測序列 的方差;計算似然值之前需要首先獲取斷面的實測污染物濃度值序列，結(jié)合預(yù)測值，代入公 式(4)即可得到不同參數(shù)組Θ,所對應(yīng)的似然值。6. 根據(jù)權(quán)利要求5所述的方法，其特征在于，按照似然值的大小對不同參數(shù)組Θ,的預(yù)測 結(jié)果進(jìn)行統(tǒng)計加權(quán)，即按照式3統(tǒng)計時，Quantity(Yj)不再代表落入某個取值區(qū)間的結(jié)果數(shù) 目，而是代表計算結(jié)果落入某個取值區(qū)間的似然值之和，相應(yīng)的N代表全部似然值之和，進(jìn) 而獲得各斷面各時刻的預(yù)測結(jié)果分布函數(shù);在此基礎(chǔ)上，估算出在置信度水平下模型預(yù)測 結(jié)果的不確定性范圍，擴(kuò)展到事故發(fā)生以后的整個時間維度獲得模型預(yù)測結(jié)果的不確定性 區(qū)間。7. 根據(jù)權(quán)利要求6所述的方法，其特征在于，所述的S5是利用實測數(shù)據(jù)對模擬預(yù)測結(jié)果 不斷進(jìn)行更新、校正的過程，其具體計算過程是:新的數(shù)據(jù)序列獲取以后，先按照步驟5計算 參數(shù)組Θ,的觀測似然值，然后使用貝葉斯公式對前一次計算得到的似然值結(jié)果，即先驗似 然值，進(jìn)行更新，得到后驗似然值結(jié)果，再按照步驟6進(jìn)行更新后的模擬預(yù)測，貝葉斯公式表 示如下：式中：L(Y | Θi)為后驗似然值;L(Θi | Y)為觀測似然值;L(Θi)為先驗似然值;C為歸一化因 子;Y為預(yù)測因變量(即預(yù)測結(jié)果）：9,為第i組參數(shù)組； 當(dāng)再有新的實測數(shù)據(jù)獲取時，前一次更新過的后驗似然值結(jié)果即變?yōu)楫?dāng)前計算的先驗 似然值結(jié)果，重復(fù)上述步驟，即可求得更新的預(yù)測結(jié)果。
【文檔編號】G06F19/00GK105930670SQ201610279685
【公開日】2016年9月7日
【申請日】2016年4月29日
【發(fā)明人】侯迪波, 許樂, 劉勛, 王柯, 劉景明, 黃平捷, 張光新, 張宏建
【申請人】浙江大學(xué)