基于徑向基神經(jīng)網(wǎng)絡的水下自主航行器尋源算法研究

吳宗秀，吳 超

(1. 上海交通大學 船舶海洋與建筑工程學院，上海 200240; 2. 上海交通大學 海洋工程國家重點實驗室，上海 200240)

水下尋源著眼于解決利用單個或編隊水下航行器對未知分布的信號場進行信號源搜索問題。水下工程領域的許多問題都可以歸結為水下尋源問題。包括：失事飛行器黑匣子的尋找、天然氣或者水下礦藏的定位、熱液噴口搜索[1]等。在這些問題中，未知信號場則相應地可以具體化為聲學信號場、化學物質(zhì)濃度信號場及溫度信號場等。水下尋源問題的解決方案對于擴展水下航行器的應用而言有著重大的意義，該領域的解決方法一般有瞬態(tài)決策方法、極值搜索算法、高斯過程回歸以及神經(jīng)網(wǎng)絡等。

Holland等[2]、Russell[3]分別就未知信號場的尋源問題提出了六邊形和Z字型前進路線的搜索策略，通過離散的信號測量以及簡單的信號強度對比進行搜索路徑的規(guī)劃。Russell[4]后續(xù)將工作拓展到三維場源搜索問題中，瞬態(tài)極值策略可以作為尋源問題一種相對穩(wěn)定的解決方案,但是其決策得出的方向精度不足。自動控制領域的方法在場源搜索問題中亦得到了廣泛的應用，Mellucci等[5]、Matveev等[6-7]均運用了滑模控制方法來控制水下航行器的首向進行信號場的極值搜索，其結果顯示，滑模控制方法生成的搜錄路徑相對較長。Cochran等[8]將極值搜索算法(extreme value search algorithm, 簡稱ESA)運用到水下尋源問題中，通過控制航行器的轉(zhuǎn)向引導航行器進行搜索，該團隊在其他研究[9]中增加前進速度為控制參數(shù)并且得到了更好的結果，但同時也指出該算法收斂速度較慢的問題。

高斯過程回歸(gaussian process regression, 簡稱GPR)具有良好的函數(shù)逼近能力，常用于解決未知信號場的擬合問題。Ai等[10]結合高斯過程回歸、Armijo步長準則、遞歸最小二乘方法和非線性模型預測控制(nonlinear model predictive control, 簡稱NMPC)進行完備的水下航行器場源搜索仿真。其研究結果表明航行器在搜索過程前期受先驗信號場影響明顯，在接近先驗信號場的最大值區(qū)域后轉(zhuǎn)向正確的路徑并收斂到真實場源。閻述學等[11]的研究著眼于利用高斯過程回歸進行水下信號場的熱點采樣，通過定義一種與均方差相關的引力模型引導航行器進行多熱點未知信號場的自適應采樣。相較于傳統(tǒng)的割草機和梳型路徑，該算法能夠以較短的路徑實現(xiàn)更準確的擬合。

徑向基神經(jīng)網(wǎng)絡(radial basis function neural network，簡稱RBFNN)具備無限精度的函數(shù)逼近能力，能夠運用于信號場的估計并且具有結構簡單、易于理解的優(yōu)點。Jeon等[12]通過推導連續(xù)時間形式的遞歸最小二乘法對RBFNN進行在線訓練，并控制裝有多個傳感器的雙連桿機械手進行二維灰度圖的極值尋找。然而利用RBFNN的尋源算法研究很少考慮神經(jīng)網(wǎng)絡的正則化，實際上正則化對于神經(jīng)網(wǎng)絡的穩(wěn)定性和泛化性非常重要。

以水下自主航行器(autonomous underwater vehicle, 簡稱AUV)在熱泉區(qū)域的硫化氫濃度分布場場源搜索為背景，使用具有增長式結構的RBFNN對未知信號場進行在線無先驗漸進擬合，并通過增量式奇異值分解算法加速神經(jīng)網(wǎng)絡的正則化計算，此外還利用動量梯度法解決AUV路徑容易陷入局部最優(yōu)區(qū)域的問題。算法的有效性通過單峰值和多峰值硫化氫濃度場中的場源搜索模擬計算進行驗證。

