吳宗秀,吳 超
(1. 上海交通大學 船舶海洋與建筑工程學院,上海 200240; 2. 上海交通大學 海洋工程國家重點實驗室,上海 200240)
水下尋源著眼于解決利用單個或編隊水下航行器對未知分布的信號場進行信號源搜索問題。水下工程領域的許多問題都可以歸結為水下尋源問題。包括:失事飛行器黑匣子的尋找、天然氣或者水下礦藏的定位、熱液噴口搜索[1]等。在這些問題中,未知信號場則相應地可以具體化為聲學信號場、化學物質濃度信號場及溫度信號場等。水下尋源問題的解決方案對于擴展水下航行器的應用而言有著重大的意義,該領域的解決方法一般有瞬態決策方法、極值搜索算法、高斯過程回歸以及神經網絡等。
Holland等[2]、Russell[3]分別就未知信號場的尋源問題提出了六邊形和Z字型前進路線的搜索策略,通過離散的信號測量以及簡單的信號強度對比進行搜索路徑的規劃。Russell[4]后續將工作拓展到三維場源搜索問題中,瞬態極值策略可以作為尋源問題一種相對穩定的解決方案,但是其決策得出的方向精度不足。自動控制領域的方法在場源搜索問題中亦得到了廣泛的應用,Mellucci等[5]、Matveev等[6-7]均運用了滑??刂品椒▉砜刂扑潞叫衅鞯氖紫蜻M行信號場的極值搜索,其結果顯示,滑??刂品椒ㄉ傻乃唁浡窂较鄬^長。Cochran等[8]將極值搜索算法(extreme value search algorithm, 簡稱ESA)運用到水下尋源問題中,通過控制航行器的轉向引導航行器進行搜索,該團隊在其他研究[9]中增加前進速度為控制參數并且得到了更好的結果,但同時也指出該算法收斂速度較慢的問題。
高斯過程回歸(gaussian process regression, 簡稱GPR)具有良好的函數逼近能力,常用于解決未知信號場的擬合問題。Ai等[10]結合高斯過程回歸、Armijo步長準則、遞歸最小二乘方法和非線性模型預測控制(nonlinear model predictive control, 簡稱NMPC)進行完備的水下航行器場源搜索仿真。其研究結果表明航行器在搜索過程前期受先驗信號場影響明顯,在接近先驗信號場的最大值區域后轉向正確的路徑并收斂到真實場源。閻述學等[11]的研究著眼于利用高斯過程回歸進行水下信號場的熱點采樣,通過定義一種與均方差相關的引力模型引導航行器進行多熱點未知信號場的自適應采樣。相較于傳統的割草機和梳型路徑,該算法能夠以較短的路徑實現更準確的擬合。
徑向基神經網絡(radial basis function neural network,簡稱RBFNN)具備無限精度的函數逼近能力,能夠運用于信號場的估計并且具有結構簡單、易于理解的優點。Jeon等[12]通過推導連續時間形式的遞歸最小二乘法對RBFNN進行在線訓練,并控制裝有多個傳感器的雙連桿機械手進行二維灰度圖的極值尋找。然而利用RBFNN的尋源算法研究很少考慮神經網絡的正則化,實際上正則化對于神經網絡的穩定性和泛化性非常重要。
以水下自主航行器(autonomous underwater vehicle, 簡稱AUV)在熱泉區域的硫化氫濃度分布場場源搜索為背景,使用具有增長式結構的RBFNN對未知信號場進行在線無先驗漸進擬合,并通過增量式奇異值分解算法加速神經網絡的正則化計算,此外還利用動量梯度法解決AUV路徑容易陷入局部最優區域的問題。算法的有效性通過單峰值和多峰值硫化氫濃度場中的場源搜索模擬計算進行驗證。

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

圖1 RBFNN結構示意Fig. 1 View of RBFNN’s construction
RBFNN在結構上是嚴格的三層網絡,分別為輸入層、隱藏層和輸出層。輸入層負責信號傳遞,隱藏層通過徑向基函數(radial basis function, 簡稱RBF)進行從輸入層到隱藏層空間的非線性變換,輸出層對隱藏層的響應進行加權求和。
選用高斯型徑向基函數:

(1)
其中,ωj為權重系數,cj為核函數中心,rj為高斯函數半徑。
設計矩陣H,其元素Hi,j=hj(xi),行數記作p,列數記作m:
(2)
權重向量W可表示為:
(3)


(4)
權重向量可以通過計算偽逆或最小二乘法求取。為了避免矩陣求解中的病態問題以及神經網絡的過擬合問題,引入全局正則化操作,對應權重向量求解公式如下:
(5)
其中,參數λ為全局正則化系數,該參數能夠提高神經網絡的泛化性和計算穩定性,λ的數值需要在網絡結構和數據集變化的時候進行優化。
綜上,RBFNN對任意輸入x∈R2的信號值預測為:
(6)
徑向基神經網絡的預測性能取決于RBF中心{C1C2……Cm}、權重向量W、全局正則化參數λ等。為了實現RBFNN的在線訓練,算法應當能夠在尋源過程中實時進行以上參數的優化。
1.3.1 徑向基函數分配
通過資源分配網絡算法(resource-allocating network, 簡稱RAN)[13]進行徑向基函數中心的分配。RAN算法開始時具有較少徑向基函數中心,根據新增樣本的新穎性決定是否添加徑向基函數進行誤差消減。

條件1:當前樣本位置xp+1與最近的徑向基函數中心Cnear的歐式距離超過閾值δ。
‖xp+1-Cnear‖>δ
(7)
條件2:神經網絡對當前輸入樣本xp+1的預測與實測數值偏差大于閾值ε。
(8)
添加徑向基函數后神經網絡發生變化,需要通過更新正則化參數λ和計算權重向量W進行神經網絡的優化。
1.3.2 在線正則化參數迭代計算
神經網絡的性能取決于預測的準確度和泛化性,正則化參數λ的引入能夠在一定程度上實現二者的平衡,而λ的取值可以通過廣義交叉驗證誤差(generalized cross-validation, 簡稱GCV)的最小化來決定。首先給出徑向基函數神經網絡的廣義交叉驗證誤差表達式:
(9)
其中,p為當前路徑點個數,P為RBFNN中的投影矩陣:
P=I-HA-1HT
(10)

(11)
值得注意的是若直接使用公式(11)會導致每次迭代計算都需進行多次復雜度為O(m3)的矩陣逆運算。觀察與矩陣A相關的表達式(5)可以發現,通過設計矩陣H的奇異值分解能夠有效簡化A的求逆過程進而簡化公式(11)各部分的計算復雜度。
首先設截斷奇異值分解(truncated SVD)H=USVT,其中U=[u1u2…ur]∈Rp×r、V∈Rm×r為正交矩陣,S∈Rr×r且有如下形式:
(12)
式(11)各部分簡化結果如下[14]:
(13)
通過以上簡化計算,能夠將式(11)的復雜度從O(m3)降低到O(p2),并且只需要在迭代計算前進行一次復雜度為O(p3)的奇異值分解計算,對于重復迭代的計算而言能夠節省大量算力,算法的偽代碼如下:

Algorithm 1 基于SVD的正則化系數求解Require:SVD分解USVT=H,函數值Y^Ensure:正則化系數λ,權重向量: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
對不同規模的設計矩陣H利用原版迭代算法和基于SVD改進的迭代算法進行正則化參數求解,設定矩陣H為方陣。設置初始λ=1,并進行50次的迭代。統計改進算法中,總耗時、SVD分解耗時以及迭代耗時與原算法耗時的比值,結果如圖2所示。

