結(jié)合空間插值的水污染濃度邊界跟蹤*

張謙，裴海龍，鄔依林

1. 廣東第二師范學(xué)院計算機(jī)科學(xué)系，廣東 廣州 510303

2. 華南理工大學(xué)/自主系統(tǒng)與網(wǎng)絡(luò)控制教育部重點實驗室/廣東省無人機(jī)系統(tǒng)工程技術(shù)研究中心，廣東 廣州 510640

水體污染是水資源保護(hù)的主要問題之一，在常規(guī)水環(huán)境保護(hù)上，各國在廢水排放、水質(zhì)監(jiān)測以及質(zhì)量標(biāo)準(zhǔn)等方面均已建立其行之有效的水環(huán)境監(jiān)控保護(hù)體系，可以有效保護(hù)水資源［1］。然而隨著經(jīng)濟(jì)社會的發(fā)展，工業(yè)排放、水上作業(yè)、農(nóng)業(yè)生產(chǎn)等活動，均有可能造成一些不可預(yù)見的突發(fā)性水污染事故，進(jìn)而對污染監(jiān)控提出了更高的要求，需要快速、準(zhǔn)確以及廉價的水體污染監(jiān)控方法［2］。一般而言，對于水體污染監(jiān)測，可以采用衛(wèi)星遙感［2-4］、傳感器網(wǎng)絡(luò)、數(shù)值模擬等方式進(jìn)行監(jiān)測。然而，對于突發(fā)性水污染事件，一般具有局部性的特點，而針對突發(fā)性水污染擴(kuò)散的數(shù)值模擬，是水污染事故應(yīng)急處理中的重要手段［5］?？紤]到在不同水體環(huán)境下水污染濃度場的數(shù)值模擬條件，文獻(xiàn)［6］模擬了急流條件下突發(fā)水污染的污染物二維水流-輸運(yùn)模型；文獻(xiàn)［7］結(jié)合深度平均流速的橫向指數(shù)分布和橫向擴(kuò)散系數(shù)的二維變化關(guān)系式，提出岸邊排放污染混合區(qū)的形狀特征與分類；文獻(xiàn)［8］提出了傾斜岸水面污染源下角形域污染物二維擴(kuò)散濃度分布；這些模型為水污染控制與監(jiān)測以及應(yīng)急處理中的快速決策提供了很好的依據(jù)。但在實際的水污染突發(fā)事件的應(yīng)急處置中，更關(guān)注的是污染物擴(kuò)散的邊界、濃度分布及其運(yùn)動過程［8］。不同于常規(guī)的水環(huán)境監(jiān)測，水利部門可以通過固定觀測點構(gòu)建監(jiān)測網(wǎng)絡(luò)，然而對于突發(fā)性水污染區(qū)域的監(jiān)測，因為污染源位置的不確定性，更方便的是使用移動設(shè)備對污染區(qū)域進(jìn)行監(jiān)測，建立污染區(qū)域的濃度場分布，進(jìn)而確定濃度邊界。然而，考慮到地形等因素的影響，很難使用移動監(jiān)測設(shè)備對突發(fā)性水污染擴(kuò)散濃度場進(jìn)行完全的實時測量，因此依靠移動采樣設(shè)備獲取測量樣本，采用空間插值重建濃度場是一種有效的手段。地統(tǒng)計空間插值方可以估計任何地點的數(shù)據(jù)，其中克里金（Kriging）插值法是最常用的空間插值法［9-11］。如果能夠通過移動檢測設(shè)備采樣并使用克里金插值法，實現(xiàn)水污染擴(kuò)散區(qū)域的濃度場重構(gòu)，那么將為移動機(jī)器人進(jìn)行水污染濃度場邊界的確立提供便利。濃度場邊界追蹤屬于機(jī)器人軌跡跟蹤問題，在這方面很多學(xué)者使用模型預(yù)測控制方法開展研究［12-16］，例如文獻(xiàn)［12］使用模型預(yù)測控制方法設(shè)計了一個水下機(jī)器人軌跡跟蹤系統(tǒng)；文獻(xiàn)［13］考慮實際工作約束和限制，提出了一種基于模型預(yù)測控制輪式機(jī)器人軌跡規(guī)劃方法；文獻(xiàn)［14］針對水下動態(tài)目標(biāo)跟蹤，提出了基于模型預(yù)測控制的水下機(jī)器人跟蹤控制方法等。然而，這些跟蹤的方法，一般都是由狀態(tài)方程直接收斂到目標(biāo)軌線。在水污染濃度邊界跟蹤上，因為無法事先獲取目標(biāo)點的位置，必須依靠移動機(jī)器人在移動過程中的采樣，來預(yù)測目標(biāo)點的位置，同時決定移動機(jī)器人的移動方向，達(dá)到跟蹤的目的。

