本發(fā)明涉及河口地區(qū)污染物溯源領域，具體提出一種專用于河口地區(qū)突發(fā)性水污染事件點源溯源方法。

背景技術：

海洋運輸是國際間商品交換中最重要的運輸方式之一，其貨物運輸量占全部國際貨物運輸量的比例大約在80％以上。由于我國東部和南部大陸海岸線有1.8萬多千米，內(nèi)海和邊海的水域面積約為470多萬平方千米，沿海城市較多，再加上近年來我國經(jīng)濟的騰飛以及綜合實力的不斷提升，我國的進出口貨物運輸總量的80％～90％是通過海洋運輸進行的。隨著國際航運業(yè)的發(fā)展，超巨型油輪已達60多萬噸，第五代集裝箱船的載箱能力已超過5000teu。但是隨著我國社會生產(chǎn)力和科學技術的迅猛發(fā)展，沿海城市工業(yè)生產(chǎn)與海上運輸業(yè)務更加繁忙，使海洋受到了不同程度的污染和破壞。

海岸河口其特殊的地理位置，由于受到潮汐作用影響，使得無論來自灣內(nèi)或外海水體中的污染物都會影響到該河口及其附近區(qū)域的水質。目前對海岸河口地區(qū)污染物溯源問題研究尚不成熟，其較河道有更為復雜的流場，缺乏對污染源位置及強度進行識別的一種簡單、可靠方法。

技術實現(xiàn)要素：

本發(fā)明是為了解決河口地區(qū)溯源問題而進行的，目的是針對河口地區(qū)發(fā)生突發(fā)性水污染事件后能夠迅速實施應急措施的需要，提供一種快速、精確的污染物溯源方法對點污染源的位置及強度進行識別。

本發(fā)明的上述技術目的是通過以下技術方案得以實現(xiàn)的：

一種河口點源突發(fā)性水污染事件溯源方法，包括以下步驟：

步驟1，選擇任意一個河口地區(qū)作為模型模擬區(qū)域并選擇模型模擬區(qū)域內(nèi)的任一位置作為污染物濃度觀測站，通過觀測站對模型模擬區(qū)域的流場進行監(jiān)測，建立實時河口流場數(shù)據(jù)庫；

步驟2，從觀測站監(jiān)測到污染物濃度開始，每隔一段時間監(jiān)測一次觀測點的污染物濃度過程數(shù)據(jù)并記錄；

步驟3，根據(jù)觀測站監(jiān)測時段內(nèi)的污染物濃度過程數(shù)據(jù)，結合河口地區(qū)半日潮特性，計算出污染物初始投放時間；

步驟4，根據(jù)已有流場數(shù)據(jù)庫的數(shù)據(jù)及的污染物初始投放時間，計算出污染物投放位置；

步驟5，建立污染物遷移模型，求解出非恒定流條件下的物質濃度；

步驟6，將污染源位置和強度解耦；

步驟7，利用相關系數(shù)構建污染物溯源位置優(yōu)化模型的目標函數(shù)，并通過步驟4中計算出的污染源位置計算目標函數(shù)值；

步驟8，對傳統(tǒng)遺傳算法按照改進初始種群、變異和選擇三個步驟進行改進，形成改進遺傳算法；

步驟9，構建污染物溯源強度優(yōu)化模型，通過位置優(yōu)化模型求解結果得到污染源位置；

步驟10，利用改進遺傳算法對污染物溯源強度優(yōu)化模型進行求解。

綜上所述，本發(fā)明具有以下有益效果：本發(fā)明針對海岸河口地區(qū)岸邊、來往船舶等污水未經(jīng)處理排放的現(xiàn)象及海洋開采、運輸過程造成的泄漏事故，研究了河口地區(qū)潮流特點，以水文水動力計算為基礎，根據(jù)當前水文情景及實測污染物濃度變化，提出了一套適用于河口地區(qū)突發(fā)性水污染事故的溯源模型，并且改進了遺傳算法對模型的求解，從而快速、準確地確定污染源位置及強度，對環(huán)保部門短時間內(nèi)掌握污染源信息并對污染事件進行應急處理提供了重大的指導意義。

附圖說明

圖1為河口點源突發(fā)性水污染事件溯源方法的工作流程圖；

圖2為改進遺傳算法計算步驟；

圖3為實施例目標區(qū)域、污染源和觀測點位置、污染源位置先驗范圍及初始種群示意圖；

圖4為觀測點污染物濃度觀測數(shù)據(jù)；