1 未知場在線估計

1.1 水下尋源問題描述

1.2 RBFNN介紹

徑向基函數(shù)神經(jīng)網(wǎng)絡(RBFNN)，具有結構簡單、學習速度快、擬合精度高、泛化能力較強和不易陷入?yún)?shù)局部最優(yōu)等優(yōu)點，被廣泛應用于函數(shù)逼近、分類、時間序列預測等諸多領域，其網(wǎng)絡結構示意如圖1所示。

圖1 RBFNN結構示意Fig. 1 View of RBFNN’s construction

RBFNN在結構上是嚴格的三層網(wǎng)絡，分別為輸入層、隱藏層和輸出層。輸入層負責信號傳遞，隱藏層通過徑向基函數(shù)(radial basis function, 簡稱RBF)進行從輸入層到隱藏層空間的非線性變換，輸出層對隱藏層的響應進行加權求和。

選用高斯型徑向基函數(shù)：

(1)

其中，ωj為權重系數(shù)，cj為核函數(shù)中心，rj為高斯函數(shù)半徑。

設計矩陣H，其元素Hi,j=hj(xi)，行數(shù)記作p，列數(shù)記作m：

(2)

權重向量W可表示為：

(3)

(4)

權重向量可以通過計算偽逆或最小二乘法求取。為了避免矩陣求解中的病態(tài)問題以及神經(jīng)網(wǎng)絡的過擬合問題，引入全局正則化操作，對應權重向量求解公式如下：

(5)

其中，參數(shù)λ為全局正則化系數(shù)，該參數(shù)能夠提高神經(jīng)網(wǎng)絡的泛化性和計算穩(wěn)定性，λ的數(shù)值需要在網(wǎng)絡結構和數(shù)據(jù)集變化的時候進行優(yōu)化。

綜上，RBFNN對任意輸入x∈R2的信號值預測為：

(6)

1.3 RBFNN在線訓練

徑向基神經(jīng)網(wǎng)絡的預測性能取決于RBF中心{C1C2……Cm}、權重向量W、全局正則化參數(shù)λ等。為了實現(xiàn)RBFNN的在線訓練，算法應當能夠在尋源過程中實時進行以上參數(shù)的優(yōu)化。

1.3.1 徑向基函數(shù)分配

通過資源分配網(wǎng)絡算法(resource-allocating network, 簡稱RAN)[13]進行徑向基函數(shù)中心的分配。RAN算法開始時具有較少徑向基函數(shù)中心，根據(jù)新增樣本的新穎性決定是否添加徑向基函數(shù)進行誤差消減。

條件1：當前樣本位置xp+1與最近的徑向基函數(shù)中心Cnear的歐式距離超過閾值δ。

‖xp+1-Cnear‖>δ

(7)

條件2：神經(jīng)網(wǎng)絡對當前輸入樣本xp+1的預測與實測數(shù)值偏差大于閾值ε。

(8)

添加徑向基函數(shù)后神經(jīng)網(wǎng)絡發(fā)生變化，需要通過更新正則化參數(shù)λ和計算權重向量W進行神經(jīng)網(wǎng)絡的優(yōu)化。

1.3.2 在線正則化參數(shù)迭代計算

神經(jīng)網(wǎng)絡的性能取決于預測的準確度和泛化性，正則化參數(shù)λ的引入能夠在一定程度上實現(xiàn)二者的平衡，而λ的取值可以通過廣義交叉驗證誤差(generalized cross-validation, 簡稱GCV)的最小化來決定。首先給出徑向基函數(shù)神經(jīng)網(wǎng)絡的廣義交叉驗證誤差表達式：

(9)

其中，p為當前路徑點個數(shù)，P為RBFNN中的投影矩陣：

P=I-HA-1HT

(10)

(11)