因此，針對上述問題，提出了使用模型預(yù)測控制與空間插值算法相結(jié)合，進(jìn)行水污染擴(kuò)散濃度邊界跟蹤的方法，通過空間插值重構(gòu)的水污染擴(kuò)散濃度場分布結(jié)合系統(tǒng)的狀態(tài)模型，再加入符合設(shè)定環(huán)境的控制量約束，尋找局部最優(yōu)點作為移動機(jī)器人跟蹤的濃度場邊界最優(yōu)軌跡，從而確立污染區(qū)域濃度場邊界。

1 濃度邊界跟蹤方案

如圖1 所示框架，為結(jié)合空間插值與模型預(yù)測控制的濃度場邊界跟蹤方案。為了驗證方案的可行性，首先對水污染擴(kuò)散區(qū)域的濃度場進(jìn)行了數(shù)值模擬，并將數(shù)值模擬的結(jié)果作為假設(shè)真實發(fā)生水污染的濃度場，然后通過對真實濃度場的隨機(jī)采樣，對比不同條件下的空間插值結(jié)果，選擇合適的參數(shù)擬合結(jié)果進(jìn)行污染區(qū)域濃度場的重構(gòu)。最后根據(jù)重構(gòu)的濃度場以及擬合參數(shù)，結(jié)合水面移動機(jī)器人狀態(tài)模型，對移動機(jī)器人的移動位置濃度值進(jìn)行預(yù)測，并以模型預(yù)測控制方法對移動機(jī)器人所在區(qū)域進(jìn)行局部尋優(yōu)，尋找并確定下一步的移動目標(biāo)，以此迭代實現(xiàn)對濃度場邊界的跟蹤并最終確立污染擴(kuò)散區(qū)域的濃度邊界。

圖1 濃度邊界跟蹤方案Fig.1 Scheme of boundary tracking for water pollution concentration field

2 污染區(qū)域濃度場分布與重構(gòu)

2.1 污染源擴(kuò)散模型

污染物在水體環(huán)境下的濃度擴(kuò)散，是一個非常復(fù)雜的過程，需要綜合考慮水流速度以及擴(kuò)散系數(shù)變化因素。近年來，國內(nèi)外很多學(xué)者對污染物擴(kuò)散濃度場進(jìn)行了數(shù)值模擬［6-8］，其中文獻(xiàn)［7］提出考慮河流流速和橫向擴(kuò)散系數(shù)變化的河流岸邊排放污染混合區(qū)的估算方法，描述了河流岸邊排放污染區(qū)的形狀特征，較好地表征了排污等值線形狀，對水污染排放區(qū)域的數(shù)值描述有較強(qiáng)的代表性。本文主要研究的是對污染源在水面擴(kuò)散時的濃度場邊界追蹤問題，濃度場數(shù)值模型［7］為

式（1）中x為自排污口沿河流流向的縱向坐標(biāo)；y為垂直于x軸從排污口指向河心的橫向坐標(biāo)，u1表示離岸橫坐標(biāo)y1的流速，M?表示岸邊污染源單位時間的排放質(zhì)量，H表示河道平均水深，m和n為正常數(shù)指數(shù)，γ和α為正常數(shù)，p= 2 +m-n＞1，φ=(1+m)/p，Γ(φ)為完全伽馬函數(shù)。

