肖力旗，曹 民，張柳明

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

0 引言

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

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

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

針對(duì)系泊船舶裝卸作業(yè)特殊的工況及穩(wěn)性監(jiān)測(cè)實(shí)時(shí)性問題，本研究通過分析橫搖周期與初穩(wěn)性高的關(guān)系，利用改進(jìn)的螢火蟲算法優(yōu)化船舶橫搖運(yùn)動(dòng)方程模型，預(yù)測(cè)船舶橫搖運(yùn)動(dòng)狀態(tài)，計(jì)算出船舶初穩(wěn)性高，從而分析系泊船舶穩(wěn)性。采用仿真實(shí)例驗(yà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)性高不應(yīng)低于0.15m［13］。船舶在靜水中的初穩(wěn)性高與橫搖周期存在如下關(guān)系［14］：

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

1.2 船舶橫搖運(yùn)動(dòng)方程

對(duì)于船舶的橫搖運(yùn)動(dòng)，若不考慮風(fēng)浪流等外界干擾，船舶的橫搖運(yùn)動(dòng)力矩方程［15］為：

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

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

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

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

1.3 傾角序列前期處理

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

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

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

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

2 改進(jìn)螢火蟲算法

2.1 標(biāo)準(zhǔn)螢火蟲算法

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

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

（2）相對(duì)亮度與距離平方成反比，當(dāng)距離過大時(shí)無吸引力。

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

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

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

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

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

2.2 改進(jìn)方法

標(biāo)準(zhǔn)的FA 在實(shí)際應(yīng)用中存在以下問題：

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

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

針對(duì)以上問題，本文在標(biāo)準(zhǔn)FA 上進(jìn)行改進(jìn)。

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

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

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

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

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

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

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

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

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

3 實(shí)例應(yīng)用及結(jié)果分析

3.1 準(zhǔn)備工作

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

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

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

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

本文實(shí)驗(yàn)仿真內(nèi)容如下：

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

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

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

3.2 算法可行性分析

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

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

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

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

Fig.3 Fitting results of normal FA圖3 標(biāo)準(zhǔn)FA 算法擬合結(jié)果

Fig.4 Fitting results of MFA圖4 MFA 算法擬合結(jié)果

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

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

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

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

3.4 實(shí)測(cè)數(shù)據(jù)量對(duì)優(yōu)化結(jié)果影響分析

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

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

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

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

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

4 結(jié)語

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