圖5為河口污染物運動示意圖；

圖6為傳統(tǒng)遺傳算法與改進遺傳算法計算比較圖；

圖7為污染源強度計算結果示意圖。

具體實施方式

本發(fā)明是江蘇省研究生教育教學改革課題中的一個重要研究項目。以下結合附圖對本發(fā)明涉及的河口突發(fā)性水污染事件點源溯源方法進行詳細地說明，但該實施例不應理解為對本發(fā)明的限制。

在本實施例中，目標區(qū)域選擇于福建泉州灣，該地區(qū)地形如附圖3所示，位置坐標采用平面坐標，污染源位置位于(367463，2751257)，在該污染源位置處傾倒質量為36000kg的污染物，在該案例中為了簡化計算，忽略物質輸移過程中的擴散項。在觀測點位置(372072，2749507)上午9點監(jiān)測到污染物開始，每隔半小時觀測一次，連續(xù)觀測24h，觀測得到的濃度數(shù)據(jù)如附圖4。污染源和觀測位置分別在圖3中以位于下方的十字和空心圓標出。

具體地，本實施例的河口突發(fā)性水污染事件點源溯源方法包括以下步驟：

步驟1，選擇任意一個河口地區(qū)作為監(jiān)測區(qū)域并在監(jiān)測區(qū)域建立測站，通過測站對監(jiān)測區(qū)域的流場進行監(jiān)測，建立實時河口流場數(shù)據(jù)庫；

步驟2，選擇該河口地區(qū)內(nèi)的任一位置作為觀測點，從觀測點上午9點監(jiān)測到污染物濃度開始，每隔0.5h監(jiān)測一次該地點污染物濃度并記錄，持續(xù)觀測24h，污染物濃度隨時間變化的坐標系如圖4；

步驟3，圖5為污染物在投放點及觀測點之間運動示意圖。假設污染物從a點投入，觀測點為c點，b點為污染物到達的最遠位置。由于受到潮汐潮流漲落影響，污染物從投入點隨落潮流運動到b點后會隨漲潮流重新回到a點。圖5中，時間t1為污染物從a點投入后隨落潮流運動到c點的時間，t2為污染物從a點隨落潮流運動到b點又后隨漲潮流回到c點的時間。由圖5，并且結合半日潮流性質可得下式：

一般情況下，根據(jù)所測得的數(shù)據(jù)，t1、t2的實際數(shù)值是未知的。圖4，令第一個峰值時間為t1，第二個峰值時間為t2，則t2-t1是可知的，大約為6h，由式(1)可得出t1的具體數(shù)值約為3h，其與圖2中第一個峰值所在時間點一致，即從污染開始排放到第一次運動至觀測點位置時的時間為3h。這說明可以通過觀測點一系列的觀測數(shù)據(jù)以及半日潮雙峰特性來計算出污染物初始釋放時間，該實施例中為上午6點。

步驟4，根據(jù)已有流場數(shù)據(jù)及步驟3中的污染物初始投放時間，利用公式(2)、(3)，可以計算出污染源投放位置；

其中：δt為流場測量時間間隔；i為時間間隔個數(shù)；u(i)[lt^-1]和v(i)[lt^-1]分別為x、y方向上速度；x(0)、y(0)和x(t)、y(t)分別為時間t＝0和t＝t時質點位置。

在該實施例中δt取為100s，(x(t)，y(t))為觀測點位置(372072，2749507)，則反算出的(x(0)、y(0))為(367682，2752337)，其位置如圖3中位于下方的十字形。

步驟5，根據(jù)下式建立污染物遷移模型；

式(4)中：h為水深[l]；c為等深的物質濃度[ml^-3]；t為時間[t]；u[lt^-1]和v[lt^-1]分別為x、y方向上速度；kxx，kxy，kyx以及kyy為二維擴散系數(shù)張量[l²t^-1]。

對于式(4)，采用有限體積法進行差分求解，利用算子分裂方法，將上式分成對流項和擴散項兩部分進行分步計算，可以求解出非恒定流條件下的污染物濃度。

步驟6，將污染源位置和強度解耦。由污染物遷移模型，根據(jù)步驟3推導出的污染物初始排放時間，結合式(5)相關系數(shù)r，假設污染源位置為(x′0，y′0)，由于污染源強度值對式(5)相關系數(shù)無影響，令污染源強度為一定值m′0：