如圖2 所示為簡化的模擬岸邊污染源擴(kuò)散濃度場情況，圖2（a）為濃度場擴(kuò)散三維視圖，圖2（b）為濃度場等高線，仿真參數(shù)如表1 所示，這里假設(shè)單位時間排放量M? = 1000 kg/s，平均水深H= 5 m，河流流速為常數(shù)0.1 m/s，橫向分布指數(shù)m和n都為0，p= 2+m-n，γ= 5.96，α= 0.01，φ=(1+m)/p，Γ(φ)= 1.

圖2 濃度場模擬Fig.2 Concentration field simulation

表1 水污染濃度場數(shù)值模擬參數(shù)Table1 Parameters for concentration field simulation

從圖2（b）可以看出，在近排放點位置的地方，濃度值比較大，等濃度線也比較密集，隨著位置逐漸變化，濃度場變化符合文獻(xiàn)［7］中所歸納的岸邊排放污染源混合區(qū)標(biāo)準(zhǔn)型類別，其形狀特征近似于半橢圓形。

2.2 濃度場重建

通過式（1）可以模擬出污染源擴(kuò)散的濃度場情況，但在實際應(yīng)用過程中，采用的是單一移動采樣設(shè)備對污染源的濃度場邊界進(jìn)行追蹤，也很難對濃度場分布進(jìn)行建模后得到精確的分布情況后再進(jìn)行追蹤，因此，更直接的方式是通過采樣設(shè)備對濃度場分布范圍內(nèi)進(jìn)行隨機(jī)濃度采樣，然后通過地統(tǒng)計學(xué)方式，對采樣得到的濃度值進(jìn)行空間插值，重構(gòu)污染區(qū)域的濃度場。在地統(tǒng)計學(xué)中，最常用的方法是Kriging 插值法。Kriging 方法從變量相關(guān)性和變異性出發(fā)，在有限區(qū)域內(nèi)對區(qū)域化變量的取值進(jìn)行無偏、最優(yōu)估計的一種方法。很多文獻(xiàn)對Kriging插值法進(jìn)行了詳細(xì)闡述［9-11］，其基本算法思路如下所述。

假設(shè)x是濃度場的任一點，Z(x)是該點的測量值，如果在整個濃度場內(nèi)總共有n個觀測點x1，x2，…，xn，那么對于任意待估計點的實測濃度值Zv(x)，其估計的濃度值Z*v(x)，可以通過該待估點影響范圍內(nèi)的n個有效樣本值Zv(xi)(i= 1，2，…，n)的線性組合來表示，即

為了擬合觀測數(shù)據(jù)值與空間距離之間的關(guān)系，一般變異函數(shù)有多種模型，常見的有高斯模型、指數(shù)模型、球狀模型等。例如，高斯模型可以表示為

綜上所述，在進(jìn)行污染源擴(kuò)散的濃度場重建的時候，就是通過采樣得到的觀測值，進(jìn)行變異函數(shù)擬合，再求解式（8）得到權(quán)重系數(shù)λi，即可對不同位置點進(jìn)行估計，從而最終重建出整個濃度場。

為了驗證使用Kriging 插值重構(gòu)岸邊污染源擴(kuò)散濃度場的效果，圖3 以不同的隨機(jī)樣本數(shù)擬合球形變差函數(shù)模型，對圖2模擬的濃度場進(jìn)行了插值測仿真，以不同的采樣樣本數(shù)，擬合球形變差函數(shù)的結(jié)果如表2所示，仿真結(jié)果如圖3所示。

表2 擬合球形變差函數(shù)結(jié)果Table 2 Fitness statistics of spherical variation function models

圖3 不同樣本數(shù)擬合球形變差函數(shù)的濃度場重構(gòu)Fig.3 Reconstruction of concentration field by fitting spherical variation function with different sample numbers