值得注意的是若直接使用公式(11)會導致每次迭代計算都需進行多次復雜度為O(m3)的矩陣逆運算。觀察與矩陣A相關的表達式(5)可以發(fā)現(xiàn)，通過設計矩陣H的奇異值分解能夠有效簡化A的求逆過程進而簡化公式(11)各部分的計算復雜度。

首先設截斷奇異值分解(truncated SVD)H=USVT，其中U=[u1u2…ur]∈Rp×r、V∈Rm×r為正交矩陣，S∈Rr×r且有如下形式：

(12)

式(11)各部分簡化結果如下[14]：

(13)

通過以上簡化計算，能夠?qū)⑹?11)的復雜度從O(m3)降低到O(p2),并且只需要在迭代計算前進行一次復雜度為O(p3)的奇異值分解計算，對于重復迭代的計算而言能夠節(jié)省大量算力，算法的偽代碼如下：

Algorithm 1 基于SVD的正則化系數(shù)求解Require:SVD分解USVT=H,函數(shù)值Y^Ensure:正則化系數(shù)λ,權重向量:W1:function REGULARIZE(U,S,V,Y^)2:λ0,λp1,μdiag(S),yUTY^3:while(abs(λ-λp)>0.000 1)do4. λpλ5:λ∑pi=1uiy2i(ui+λp)3∑pi=1λpui+λp∑pi=1λ2py2i(ui+λp)2∑pi=1ui(ui+λp)26:end while7:WV(STS+λI)-1Sy8:return[λ W]9:end function

對不同規(guī)模的設計矩陣H利用原版迭代算法和基于SVD改進的迭代算法進行正則化參數(shù)求解，設定矩陣H為方陣。設置初始λ=1，并進行50次的迭代。統(tǒng)計改進算法中，總耗時、SVD分解耗時以及迭代耗時與原算法耗時的比值，結果如圖2所示。

圖2 迭代公式計算耗時對比Fig. 2 Time cost of different re-estimation method

從結果中可以觀察到，改進算法中迭代計算的耗時與原算法的耗時比值隨著神經(jīng)網(wǎng)絡規(guī)模的增大而減小，改進算法中的SVD分解耗時與原算法的耗時比值卻維持在1/50附近。整體而言改進算法的耗時遠小于原算法，但是其中SVD分解耗時的復雜度階數(shù)與原算法相似，需要進一步簡化。

1.3.3 設計矩陣H增量式SVD分解

對設計矩陣H進行奇異值分解能夠顯著降低正則化系數(shù)求解的復雜度，卻引入了SVD分解的計算復雜度隨神經(jīng)網(wǎng)絡階數(shù)增加而增加的問題。因此奇異值分解的簡化變得至關重要。觀察式(2)中設計矩陣H的結構能夠發(fā)現(xiàn)，增加樣本數(shù)據(jù)和增加徑向基函數(shù)實質(zhì)上是在增加設計矩陣H的行數(shù)或者列數(shù)。這意味著每次神經(jīng)網(wǎng)絡結構變化時H的主體部分不發(fā)生變動。結合Brand的工作[15]利用Hp,m的奇異值分解結果對Hp,m+1或Hp+1,m進行增量式SVD分解，定義Hp,m+1和Hp+1,m分別為增加徑向基函數(shù)和增加樣本數(shù)據(jù)后的設計矩陣。以下以Hp,m+1為例簡述增量式SVD的計算過程。

對設計矩陣Hm有截斷奇異值分解，其中Sm階數(shù)為r：

Hp,m=USVT

(14)

Hp,m+1與Hp,m之間存在如下關系：

Hp,m+1=[Hp,mhm+1]，hm+1=[hm+1(x1) ……h(huán)m+1(xp)]T

(15)

計算：

(16)

對M進行歸一化處理：

J=M/K,K=‖M‖

(17)

設計矩陣運算如下：

(18)

Q=U′S′V′T

(19)

可以得到更新后的奇異值分解：

Hp,m+1=U″S″V″T