圖2 迭代公式計算耗時對比Fig. 2 Time cost of different re-estimation method
從結果中可以觀察到,改進算法中迭代計算的耗時與原算法的耗時比值隨著神經網絡規模的增大而減小,改進算法中的SVD分解耗時與原算法的耗時比值卻維持在1/50附近。整體而言改進算法的耗時遠小于原算法,但是其中SVD分解耗時的復雜度階數與原算法相似,需要進一步簡化。
1.3.3 設計矩陣H增量式SVD分解
對設計矩陣H進行奇異值分解能夠顯著降低正則化系數求解的復雜度,卻引入了SVD分解的計算復雜度隨神經網絡階數增加而增加的問題。因此奇異值分解的簡化變得至關重要。觀察式(2)中設計矩陣H的結構能夠發現,增加樣本數據和增加徑向基函數實質上是在增加設計矩陣H的行數或者列數。這意味著每次神經網絡結構變化時H的主體部分不發生變動。結合Brand的工作[15]利用Hp,m的奇異值分解結果對Hp,m+1或Hp+1,m進行增量式SVD分解,定義Hp,m+1和Hp+1,m分別為增加徑向基函數和增加樣本數據后的設計矩陣。以下以Hp,m+1為例簡述增量式SVD的計算過程。
對設計矩陣Hm有截斷奇異值分解,其中Sm階數為r:
Hp,m=USVT
(14)
Hp,m+1與Hp,m之間存在如下關系:
Hp,m+1=[Hp,mhm+1],hm+1=[hm+1(x1) ……hm+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)
通過以上的分解能夠將Hp,m的奇異值分解轉化為Qr+1,r+1的奇異值分解。并且值得注意的是Q實際上為邊對角矩陣,其奇異值分解的計算復雜度可以通過O″ Leary[16]的工作從O(r3)簡化到O(r2)。通過以上的簡化計算能夠將設計矩陣奇異值分解的計算復雜度降低到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
對不同階數的設計矩陣增加列數據并求取新設計矩陣的奇異值分解。計算結果如圖3所示,能夠直觀地看到增量式分解顯著降低了計算耗時并且效果隨著神經網絡規模的增大而增大,計算過程在Matlab2019 中進行,考慮到數學計算軟件對于計算過程的優化,結果中呈現的計算耗時比例未必能夠完全與理論上的計算復雜度分析完全對應。

圖3 SVD分解耗時對比Fig. 3 Time cost of different SVD method
第1節的工作能夠支持航行器進行在線的信號場擬合,如何規劃航行器的航行路徑則是本節重點。算法的主體思想實際上是代理優化,通過擬合場的梯度信息進行未知場的探索,首先應該進行擬合場的梯度求解。
由函數之和的求導法則可知,路徑點梯度為所有組成函數在該路徑點梯度之和:

(22)
已知徑向基函數形式,可得徑向基函數hj在當前路徑點梯度:

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

(24)
在AUV自主尋優的研究中,航行器陷入局部最優區域是難以避免的問題。類似于遺傳算法、模擬退火、粒子群算法等群體智能優化算法可能適用于多航行器的尋優問題,而對于單航行器的尋優問題而言,考慮到航行路徑、搜索時間等約束,應用以上群體優化算法是不切實際的。
采用RBFNN對未知分布場進行“擬合估計—最大梯度搜索—重新擬合”的迭代擬合搜索方法,實際上是一種局部鄰域搜索算法,相應的存在著拓展算法如禁忌搜索算法和動量梯度算法。禁忌搜索[17]算法的關鍵思想是通過禁忌表標記已經搜索過的區域并在后續步驟中避免訪問,通過限制訪問區域來迫使算法接受次優解,并有一定幾率訪問到局部最優解之外的區域并借此脫離局部最大值區域。該算法的以上特性會導致AUV路徑在局部最優區域附近滯留,而且脫離局部區域的路徑具有隨機性,不適用于文中研究。
動量梯度法作為梯度法的拓展[18],其算法核心是當前路徑點xi并非直接采用梯度方向作為前進方向,而是選擇計算所有歷史路徑點梯度的指數加權平均數作為前進方向。記當前路徑點梯度為gi=?fpre(xi,D),則求取下一路徑點步驟可由式(25)~(26)表示。