如圖3 所示，圖3（a）使用100 個樣本、圖3（b）使用200 個樣本、圖3（c）使用500 個樣本，分別對圖2所示模擬的濃度場進(jìn)行Kriging 插值所得的濃度場預(yù)測等高線、變差函數(shù)擬合效果以及預(yù)測結(jié)果與實際分布之間的誤差，從表2 以及圖3 可知，樣本數(shù)越多，變差函數(shù)的擬合效果越好，對濃度場的重建效果越好，此外，預(yù)測誤差分布主要集中在離污染源排放位置較近的區(qū)域，越遠(yuǎn)離排放源，誤差越接近于0。因為本文主要是對污染源擴(kuò)散邊界進(jìn)行跟蹤，跟蹤的軌線離排放源較遠(yuǎn)，因此經(jīng)過插值重構(gòu)的污染源濃度場邊界，可以反映真實的濃度擴(kuò)散，是可以用來作為追蹤軌線的。

圖4為采用500個采樣點時，變差函數(shù)分別擬合高斯模型和指數(shù)模型時重建濃度場分布的效果。從圖4（a）擬合高斯模型的效果來看，重構(gòu)的濃度場邊界不能完成呈現(xiàn)半橢圓形狀，且重建的濃度場與原濃度場之間的誤差有較大的波動，尤其是離排放源較遠(yuǎn)的地方，誤差變動仍然較大，不利于后續(xù)開展邊界追蹤操作。圖4（b）為擬合指數(shù)模型效果，從重構(gòu)的濃度場分布情況看，等濃度線重構(gòu)效果優(yōu)于圖3（b）采用200 采樣點的球狀模型，但在x方向較小時，y方向曲線效果不如圖3（c），誤差波動比較大，且變差函數(shù)擬合的效果不理想。因此，綜上分析，在進(jìn)行邊界追蹤重建濃度場分布的時候，將采用500個采樣點擬合球狀變差函數(shù)模型進(jìn)行空間插值。

圖4 變差函數(shù)擬合高斯模型和指數(shù)模型的濃度場重構(gòu)Fig.4 Reconstruction of concentration field by fitting Gaussian and exponential variation function

3 濃度邊界追蹤

3.1 軌跡跟蹤

對污染物擴(kuò)散濃度邊界進(jìn)行追蹤，事實上是控制移動機(jī)器人的位置和方向，使得移動機(jī)器人自起始位置開始，其移動軌跡快速靠近濃度邊界等濃度線，并沿著等濃度線進(jìn)行移動。假設(shè)移動機(jī)器人是一個在水面移動的采樣設(shè)備，那么它的3自由度動態(tài)模型如下式所示［13-14］

式（10）中，狀態(tài)向量[x，y，θ]T分別表示位置和移動的方向角，v表示移動的速度，可以設(shè)為常數(shù)。移動機(jī)器人在進(jìn)行濃度邊界追蹤時，移動機(jī)器人首先要根據(jù)隨機(jī)采樣得到樣本進(jìn)行當(dāng)前位置的Kriging 插值濃度預(yù)測，根據(jù)預(yù)測值與重構(gòu)的邊界濃度值進(jìn)行比較，尋找局部最優(yōu)點，從而控制移動機(jī)器人的移動位置和方向，如此循環(huán)操作達(dá)到邊界最終的目的，這個過程與模型預(yù)測控制進(jìn)行移動機(jī)器人軌跡跟蹤相匹配。因此，將以模型預(yù)測控制的方式對濃度邊界進(jìn)行追蹤。

在模型預(yù)測控制中，對于連續(xù)的系統(tǒng)狀態(tài)方程，通過選取合適的采樣頻率，將系統(tǒng)離散化為

式（17）中，Np和Nu分別表示預(yù)測步長與控制步長。K(x(k+i)，y(k+i))表示根據(jù)移動機(jī)器人預(yù)測所在位置，結(jié)合擬合參數(shù)進(jìn)行Kriging 插值所得到的預(yù)測位置的污染濃度值，cr表示追蹤的參考濃度值，Δω為控制量的變化，Q= diag[q1，q2，…，qNp]和R表示懲罰權(quán)重。通過最小化目標(biāo)函數(shù)J，就可以得到控制量ω。