式(5)中：cj(j＝1，2……k)，(k≤n)為k個時刻污染物濃度測量數(shù)據(jù)，cj′(j＝1，2……k)，(k≤n)為k個時刻污染物濃度計算數(shù)據(jù)。和分別為數(shù)據(jù)系列cj和cj′的算術平均值。

步驟7，構建污染物溯源位置優(yōu)化模型的目標函數(shù)，如式：

f(x，y)＝min(abs(1-r))(6)

參數(shù)x，y的先驗范圍可以根據(jù)步驟4中計算出的x(0)、y(0)值給定，如式：

x0min≤x≤x0max；y0min≤y≤y0max(7)

式(7)中，x0min＝364890，x0max＝369304，y0min＝2749569，y0max＝2752827。該范圍如圖3中矩形框所示。根據(jù)步驟6假定的污染源位置(x0′，y0′)，計算出目標函數(shù)值，目標函數(shù)值越小，則該假定的污染源位置越接近于真實源位置。

步驟8，對傳統(tǒng)遺傳算法進行改進，得到改進遺傳算法，改進后遺傳算法步驟如下：

步驟8a、改進初始種群。初始種群的生成帶有一定的隨機性，為了加快收斂速度，我們對初始種群進行優(yōu)選。初始種群規(guī)模為n，按照傳統(tǒng)遺傳算法的規(guī)則產(chǎn)生2n的個體，并將2n的個體按照適應度函數(shù)的值由小到大排序，選取前n的個體作為初始種群，中i為種群個體編號并且為整數(shù)，gen為進化代數(shù)。

步驟8b、變異。變異帶有隨機性，父代的優(yōu)良性能經(jīng)過變異以后可能并未遺傳給子代。為了挑選出適應性能更好的變異個體，改進后的算法產(chǎn)生變異個體的方式做如下修改。首先，按照傳統(tǒng)遺傳算法產(chǎn)生中間變異個體vi′^gen，將其按適應度函數(shù)f(vi′^gen)的值由小到大排序。變異個體產(chǎn)生方式如下，下式中α為權重系數(shù)：

步驟8c、交叉。按照傳統(tǒng)遺傳算法產(chǎn)生中間變異個體ui′^gen，將其按適應度函數(shù)f(ui′^gen)的值由小到大排序。交叉?zhèn)€體產(chǎn)生方式如下：

步驟8d、選擇。按照傳統(tǒng)遺傳算法產(chǎn)生中間選擇個體xi′^gen+1，將其按適應度函數(shù)f(xi′^gen+1)的值由小到大排序。選擇個體產(chǎn)生方式如下：

在步驟7中的先驗信息范圍內(nèi)生成初始種群，其具體分布如圖3中若干菱形所示。分別用遺傳算法迭代5次和10次，與傳統(tǒng)遺傳算法進行比較，結果如圖6，(a)、(b)分別為傳統(tǒng)遺傳算法迭代5次和10次后結果，(a)、(b)中叉形表示其迭代后種群；(c)、(d)分別為改進后遺傳算法迭代5次和10次后結果，(c)、(d)中叉形表示其迭代后種群。

步驟9，構建污染物溯源強度優(yōu)化模型。通過位置優(yōu)化模型求解結果得到污染源位置，固定該位置不變，假定眾多m′0，分別用式(4)污染物遷移模型計算c′1，c′2，c′3……c′k。其先驗范圍可以由下式給出：

根據(jù)上式可以構建優(yōu)化模型進一步推算污染源強度。

目標函數(shù)：m′0＝min(∑ωi(cj′-cj)²)(12)

目標函數(shù)：m0min≤m′0≤m0max(13)

式中：m′0的約束范圍可以依據(jù)式(7)進行一定比例放大縮小后確定；系數(shù)ωi＝1/(cj+1.0)²。

步驟10，利用改進遺傳算法對污染物溯源強度優(yōu)化模型進行求解，為了減小模型隨機性對結果的影響，重復獨立計算20次，求解結果如圖7。

根據(jù)所計算出的污染物源位置與強度，可以給相關部門快速提供可靠的污染源信息，便于環(huán)保主管部門針對該污染源快速制定相應的應急措施，最大限度地減小污染所造成的損失。

本具體實施例僅僅是對本發(fā)明的解釋，其并不是對本發(fā)明的限制，本領域技術人員在閱讀完本說明書后可以根據(jù)需要對本實施例做出沒有創(chuàng)造性貢獻的修改，但只要在本發(fā)明的權利要求范圍內(nèi)都受到專利法的保護。