(20)

其中，

(21)

通過以上的分解能夠?qū)p,m的奇異值分解轉(zhuǎn)化為Qr+1,r+1的奇異值分解。并且值得注意的是Q實際上為邊對角矩陣，其奇異值分解的計算復雜度可以通過O″ Leary[16]的工作從O(r3)簡化到O(r2)。通過以上的簡化計算能夠?qū)⒃O計矩陣奇異值分解的計算復雜度降低到O((m+p)r2)，在計算中設定比值r/p=0.25。

Algorithm 2 增量式SVD分解Require:原分解:USVT=H新增列:hm+1Ensure:增量SVD: U″,S″,V″,h″T=[H h]1:function INC SVD(U,S,V,h)2:[U'S'V']SVDSUThS|h-UUTh| 3:U″[U(h-UUTL)/|h-UUTh|]U'4:S″S'5:V″V001 V'6:return[U″,S″,V″]7:end function

對不同階數(shù)的設計矩陣增加列數(shù)據(jù)并求取新設計矩陣的奇異值分解。計算結果如圖3所示，能夠直觀地看到增量式分解顯著降低了計算耗時并且效果隨著神經(jīng)網(wǎng)絡規(guī)模的增大而增大，計算過程在Matlab2019 中進行，考慮到數(shù)學計算軟件對于計算過程的優(yōu)化，結果中呈現(xiàn)的計算耗時比例未必能夠完全與理論上的計算復雜度分析完全對應。

圖3 SVD分解耗時對比Fig. 3 Time cost of different SVD method

2 自主尋源算法

第1節(jié)的工作能夠支持航行器進行在線的信號場擬合，如何規(guī)劃航行器的航行路徑則是本節(jié)重點。算法的主體思想實際上是代理優(yōu)化，通過擬合場的梯度信息進行未知場的探索，首先應該進行擬合場的梯度求解。

2.1 擬合函數(shù)梯度求解

由函數(shù)之和的求導法則可知，路徑點梯度為所有組成函數(shù)在該路徑點梯度之和：

(22)

已知徑向基函數(shù)形式，可得徑向基函數(shù)hj在當前路徑點梯度：

(23)

對式(23)進行求和可以得擬合函數(shù)在路徑點xi的梯度gi為：

(24)

2.2 脫離局部最優(yōu)區(qū)域

在AUV自主尋優(yōu)的研究中，航行器陷入局部最優(yōu)區(qū)域是難以避免的問題。類似于遺傳算法、模擬退火、粒子群算法等群體智能優(yōu)化算法可能適用于多航行器的尋優(yōu)問題，而對于單航行器的尋優(yōu)問題而言，考慮到航行路徑、搜索時間等約束，應用以上群體優(yōu)化算法是不切實際的。

采用RBFNN對未知分布場進行“擬合估計—最大梯度搜索—重新擬合”的迭代擬合搜索方法，實際上是一種局部鄰域搜索算法，相應的存在著拓展算法如禁忌搜索算法和動量梯度算法。禁忌搜索[17]算法的關鍵思想是通過禁忌表標記已經(jīng)搜索過的區(qū)域并在后續(xù)步驟中避免訪問，通過限制訪問區(qū)域來迫使算法接受次優(yōu)解，并有一定幾率訪問到局部最優(yōu)解之外的區(qū)域并借此脫離局部最大值區(qū)域。該算法的以上特性會導致AUV路徑在局部最優(yōu)區(qū)域附近滯留，而且脫離局部區(qū)域的路徑具有隨機性，不適用于文中研究。

動量梯度法作為梯度法的拓展[18]，其算法核心是當前路徑點xi并非直接采用梯度方向作為前進方向，而是選擇計算所有歷史路徑點梯度的指數(shù)加權平均數(shù)作為前進方向。記當前路徑點梯度為gi=?fpre(xi,D)，則求取下一路徑點步驟可由式(25)～(26)表示。

(25)

xi+1=xi+hDiri

(26)

其中，μ=0.9為衰減系數(shù)，h=10為航行器航行步長，Diri為路徑點對應前進方向。

動量梯度算法借鑒了物理學的概念，使得AUV在搜索過程中通過計算梯度指數(shù)加權平均數(shù)積攢“能量”。其中衰減系數(shù)μ起著至關重要的作用，當衰減系數(shù)μ=0時，動量梯度法和正常的梯度法相同，當μ=1時，情況類似于無摩擦運動，路徑難以在極大值附近收斂。相比于梯度法，動量梯度法具有易于逃脫局部最優(yōu)區(qū)域和更好通過信號場高原區(qū)的優(yōu)點。

2.3 算法流程

文中算法主要內(nèi)容形成流程如圖4所示，首先進行神經(jīng)網(wǎng)絡的初始化：選擇出發(fā)點坐標，在出發(fā)點附近生成一個徑向基函數(shù)；然后進行流程的循環(huán)直至航行器的二維坐標收斂：網(wǎng)絡正則化迭代—權重向量計算—當前路徑點梯度計算—使用動量梯度法生成航行器下一步的前進方向—獲取下一路徑點信號值—判定新樣本的新穎性并決定是否添加徑向基函數(shù)。

圖4 尋源算法流程Fig. 4 Flow diagram of source-seeking algorithm

3 場源搜索過程模擬計算

3.1 單峰場場源搜索模擬

為驗證尋源算法的有效性，模擬了一個熱液噴口附近硫化氫濃度的單極值信號場S，如圖5(a)所示，該分布場坐標范圍為1 000 m×1 000 m，最大值為Smax=25.06 mmol/kg，位于坐標點Pglomax(629,327)。海洋環(huán)境中許多信號值隨距離增加有極大的衰減[19]，在遠離場源的區(qū)域信號場梯度較小，文中設計的硫化氫場便具有如此的性質(zhì)，對該信號場求取梯度，得場梯度分布如圖5(b)所示。

圖5 單峰硫化氫分布場Fig. 5 Single-peak field of H2S distribution

設置出發(fā)點分別為P1(100,100)、P2(100,500)、P3(100,900)、P4(500,900)、P5(900,900)、P6(900,500)、P7(900,100)、P8(500,100)進行尋源過程的模擬計算，為了驗證算法穩(wěn)定性，在真實信號中添加均值為0，方差為0.1的高斯噪聲。每次出發(fā)點進行1 000次仿真并計算平均路徑長度、最大值誤差、最大值坐標誤差、神經(jīng)網(wǎng)絡的GCV誤差及全場均方誤差MSE等，計算結果如表1所示。注意到，不同出發(fā)點的最大值以及最大值坐標誤差基本相同，說明該算法的收斂結果與出發(fā)點無關，收斂的精度實際上與航行器的步長相關。徑向基函數(shù)神經(jīng)網(wǎng)絡無法對遠離樣本的區(qū)域進行準確預測，因此盡管GCV誤差僅為0.01的水平，但是全場的均方誤差卻相對要高。由于神經(jīng)網(wǎng)絡能夠?qū)颖緟^(qū)域進行精準地預測，該算法對于信號場路徑點上的梯度估計有著誤差低于20°的精度，能夠為航行器前進方向的計算提供良好基礎。