3.2 邊界追蹤仿真

根據(jù)式（10）和式（17），結(jié)合第二節(jié)所述濃度場重構(gòu)，水面移動機(jī)器人跟蹤濃度邊界進(jìn)行跟蹤仿真。由式（17），假設(shè)參考軌跡為污染源濃度場分布中濃度值為50 的等濃度線，軌跡預(yù)測值為當(dāng)前位置下，通過500個采樣點擬合球形變差函數(shù)后的Kriging插值預(yù)測值，預(yù)測步長和控制步長均設(shè)置為3，系統(tǒng)的控制量u=ω，移動速度設(shè)置為常數(shù)，僅對角速度變化進(jìn)行約束，約束條件為|ω|＜π/2，也就是角度θ=ωt不能在單位時間內(nèi)進(jìn)行大方向的轉(zhuǎn)動。

如圖5 所示，圖5（a）所模擬的真實的濃度場分布數(shù)值設(shè)置參數(shù)與圖2（b）一致，圖中黑色的點為在濃度場范圍內(nèi)隨機(jī)采樣點的位置。圖5（b）是變差函數(shù)擬合球狀模型的結(jié)果，從圖上可知，基本能夠完全擬合。圖5（c）是通過采樣模擬的真實濃度場，進(jìn)行空間插值而重構(gòu)的濃度場，這個重構(gòu)的濃度場就是將要進(jìn)行軌跡跟蹤的邊界。圖5（d）即為濃度邊界的軌跡跟蹤結(jié)果，紅色點表示移動機(jī)器人從起始位置向指定濃度值所在區(qū)域靠攏時的每一次計算所得到的移動軌跡，紅色的星形圖標(biāo)表示移動機(jī)器人的位置移動到了預(yù)定濃度邊界，并將從該位置開始進(jìn)行等濃度線邊界追蹤，當(dāng)?shù)竭_(dá)岸邊時跟蹤停止。圖中假設(shè)移動機(jī)器人起始位置為任意設(shè)置，并將跟蹤濃度值為50 的等濃度線，移動機(jī)器人的速度v= 1m/s 設(shè)置為常數(shù)。對于移動機(jī)器人的每一次迭代位置，均需要通過Kriging 插值計算預(yù)測步長決定的移動機(jī)器人所在位置的濃度值，并根據(jù)預(yù)測的濃度值與等濃度線的差，以及控制量增量，對最小化式（17）的目標(biāo)函數(shù)求最優(yōu)值，從而得到所需的最優(yōu)控制量ω，以此控制量修正移動機(jī)器人的位置，當(dāng)移動到下一個位置時，再重復(fù)上述預(yù)測求最優(yōu)的步驟，從而最終實現(xiàn)水面移動機(jī)器人跟蹤污染物擴(kuò)散濃度邊界。從圖5（d）可以看到，濃度邊界跟蹤的效果非常好，能沿著等濃度線進(jìn)行移動，這也驗證了通過空間插值預(yù)測的模型預(yù)測控制濃度邊界軌跡跟蹤方法可行。

圖5 濃度邊界追蹤Fig.5 Boundary tracking for water pollution concentration field

4 結(jié) 論

本文通過數(shù)值模擬河流岸邊排放污染區(qū)濃度擴(kuò)散情況，并比較不同的Kriging 空間插值變差函數(shù)擬合方法，重構(gòu)污染區(qū)域的濃度場分布，通過模型預(yù)測控制方式，實現(xiàn)使用水面移動機(jī)器人設(shè)備跟蹤濃度擴(kuò)散邊界軌跡的目的。仿真實驗結(jié)果表明，該方法能夠穩(wěn)定、快速、準(zhǔn)確地進(jìn)行濃度邊界跟蹤，可以為實際應(yīng)用過程中的區(qū)域跟蹤提供有價值的借鑒。