基于改進螢火蟲算法的船舶裝卸作業(yè)穩(wěn)性研究

肖力旗，曹 民，張柳明

（1.上海理工大學光電信息與計算機工程學院，上海 200082；2.上海交通大學電子信息與電氣工程學院，上海 200240）

0 引言

在經(jīng)濟全球化時代背景下，海洋運輸在國際貿易、大型貨物運輸中擁有無可替代的地位，船舶是海洋運輸最主要的交通工具，船舶穩(wěn)性［1-2］是保障航運安全的重要指標。國際海事組織（International Maritime Organization，IMO）在第二代完整穩(wěn)性衡準中確定了5 種失效模型［1］，以此規(guī)范不同狀態(tài)船舶的穩(wěn)性問題。在碼頭裝卸作業(yè)工況下，船舶穩(wěn)性易受諸多不確定因素影響［3-4］，研究該工況下的船舶穩(wěn)性對碼頭作業(yè)安全具有重要意義。

針對船舶穩(wěn)性的研究主要集中在船舶耐波性，研究方向主要包括兩個部分：①在時域上研究特定波浪對船舶運動狀態(tài)影響［5］；②基于統(tǒng)計思想，利用譜分析理論分析船舶在不規(guī)則波中的頻域響應。Zhang 等［6-7］研究了規(guī)則和非規(guī)則波中參數(shù)橫搖的數(shù)學模型，分析了穩(wěn)性高度對船舶穩(wěn)性的影響；朱軍等［8］用搖晃耦合切片法分析橫搖運動規(guī)律；Chakrabarti［9］對非線性橫搖阻尼影響進行評估，結合實驗提出在改進的經(jīng)驗公式；Yin 等［10］提出基于自適應滑動數(shù)據(jù)窗的船舶橫搖預報方案，實現(xiàn)船舶橫搖動態(tài)在線預測，這些研究為船舶橫搖和穩(wěn)性分析提供了研究基礎。然而在港口風浪較小情況下，系泊裝卸作業(yè)工況下船舶失穩(wěn)主要受重物激勵和重心變化影響，利用譜分析理論分析船舶在該工況下的穩(wěn)性并不適用。

近年來，智能控制算法發(fā)展迅猛，算法不斷優(yōu)化，逐漸深入到各工業(yè)領域。Wu 等［11］利用模糊算法研究了船舶失穩(wěn)各因素影響情況，但研究對象仍為航行中的船舶，利用智能控制算法研究裝卸作業(yè)的船舶穩(wěn)性較少。

針對系泊船舶裝卸作業(yè)特殊的工況及穩(wěn)性監(jiān)測實時性問題，本研究通過分析橫搖周期與初穩(wěn)性高的關系，利用改進的螢火蟲算法優(yōu)化船舶橫搖運動方程模型，預測船舶橫搖運動狀態(tài)，計算出船舶初穩(wěn)性高，從而分析系泊船舶穩(wěn)性。采用仿真實例驗證了該方法的可行性和有效性。

1 模型理論依據(jù)

1.1 小傾角穩(wěn)性

小傾角穩(wěn)性也稱初穩(wěn)性，一般情況下小傾角指傾斜角小于10°～15°［12］。在正常氣象狀況下，港區(qū)船舶屬于小傾角穩(wěn)性問題，通常用初穩(wěn)性高來描述該狀態(tài)下的船舶穩(wěn)性。2008 年IMO 船舶穩(wěn)性規(guī)則中規(guī)定船舶初穩(wěn)性高不應低于0.15m［13］。船舶在靜水中的初穩(wěn)性高與橫搖周期存在如下關系［14］：

式中，GM為初穩(wěn)性高，T為橫搖周期，B為型寬，d為平均吃水，Lwl為設計水線長。

1.2 船舶橫搖運動方程

對于船舶的橫搖運動，若不考慮風浪流等外界干擾，船舶的橫搖運動力矩方程［15］為：