表1 單峰場尋源算法仿真結果Tab. 1 Results of source-seeking simulation in single-peak field

文中將文獻[10，20]中的尋源算法進行相同出發(fā)點的尋源過程模擬并對比其結果。文獻[20]通過控制航行器執(zhí)行圓形路徑機動進行梯度的最小二乘估計并尋找場源。文獻[10]方法基于考慮先驗場的高斯過程回歸(GPR)，進行場源搜索仿真。將3種尋源算法的路徑點分布(出發(fā)點3)、尋源過程路徑長度、路徑點梯度估計誤差以及GPR方法和RBF方法的全場均方誤差MSE進行對比，如圖6～9所示。

圖6 單峰值信號場不同算法尋源過程路徑點分布(P3)Fig. 6 Way point distribution of different source-seeking algorithms in single-peak field (P3)

圖7 均方誤差MSEFig. 7 Mean square error

圖8 路徑點梯度方向預測誤差Fig. 8 Average gradient direction error of track

圖9 尋源過程平均路徑Fig. 9 Average length of seeking procedure

分析以上數(shù)據(jù)可以發(fā)現(xiàn)，圓形機動估計梯度算法的路徑點梯度估計誤差最大，GPR算法由于存在先驗信號場的影響，梯度估計不如使用無先驗信息的RBFNN算法準確，圖6中顯示基于RBFNN算法產(chǎn)生的路徑集中于連接出發(fā)點和極值點的直線上，并且由于采用了動量梯度法，其路徑往往會在經(jīng)過極值點后向前延伸一段；GPR算法生成的路徑大部分會經(jīng)過先驗信號場的極值區(qū)域，這導致其路徑相對較長并且均方誤差大于基于RBFNN的算法。而基于圓形估計梯度的算法生成的路徑散亂，在梯度較小的區(qū)域由于信號噪聲的存在不能夠準確地估計梯度，又由于圓形路徑估計梯度的過程存在，其路徑遠大于其他兩個算法。