(25)
xi+1=xi+hDiri
(26)
其中,μ=0.9為衰減系數,h=10為航行器航行步長,Diri為路徑點對應前進方向。
動量梯度算法借鑒了物理學的概念,使得AUV在搜索過程中通過計算梯度指數加權平均數積攢“能量”。其中衰減系數μ起著至關重要的作用,當衰減系數μ=0時,動量梯度法和正常的梯度法相同,當μ=1時,情況類似于無摩擦運動,路徑難以在極大值附近收斂。相比于梯度法,動量梯度法具有易于逃脫局部最優區域和更好通過信號場高原區的優點。
文中算法主要內容形成流程如圖4所示,首先進行神經網絡的初始化:選擇出發點坐標,在出發點附近生成一個徑向基函數;然后進行流程的循環直至航行器的二維坐標收斂:網絡正則化迭代—權重向量計算—當前路徑點梯度計算—使用動量梯度法生成航行器下一步的前進方向—獲取下一路徑點信號值—判定新樣本的新穎性并決定是否添加徑向基函數。

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

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

表1 單峰場尋源算法仿真結果Tab. 1 Results of source-seeking simulation in single-peak field
文中將文獻[10,20]中的尋源算法進行相同出發點的尋源過程模擬并對比其結果。文獻[20]通過控制航行器執行圓形路徑機動進行梯度的最小二乘估計并尋找場源。文獻[10]方法基于考慮先驗場的高斯過程回歸(GPR),進行場源搜索仿真。將3種尋源算法的路徑點分布(出發點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
分析以上數據可以發現,圓形機動估計梯度算法的路徑點梯度估計誤差最大,GPR算法由于存在先驗信號場的影響,梯度估計不如使用無先驗信息的RBFNN算法準確,圖6中顯示基于RBFNN算法產生的路徑集中于連接出發點和極值點的直線上,并且由于采用了動量梯度法,其路徑往往會在經過極值點后向前延伸一段;GPR算法生成的路徑大部分會經過先驗信號場的極值區域,這導致其路徑相對較長并且均方誤差大于基于RBFNN的算法。而基于圓形估計梯度的算法生成的路徑散亂,在梯度較小的區域由于信號噪聲的存在不能夠準確地估計梯度,又由于圓形路徑估計梯度的過程存在,其路徑遠大于其他兩個算法。
改變原H2S濃度分布場,在原有的H2S濃度分布場S上添加信號場S1:
(27)
該分布場(圖10)中局部最大值區域位于Qlocmax(450,750)附近,全局極大值位于Qglomax(629,329)。考慮到本小節旨在測試尋源算法脫離局部最優區域的能力,設計可能陷入局部最大值區域的出發點Q1(200,900)、Q2(300,900)、Q3(400,900)以及Q4(500,900),對其尋源路徑進行多次模擬計算,直至每個出發點的尋源成功率收斂。最后得出不同出發點的平均路徑長度、脫離局部最優區域的成功率等數據如表2所示,并繪制出發點1和出發點4的尋源過程路徑分布如圖11所示。結果顯示算法能夠以較高的成功率通過局部最大值區域,但是由于局部最大值區域梯度的干擾,AUV的路徑可能在該區域彎曲,導致整個尋源過程總路徑長度的增加。

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

圖11 多峰值信號場不同出發點尋源路徑點分布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
針對水下航行器在未知環境中的場源搜索問題進行了研究,利用徑向基神經網絡對未知信號場進行擬合,通過引入正則化參數提高神經網絡的泛化性,以廣義交叉驗證誤差為準則優化正則化參數,并以增量式奇異值分解降低正則化參數求解過程的計算復雜度,最后通過擬合場的梯度信息以及動量梯度法規劃水下航行器的運動方向。
與其他研究中的尋源算法進行尋源模擬對比發現:文中方法能夠對路徑點的梯度進行更為準確的估計,從而生成從出發點到場源更短的搜索路徑,并且在全場范圍內的均方誤差更小。此外進行了多峰值信號場的搜索模擬計算,結果顯示文中的自主尋源算法能夠以較高的成功率脫離局部最優區域。在這里僅考慮單航行器對未知信號場的搜索問題,后續的工作將基于現有的神經網絡在線預測算法針對多航行器的協同搜索問題開展。