其中，M()是慣性力矩，與角加速度相關，M()是阻尼力矩，與角速度相關，M()是回復力矩，由傾斜角度決定。對于線性阻尼，船舶橫搖運動方程可寫為［16-17］：

其 中，θm是橫搖最大傾角，vθ是橫搖衰減系數(shù)，vθ=μθnθ，μθ為無因次衰減系數(shù)，nθ為船舶橫搖固有頻率，εθ是其相位角。

船舶在靜水中受到外界一定傾角的作用下，其橫搖做衰減振蕩。由于在船舶裝載過程中很難保持船舶左右兩側的平衡，因此在船舶橫搖幅度衰減到極小時可能會出現(xiàn)穩(wěn)態(tài)傾角θ0。若考慮穩(wěn)態(tài)傾角情況，則橫搖公式可寫為：

由于船舶裝載作業(yè)是間斷性的加裝貨物，本研究需要在盡可能短的時間內根據(jù)已檢測到的船舶傾角信息，將其橫搖周期及穩(wěn)態(tài)傾角解算出來。式（5）中，已知一定長度的傾角和時間序列，而還有θm、μθ、nθ、θ0、εθ五個參數(shù)是未知的，對于5 個未知參數(shù)的非線性曲線擬合其計算量相當大，無法滿足在一次裝載間隔內計算出橫搖周期需求，因此需要根據(jù)其他特征確定部分參數(shù)。

1.3 傾角序列前期處理

實際工作環(huán)境中，船舶在港口裝載貨物時，其工作環(huán)境風浪較小，但仍然不可忽略其影響。由于風浪噪聲屬于混沌系統(tǒng)，其干擾函數(shù)難以確定，通常使用多項式擬合方法使傾角信號更加平滑。

假設采集到的傾角序列為{θ1,θ2,θ3...θj}，按一定長度截取序列進行分段多項式擬合，利用來替代原信號的值。tj為采集到θj傾角值對應的時間，傾角信號采樣頻率T=tj+1-tj，ai為多項式系數(shù)。將ai以某一初值進行迭代，當達到最小時獲得最佳多項式系數(shù)aj的值。采用梯度下降法來加快迭代速度。在多項式擬合過程中，為提高擬合速度、防止過擬合，選取的多項式階數(shù)不應過高，截取的序列長度應與多項式階數(shù)相關。

對于已經(jīng)經(jīng)過多項式擬合的數(shù)據(jù)，其波峰（波谷）處一階導數(shù)為0。在一個激勵作用下，橫搖曲線是隨時間衰減的曲線，橫搖最大傾角即為第一次出現(xiàn)的波峰（波谷）傾角θmax。將出現(xiàn)橫搖最大傾角的時間標定為0 時刻，有εθ=0（若第一個為波谷，則εθ=π/2）。將θmax、t=0 代入式中，可得θmax=θm+θ0。則式（5）可寫為：

在該式中，存在θm、μθ、nθ三個未知參數(shù)，要得到橫搖曲線參數(shù)的最優(yōu)解，本研究使用螢火蟲算法進行橫搖參數(shù)求解。

2 改進螢火蟲算法

2.1 標準螢火蟲算法

螢火蟲算法（Firefly Algorithm，F(xiàn)A）是一種基于自然啟發(fā)的群智能優(yōu)化算法，來源于螢火蟲發(fā)光并相互吸引的行為［18］，其原則可簡化為：

（1）每個螢火蟲相互吸引無性別區(qū)分。

（2）相對亮度與距離平方成反比，當距離過大時無吸引力。

（3）亮度取決于是否接近目標函數(shù)的最佳優(yōu)化。亮度越大，吸引力越大，小的螢火蟲向亮度大的移動，移動速度與吸引力大小相關。

其中，螢火蟲相對距離為笛卡爾距離。第i只螢火蟲與第j只螢火蟲的笛卡爾距離rij為：