3.2 多峰場場源搜索模擬

改變原H2S濃度分布場，在原有的H2S濃度分布場S上添加信號場S1：

(27)

該分布場(圖10)中局部最大值區(qū)域位于Qlocmax(450,750)附近，全局極大值位于Qglomax(629,329)。考慮到本小節(jié)旨在測試尋源算法脫離局部最優(yōu)區(qū)域的能力，設計可能陷入局部最大值區(qū)域的出發(fā)點Q1(200,900)、Q2(300,900)、Q3(400,900)以及Q4(500,900)，對其尋源路徑進行多次模擬計算，直至每個出發(fā)點的尋源成功率收斂。最后得出不同出發(fā)點的平均路徑長度、脫離局部最優(yōu)區(qū)域的成功率等數(shù)據(jù)如表2所示，并繪制出發(fā)點1和出發(fā)點4的尋源過程路徑分布如圖11所示。結果顯示算法能夠以較高的成功率通過局部最大值區(qū)域，但是由于局部最大值區(qū)域梯度的干擾，AUV的路徑可能在該區(qū)域彎曲，導致整個尋源過程總路徑長度的增加。

圖10 多峰場硫化氫分布場Fig. 10 Multi-peaks field of H2S distribution

圖11 多峰值信號場不同出發(fā)點尋源路徑點分布Fig. 11 Distribution of source seeking way points from different starting points in multi-peak field

表2 多峰場尋源算法仿真結果Tab. 2 Results of source-seeking simulation in multi-peak field

4 結 語

針對水下航行器在未知環(huán)境中的場源搜索問題進行了研究，利用徑向基神經(jīng)網(wǎng)絡對未知信號場進行擬合，通過引入正則化參數(shù)提高神經(jīng)網(wǎng)絡的泛化性，以廣義交叉驗證誤差為準則優(yōu)化正則化參數(shù)，并以增量式奇異值分解降低正則化參數(shù)求解過程的計算復雜度，最后通過擬合場的梯度信息以及動量梯度法規(guī)劃水下航行器的運動方向。

與其他研究中的尋源算法進行尋源模擬對比發(fā)現(xiàn)：文中方法能夠?qū)β窂近c的梯度進行更為準確的估計，從而生成從出發(fā)點到場源更短的搜索路徑，并且在全場范圍內(nèi)的均方誤差更小。此外進行了多峰值信號場的搜索模擬計算，結果顯示文中的自主尋源算法能夠以較高的成功率脫離局部最優(yōu)區(qū)域。在這里僅考慮單航行器對未知信號場的搜索問題，后續(xù)的工作將基于現(xiàn)有的神經(jīng)網(wǎng)絡在線預測算法針對多航行器的協(xié)同搜索問題開展。