張謙，裴海龍，鄔依林

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

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

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

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

1 濃度邊界跟蹤方案

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

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

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

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

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

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

如圖2 所示為簡(jiǎn)化的模擬岸邊污染源擴(kuò)散濃度場(chǎng)情況，圖2（a）為濃度場(chǎng)擴(kuò)散三維視圖，圖2（b）為濃度場(chǎng)等高線，仿真參數(shù)如表1 所示，這里假設(shè)單位時(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 濃度場(chǎng)模擬Fig.2 Concentration field simulation

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

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

2.2 濃度場(chǎng)重建

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

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

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

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

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

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

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

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

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

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

3 濃度邊界追蹤

3.1 軌跡跟蹤

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

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

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

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

3.2 邊界追蹤仿真

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

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

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

4 結(jié) 論

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