其中，xi為第i只螢火蟲空間坐標，d為螢火蟲所在空間維數(shù)。若第i只亮度高于第j只，則第i只螢火蟲對第j只螢火蟲的吸引力βj為：

其中，β0為最大吸引力，γ為光吸收強度系數(shù)。第j只螢火蟲向第i只螢火蟲移動時位置更新為［19］：

其中，t為算法迭代次數(shù)，a為步長因子，ξj為標準正態(tài)分布U（0，1）的值。

2.2 改進方法

標準的FA 在實際應用中存在以下問題：

（1）如果多維空間中各維度取值范圍和求解精度不同而依舊使用統(tǒng)一的步長因子，會導致各維度無法在同一迭代次數(shù)達到最優(yōu)解。

（2）隨著算法迭代次數(shù)增加，螢火蟲個體聚集于最優(yōu)解附近時，如果使用固定步長因子來更新位置會因為步長因子過大使得螢火蟲移動距離過大，導致螢火蟲在最優(yōu)解位置附近振蕩［20］，無法收斂解出最優(yōu)解，影響算法的效率和求解精度。

針對以上問題，本文在標準FA 上進行改進。

對于問題（1），在d維空間中，通過引入a→=(a1,a2,a3...ad)取代a，稱為步長向量。根據(jù)不同維度的范圍和精度設置不同的步長向量。

對于問題（2），本研究在標準螢火蟲算法基礎上，隨著螢火蟲逐漸接近最優(yōu)解逐步縮小步長因子。開始迭代時，將a設為較大值，執(zhí)行迭代時步長因子調整為：

這樣，步長因子會根據(jù)螢火蟲距離的縮短而減小。當算法接近最優(yōu)解時能保持有效的收斂速度，這種改進的螢火蟲算法簡稱為MFA（Modified Firefly Algorithm）。

2.3 優(yōu)化船舶橫搖曲線參數(shù)

式（6）中經(jīng)參數(shù)簡化后的船舶橫搖運動方程有θm、μθ、nθ三個未知參數(shù)，要獲得未知參數(shù)的最優(yōu)解，本文將擬合結果與實際數(shù)據(jù)的殘差平方和最小作為評估標準，即：

式中，θj為橫搖傾角實際值，為根據(jù)擬合獲得的橫搖傾角，m為橫搖傾角序列數(shù)據(jù)量。

利用MFA 求解船舶橫搖參數(shù)的主要步驟和算法流程如圖1 所示。

流程如下：①設置FA 基本參數(shù)，包括螢火蟲數(shù)目、步長因子a初始值、光吸收強度系數(shù)γ、最大吸引度β0、最大迭代次數(shù)；②初始化螢火蟲位置X={x1,x2,x3...xn}；③導入經(jīng)處理后的橫搖傾角數(shù)據(jù){θ1,θ2,θ3...θj}；④根據(jù)式（8）計算出螢火蟲間的吸引度；⑤根據(jù)式（11）更新步長因子；⑥根據(jù)式（9）更新螢火蟲位置；⑦判斷是否達到最大迭代次數(shù)，或殘差平方和小于設定值，如果是，執(zhí)行下一步，否則返回步驟④；⑧輸出θm、μθ、nθ最優(yōu)解。

Fig.1 Flow of the proposed MFA圖1 MFA 算法流程

3 實例應用及結果分析

3.1 準備工作

數(shù)據(jù)采集系統(tǒng)結構如圖2（彩圖掃OSID 碼可見，下同）所示，以某集裝箱船舶為例，使用兩軸傾角傳感器采集船舶傾角信號，通過ZigBee 無線模塊將采集信號傳輸?shù)缴衔粰C。上位機軟件使用LABVIEW 編譯環(huán)境。船舶基本信息如表1 所示。

Fig.2 Schematic diagram of the data acquisition system圖2 數(shù)據(jù)采集系統(tǒng)結構

Table 1 Basic information of one ship表1 某船舶基本信息 （m）

實驗前需要將兩軸動態(tài)傾角傳感器進行水平校準，然后將其安裝在吊具上，當?shù)蹙咧鋾r開始采集傾角信號，起箱時結束采集，所得傾角信號通過ZigBee 無線模塊發(fā)送至上位機進行處理，系統(tǒng)主程序采用LabVIEW 編程環(huán)境開發(fā)，參數(shù)優(yōu)化程序采用MATLAB 腳本節(jié)點編寫。

本文實驗仿真內容如下：

（1）利用標準螢火蟲算法FA、MFA 分別優(yōu)化不同裝載狀態(tài)下的橫搖傾角數(shù)據(jù)，證明利用參數(shù)優(yōu)化思想分析船舶橫搖的可行性。

（2）對兩種算法進行50 次優(yōu)化，分析優(yōu)化結果和算法誤差，對比分析MFA 算法的優(yōu)越性。

（3）改變橫搖傾角數(shù)據(jù)量，分析不同數(shù)據(jù)量下MFA 優(yōu)化結果精度，分析最優(yōu)數(shù)據(jù)量。

3.2 算法可行性分析

在進行算法計算之前需設定初始參數(shù)，本實驗中FA 和MFA 算法參數(shù)如表2、表3 所示。

Table 2 Parameters settings of MFA表2 MFA 算法參數(shù)設置

Table 3 Parameters settings of FA表3 FA 算法參數(shù)設置

從圖3、圖4 可以看出，F(xiàn)A 和MFA 兩種算法對不同裝載狀態(tài)下的船舶橫搖運動均能實現(xiàn)曲線回歸。半載和滿載狀態(tài)下船舶橫搖呈現(xiàn)為衰減的正弦曲線，符合船舶橫搖運動方程的描述，且滿載時橫搖周期較長，符合船舶橫搖運動原理。當FA 和MFA 算法迭代次數(shù)相同時，可以很明顯地看出MFA 擬合曲線與原數(shù)據(jù)更為貼近，F(xiàn)A 優(yōu)化曲線相較于真實值有明顯的偏移現(xiàn)象。

Fig.3 Fitting results of normal FA圖3 標準FA 算法擬合結果

Fig.4 Fitting results of MFA圖4 MFA 算法擬合結果

由此可以判斷，利用參數(shù)優(yōu)化思想結合船舶橫搖運動方程來分析船舶橫搖運動的方法是可行的。利用一定長度的橫搖傾角數(shù)據(jù)訓練船舶橫搖運動模型，當算法收斂到最佳擬合時可求解到橫搖運動曲線參數(shù)。同等迭代次數(shù)下對比FA 和MFA 結果，可以初步得出解算精度與使用的算法和算法參數(shù)相關的結論。

3.3 MFA 優(yōu)越性分析

螢火蟲算法優(yōu)化精度與迭代次數(shù)相關，隨著迭代次數(shù)的增加，解算結果收斂于最佳值，優(yōu)化精度也隨之提升。為分析MFA 相較于標準FA 的優(yōu)越性，本實驗基于3.2 節(jié)設定的參數(shù)進行50 次優(yōu)化，分析參數(shù)優(yōu)化均值和相對誤差。

表4 是兩組數(shù)據(jù)50 次優(yōu)化結果的平均值。表4 顯示，橫搖運動傾角峰值和無因次衰減系數(shù)與裝載狀態(tài)相關性不大，而振動周期隨著船舶裝滿而變大，符合船舶橫搖基本原理。當兩種算法迭代次數(shù)相同時，MFA 算法計算精度優(yōu)勢明顯。其中，半載時擬合曲線與原數(shù)據(jù)決定系數(shù)為0.972 2，滿載時決定系數(shù)為0.984 1，表明擬合結果具有極高的可靠性。而標準的FA 最大相對誤差和平均相對誤差都大于MFA，決定系數(shù)也相對較低，說明標準FA 在同等迭代步數(shù)下的優(yōu)化結果精度比MFA 更低。實際上增加迭代步數(shù)同樣可以使FA 的優(yōu)化結果達到MFA 的效果，不過這需要更大的計算量，由此對比體現(xiàn)出MFA 的優(yōu)越性。

Table 4 Optimization results of MFA and FA parameters under different loading conditions表4 不同裝載狀態(tài)MFA、FA 參數(shù)優(yōu)化結果

3.4 實測數(shù)據(jù)量對優(yōu)化結果影響分析

在港口實際工作中，吊臂具著箱過程時間短暫，所采集到的實際傾角數(shù)據(jù)量有限。為實現(xiàn)在有限的實測數(shù)據(jù)內求解最優(yōu)結果并進行橫搖運動預測，需要分析實測數(shù)據(jù)量對擬合結果的影響，在盡可能少的數(shù)據(jù)量情況下擬合出可靠性高的橫搖運動曲線。

本實驗中橫搖傾角采樣周期為0.15s，船舶為滿載狀態(tài)，橫搖周期約為18.4s，為分析實測數(shù)據(jù)量對優(yōu)化結果的影響，本實驗分別對1/4、1/2、3/4、1 個周期的數(shù)據(jù)進行模型參數(shù)優(yōu)化。

圖5 顯示，當代入算法的實測數(shù)據(jù)為T/4 時，在T/4 時間后擬合結果與實測值有很大誤差，優(yōu)化模型無法準確預測橫搖運動。代入算法的數(shù)據(jù)量為T/2、3T/4、T 時，優(yōu)化模型的預測結果與實測值有極高的吻合度。結合測量船舶實際橫搖傾角的工作環(huán)境，使用T/2 的實測數(shù)據(jù)訓練模型，優(yōu)化橫搖運動模型參數(shù)最為合適，其優(yōu)化結果預測橫搖運動有很高的可靠性。

MFA 作為一種群智能算法，優(yōu)勢在于能從維度空間各個方向收斂，相對于傳統(tǒng)的優(yōu)化算法有著更快的迭代速度［21-23］，算法在迭代開始時就能自動分為多個子群，子群聚集在局部最優(yōu)解附近，通過局部最優(yōu)解的比對尋找到全局最優(yōu)解。MFA 速度快、精度高的特點適用于優(yōu)化船舶橫搖運動模型參數(shù)，并預測橫搖運動狀態(tài)。但在實例實驗中，當局部最優(yōu)解與全局最優(yōu)解位置接近時，存在求解結果為局部最優(yōu)解情況，該問題有待進一步研究。

Fig.5 Fitting results with different data volumes圖5 不同數(shù)據(jù)量的擬合結果

4 結語

本文通過解算船舶橫搖參數(shù)求解初穩(wěn)性高的方法分析系泊船舶裝卸作業(yè)工況下的穩(wěn)性問題，并針對橫搖方程提出基于改進螢火蟲算法的優(yōu)化方法。該方法提高了求解效率，改進了步長向量和步長更新算法。通過實驗驗證了該方法的可行性。對比迭代30 次情況下MFA 相較于FA的優(yōu)越性。實驗表明，MFA 算法能夠在較短時間內較好擬合不同裝載狀態(tài)下的船舶橫搖運動曲線。在半載情況下，擬合結果與原數(shù)據(jù)決定系數(shù)為0.972 2，滿載情況下為0.984 1，擬合結果具有很強的可靠性。結合實際工況分析了訓練模型的實測數(shù)據(jù)最合適的數(shù)據(jù)量為T/2。該方法將智能群優(yōu)化算法運用到船舶橫搖運動狀態(tài)預測，為系泊船舶裝卸作業(yè)的穩(wěn)性分析和預測提供了解決方案。