周志剛，黃建偉，董根金

（武漢第二船舶設(shè)計(jì)研究所，湖北武漢 430064）

0 引言

潛艇在應(yīng)急上浮過(guò)程中，容易出現(xiàn)大攻角運(yùn)動(dòng)狀態(tài)。特別是當(dāng)艙室進(jìn)水和舵卡事故，原重力與浮力平衡、縱傾力矩平衡都受到破壞時(shí)，潛艇垂向速度與前進(jìn)速度相比并非小量，艇體運(yùn)動(dòng)相對(duì)水流形成大攻角。大攻角情況下的潛艇水動(dòng)力非常復(fù)雜，表現(xiàn)出強(qiáng)烈的非線(xiàn)性特征?？紤]到垂直面的運(yùn)動(dòng)參數(shù)是造成潛艇大機(jī)動(dòng)水動(dòng)力非線(xiàn)性的主要原因，因此可以將大攻角所產(chǎn)生的非線(xiàn)性水動(dòng)力看作是潛艇垂向速度w和縱傾角速度q的函數(shù)。當(dāng)前國(guó)內(nèi)相關(guān)研究［1］的思路是采用低階多項(xiàng)式或樣條函數(shù)對(duì)潛艇大攻角機(jī)動(dòng)運(yùn)動(dòng)水動(dòng)力進(jìn)行擬合，但其擬合精度或?qū)嵱眯缘确矫嫱⒉荒茏屓藵M(mǎn)意。采用合適的方法準(zhǔn)確表達(dá)潛艇大攻角機(jī)動(dòng)運(yùn)動(dòng)的水動(dòng)力，建立用于預(yù)報(bào)其大攻角機(jī)動(dòng)運(yùn)動(dòng)狀態(tài)變化的數(shù)學(xué)模型，對(duì)潛艇的動(dòng)力抗沉等有重大的指導(dǎo)意義。

本文根據(jù)已有的潛艇大攻角機(jī)動(dòng)運(yùn)動(dòng)試驗(yàn)資料，結(jié)合多項(xiàng)式和樣條函數(shù)的方法擬合其水動(dòng)力的結(jié)果，指出了各自的局限性，設(shè)計(jì)了以攻角α（w'=sinα）和縱傾角速度q'作為輸入量，縱向力X'、垂向力Z'、縱傾力矩M'分別作為輸出量的3層BP神經(jīng)網(wǎng)絡(luò)，應(yīng)用MATLAB 7.1神經(jīng)網(wǎng)絡(luò)工具箱對(duì)這一系列試驗(yàn)數(shù)據(jù)的水動(dòng)力進(jìn)行了有效的擬合。

1 已有的潛艇大機(jī)動(dòng)水動(dòng)力表達(dá)形式

1.1 已有潛艇大機(jī)動(dòng)水動(dòng)力試驗(yàn)數(shù)據(jù)特點(diǎn)

潛艇大攻角狀態(tài)的水動(dòng)力是通過(guò)約束模擬試驗(yàn)獲得的無(wú)因次量，包括2個(gè)運(yùn)動(dòng)狀態(tài)參數(shù)α和q'，水動(dòng)力則為X'，Z'和M'，其中攻角α的范圍達(dá)到了±90°。將X'，Z'和M'的試驗(yàn)點(diǎn)在空間展開(kāi)，得到的是各力（矩）隨攻角和縱傾角速度變化的曲面，發(fā)現(xiàn)大攻角情況下的水動(dòng)力隨α的變化表現(xiàn)出非常強(qiáng)烈的非線(xiàn)性，而隨q'的變化則近似為線(xiàn)性關(guān)系。如圖1所示，（a）圖為某一q'值下的縱傾力矩M'隨攻角α變化的曲線(xiàn);（b）圖為某一攻角α值下的縱傾力矩M'隨q'變化的曲線(xiàn)。由圖還可以看到，試驗(yàn)數(shù)據(jù)并非是均勻分布的，在水動(dòng)力變化較大的一些攻角區(qū)間內(nèi)試驗(yàn)數(shù)據(jù)相對(duì)很少;而水動(dòng)力在隨q'變化的曲線(xiàn)上，由于受到試驗(yàn)條件的限制，有一大段區(qū)間內(nèi)的水動(dòng)力值無(wú)法通過(guò)試驗(yàn)獲得。

圖1 縱傾力矩M'試驗(yàn)數(shù)據(jù)特點(diǎn)Fig.1Characteristic of M'testing data

1.2 多項(xiàng)式的水動(dòng)力表達(dá)形式

根據(jù)相關(guān)資料的研究分析，在潛艇大攻角機(jī)動(dòng)范圍內(nèi)，宜采用二階多項(xiàng)式表達(dá)其水動(dòng)力，以避免由于項(xiàng)數(shù)增加導(dǎo)致參數(shù)估計(jì)時(shí)的較大協(xié)方差，以及高階水動(dòng)力導(dǎo)數(shù)的求取;而且這種表達(dá)形式與潛艇標(biāo)準(zhǔn)空間運(yùn)動(dòng)方程中水動(dòng)力的表達(dá)形式保持一致，能從本質(zhì)上反映潛艇水動(dòng)力的物理特性。圖2給出了在攻角±90°范圍內(nèi)，經(jīng)這種方法處理后在某一縱傾角速度下縱傾力矩隨攻角變化的擬合曲線(xiàn)。

圖2 縱傾力矩M'的多項(xiàng)式表達(dá)擬合曲線(xiàn)Fig.2Fitted curf of M'by polynomial

如圖2所示，多項(xiàng)式表達(dá)的水動(dòng)力只能大致反映其在整個(gè)攻角范圍內(nèi)的變化趨勢(shì)，但是其擬合精度顯然不足以建立準(zhǔn)確的數(shù)學(xué)仿真模型，從而滿(mǎn)足實(shí)艇操縱性預(yù)報(bào)的要求。如果采取分段低階多項(xiàng)式的形式表達(dá)水動(dòng)力，雖然擬合精度可以提高一些，但在各區(qū)間段的端點(diǎn)不能保持光滑過(guò)渡，處理試驗(yàn)數(shù)據(jù)時(shí)也會(huì)帶來(lái)相當(dāng)?shù)穆闊?。由此可?jiàn)，低階多項(xiàng)式形式不適用于在整個(gè)攻角范圍內(nèi)表達(dá)潛艇大機(jī)動(dòng)水動(dòng)力。

1.3 樣條函數(shù)的水動(dòng)力表達(dá)形式

利用樣條函數(shù)建立非線(xiàn)性水動(dòng)力數(shù)學(xué)模型的優(yōu)點(diǎn)是，可以用低階多項(xiàng)式很好地逼近各種非線(xiàn)性。樣條函數(shù)模型也是基于Bryan準(zhǔn)定常假設(shè)的飛機(jī)大攻角氣動(dòng)力常用代數(shù)模型之一。在大攻角下，飛機(jī)的縱向氣動(dòng)力系數(shù)主要是攻角α的非線(xiàn)性函數(shù)，只需采用單變量的樣條函數(shù)模型;橫向氣動(dòng)系數(shù)對(duì)攻角α和側(cè)滑角β都是非線(xiàn)性的，必須采用雙變量的樣條函數(shù)模型［4］。從前面試驗(yàn)數(shù)據(jù)的分析知道，潛艇垂直面大攻角機(jī)動(dòng)運(yùn)動(dòng)水動(dòng)力的非線(xiàn)性只是體現(xiàn)在攻角α上，實(shí)踐證明采用3次樣條函數(shù)插值即可滿(mǎn)足工程應(yīng)用要求，保證了擬合出來(lái)的水動(dòng)力在試驗(yàn)點(diǎn)處二階連續(xù)。在區(qū)間［α1，αn］上關(guān)于變量α的3次樣條函數(shù)［5］可以表示為:

樣條節(jié)點(diǎn)α1，α2，…，αn滿(mǎn)足α1＜α2＜…＜αn的順序。Ch為樣條的常數(shù)項(xiàng)。圖3給出的是基于該函數(shù)形式擬合的某一縱傾角速度下縱傾力矩隨攻角變化的擬合曲線(xiàn)。

如圖3所示，樣條函數(shù)擬合的水動(dòng)力嚴(yán)格地經(jīng)過(guò)了每1個(gè)試驗(yàn)點(diǎn)，因而具有較高的精度。但這種高精度是建立在水動(dòng)力試驗(yàn)數(shù)據(jù)的高準(zhǔn)確度上，樣條函數(shù)形式并不利于消除試驗(yàn)數(shù)據(jù)本身的誤差，且表達(dá)完整的水動(dòng)力情況需要大量的試驗(yàn)數(shù)據(jù)，一旦數(shù)據(jù)偏少或分布極其不均，擬合的曲線(xiàn)可能會(huì)出現(xiàn)不必要的“振蕩”，表現(xiàn)在空間上就是一些凸凹不平的坑，因而采用樣條函數(shù)形式表達(dá)潛艇大機(jī)動(dòng)水動(dòng)力也有局限性。

圖3 縱傾力矩M'的樣條函數(shù)表達(dá)擬合曲線(xiàn)Fig.3Fitted curf of M'by spline interpolation

綜上所述，如何應(yīng)用合理的表達(dá)方式在樣本數(shù)據(jù)較少的情況下準(zhǔn)確預(yù)報(bào)潛艇大機(jī)動(dòng)水動(dòng)力，以及保證其在沒(méi)有試驗(yàn)數(shù)據(jù)的區(qū)間具有較好的泛化能力，并且能推廣到大機(jī)動(dòng)運(yùn)動(dòng)整個(gè)攻角區(qū)間，具有非常重要的工程應(yīng)用價(jià)值，值得進(jìn)一步研究。

2 潛艇大機(jī)動(dòng)水動(dòng)力表達(dá)

2.1 BP神經(jīng)網(wǎng)絡(luò)原理及構(gòu)建過(guò)程

BP網(wǎng)絡(luò)是1種單向傳播的多層前向網(wǎng)絡(luò)，本文采用的3層BP網(wǎng)絡(luò)結(jié)構(gòu)如圖4所示。輸入層為2個(gè)神經(jīng)元，輸出層為1個(gè)神經(jīng)元，設(shè)隱層神經(jīng)元數(shù)為N，隱層傳遞函數(shù)采用f1（·），輸出層傳遞函數(shù)采用f2（·），誤差函數(shù)采用E（·），則網(wǎng)絡(luò)的學(xué)習(xí)過(guò)程可用下列函數(shù)描述:

圖43 層BP網(wǎng)絡(luò)結(jié)構(gòu)圖Fig.4Frame of BP neural net with three layers

學(xué)習(xí)過(guò)程中將誤差反向傳播計(jì)算，對(duì)網(wǎng)絡(luò)各相鄰節(jié)點(diǎn)的連接權(quán)值w和偏置θ進(jìn)行調(diào)整，使誤差函數(shù)E（·）沿負(fù)梯度方向下降，然后再轉(zhuǎn)入正向傳播過(guò)程，反復(fù)迭代，直到誤差達(dá)到穩(wěn)態(tài)或小于給定的值為止。

當(dāng)輸入、輸出變量以及隱層層數(shù)明確的情況下，構(gòu)建BP神經(jīng)網(wǎng)絡(luò)的一般步驟如下:

1）對(duì)樣本的歸一化處理，以消除各變量的量綱影響，MATLAB中命令premnmx可將樣本數(shù)據(jù)處理到區(qū)間（－1，1）上。

2）根據(jù)神經(jīng)網(wǎng)絡(luò)的用途選擇相應(yīng)的輸出層傳遞函數(shù)，根據(jù)數(shù)據(jù)的特點(diǎn)選擇合適的隱層傳遞函數(shù)。

3）設(shè)計(jì)隱層最佳單元數(shù)，一般是在經(jīng)驗(yàn)公式范圍內(nèi)采用試湊法尋找［5］。隱層神經(jīng)元的個(gè)數(shù)是神經(jīng)網(wǎng)絡(luò)訓(xùn)練成敗的關(guān)鍵，不合理的單元個(gè)數(shù)容易出現(xiàn)“欠擬合”或“過(guò)擬合”的現(xiàn)象。

4）選擇合適的訓(xùn)練函數(shù)，MATLAB提供的BP神經(jīng)網(wǎng)絡(luò)的訓(xùn)練函數(shù)有數(shù)十種。

2.2 建立潛艇大機(jī)動(dòng)水動(dòng)力表達(dá)的BP神經(jīng)網(wǎng)絡(luò)

1）樣本預(yù)處理與傳遞函數(shù)選擇本文采用更一般的轉(zhuǎn)換公式

隱層傳遞函數(shù)取反對(duì)稱(chēng)形式的雙曲正切函數(shù)，輸出層傳遞函數(shù)取純線(xiàn)性函數(shù)。

2）隱層神經(jīng)元數(shù)的確定

對(duì)于本組潛艇大機(jī)動(dòng)的試驗(yàn)數(shù)據(jù)，采用經(jīng)驗(yàn)公式設(shè)計(jì)的隱層單元數(shù)，其建立的BP神經(jīng)網(wǎng)絡(luò)模型并不能獲得滿(mǎn)意的效果。最佳神經(jīng)元個(gè)數(shù)的取值范圍必將擴(kuò)大，盲目摸索非常耗費(fèi)人力和時(shí)間。本文采用近似于窮舉法來(lái)選擇隱層神經(jīng)元的個(gè)數(shù)，即先確定1個(gè)較大的隱層神經(jīng)元數(shù)目的范圍，再?gòu)妮^少神經(jīng)元開(kāi)始訓(xùn)練，然后逐漸增加神經(jīng)元的個(gè)數(shù)，重復(fù)學(xué)習(xí)訓(xùn)練。這個(gè)過(guò)程非常類(lèi)似于人腦對(duì)事物認(rèn)知從簡(jiǎn)單到復(fù)雜的情況，每一步的訓(xùn)練過(guò)程都以上一步的訓(xùn)練結(jié)果為初值，相當(dāng)于對(duì)上次學(xué)習(xí)有所“記憶”，整個(gè)確定最佳隱層神經(jīng)元個(gè)數(shù)的過(guò)程也可以看作是循環(huán)強(qiáng)化學(xué)習(xí)的過(guò)程，這樣不僅能完全自適應(yīng)地找到最佳隱層單元數(shù)，而且還獲得了非常高的訓(xùn)練精度。這一點(diǎn)可以對(duì)比應(yīng)用最佳隱層只進(jìn)行1次訓(xùn)練的擬合曲線(xiàn)和循環(huán)訓(xùn)練后的擬合曲線(xiàn)來(lái)驗(yàn)證，如圖5所示。

圖5 不同訓(xùn)練次數(shù)擬合曲線(xiàn)對(duì)比Fig.5Compare of fitted curf by different train time

考慮到潛艇在整個(gè)大攻角范圍內(nèi)機(jī)動(dòng)運(yùn)動(dòng)時(shí)的實(shí)際應(yīng)用情形，為了保證網(wǎng)絡(luò)在正常機(jī)動(dòng)和大攻角機(jī)動(dòng)狀態(tài)下均具有良好的適用性，本文還設(shè)計(jì)了2個(gè)均方差性能參數(shù)mse±30°和mse±90°，以及擬合水動(dòng)力隨攻角變化曲線(xiàn)的拐點(diǎn)個(gè)數(shù)來(lái)控制循環(huán)次數(shù)，一旦同時(shí)達(dá)到這3個(gè)參數(shù)設(shè)定的值即認(rèn)為網(wǎng)絡(luò)已經(jīng)訓(xùn)練成熟，其中，mse±30°和mse±90°分別用于表征（－30°，30°）和（－90°，90°）區(qū)間內(nèi)的網(wǎng)絡(luò)預(yù)報(bào)值與試驗(yàn)值的均方差。

3）訓(xùn)練函數(shù)的確定及其預(yù)測(cè)誤差

BP神經(jīng)網(wǎng)絡(luò)訓(xùn)練函數(shù)實(shí)質(zhì)上基于最速下降反向傳播算法經(jīng)過(guò)改進(jìn)衍化而成。例如，在誤差反向傳播的基礎(chǔ)上，將每1個(gè)權(quán)值的變化加上1項(xiàng)正比于前次權(quán)值變化量的值，并自適應(yīng)學(xué)習(xí)速率的訓(xùn)練函數(shù)traingdx;基于數(shù)值優(yōu)化方法的算法，結(jié)合了梯度下降法和牛頓法在誤差下降搜索方向的優(yōu)點(diǎn)，由Levenberg和Marquardt提出的訓(xùn)練函數(shù)trainlm;同樣以L(fǎng)-M優(yōu)化理論為基礎(chǔ)對(duì)網(wǎng)絡(luò)權(quán)值和偏置進(jìn)行調(diào)整，而將誤差性能函數(shù)mse修改為包含有網(wǎng)絡(luò)權(quán)值和偏置均方值的誤差性能函數(shù)msereg的Bayesian正則化訓(xùn)練函數(shù)trainbr，其中:

實(shí)驗(yàn)結(jié)果顯示，網(wǎng)絡(luò)選擇這3種訓(xùn)練函數(shù)相對(duì)于其他訓(xùn)練函數(shù)在收斂精度上具有絕對(duì)的優(yōu)勢(shì)，本文將相應(yīng)的訓(xùn)練后最佳隱層神經(jīng)元個(gè)數(shù)N和設(shè)定的均方差性能參數(shù)mse±30°和mse±90°列于表1中。

從表1可見(jiàn)，在網(wǎng)絡(luò)規(guī)模相差不大的情況下，Bayesian正則化方法訓(xùn)練的BP神經(jīng)網(wǎng)絡(luò)相對(duì)于其他2種方法在擬合精度上幾乎提高了1個(gè)數(shù)量級(jí)，而且分析其誤差性能函數(shù)msereg可知，這種方法本身能使網(wǎng)絡(luò)獲得較小的權(quán)值和偏置，使得網(wǎng)絡(luò)的響應(yīng)趨于平滑，減小了過(guò)度擬合的可能性。圖6給出了采用函數(shù)trainbr訓(xùn)練后縱向力、垂向力和縱傾力矩隨攻角和縱傾角速度變化的曲面。

圖6 Bayesian正則化法訓(xùn)練的BP神經(jīng)網(wǎng)絡(luò)擬合曲面Fig.6Fitted curf of hydrodynamics based on Bayesian-BP neural network

可以看到，所建立的網(wǎng)絡(luò)在整個(gè)潛艇大機(jī)動(dòng)范圍內(nèi)顯示了良好的推廣性，沒(méi)有出現(xiàn)類(lèi)似樣條函數(shù)表達(dá)形式獲得曲面的凸凹不平的坑，而且神經(jīng)網(wǎng)絡(luò)模型計(jì)算值與數(shù)據(jù)樣本的相關(guān)系數(shù)分別達(dá)到了RX=0.995 8，RZ=0.999 6，RM=0.998 9。

2.3 潛艇大機(jī)動(dòng)水動(dòng)力不同表達(dá)方式比較

圖7所示從上到下依次為樣條函數(shù)、Bayesian正則化法BP神經(jīng)網(wǎng)絡(luò)和多項(xiàng)式擬合的縱向力隨攻角和縱傾角速度變化曲面。這里3個(gè)曲面在垂直方向的距離是為了便于比較圖形而設(shè)置的，并不代表任意兩者之間存在這樣一個(gè)差值。由圖示效果可以直觀(guān)得出如下結(jié)論:Bayesian正則化法BP神經(jīng)網(wǎng)絡(luò)對(duì)數(shù)據(jù)樣本集較小的潛艇大機(jī)動(dòng)水動(dòng)力函數(shù)逼近能力優(yōu)異，是一種優(yōu)于多項(xiàng)式和樣條函數(shù)形式表達(dá)潛艇大機(jī)動(dòng)運(yùn)動(dòng)水動(dòng)力的方法。

圖7 3種表達(dá)方式對(duì)縱向力的擬合曲面對(duì)比圖Fig.7Compare of X'fitted curf by three means

3 結(jié)語(yǔ)

1）潛艇的大攻角機(jī)動(dòng)運(yùn)動(dòng)表現(xiàn)出強(qiáng)烈的非線(xiàn)性特征，采用多項(xiàng)式、樣條函數(shù)形式表達(dá)其水動(dòng)力具有各自明顯的局限性或不足，本文應(yīng)用BP神經(jīng)網(wǎng)絡(luò)模型擬合了整個(gè)大攻角范圍內(nèi)的水動(dòng)力試驗(yàn)值，展示了其作為一種函數(shù)逼近工具的優(yōu)越性;

2）本文建立的BP神經(jīng)網(wǎng)絡(luò)采取了循環(huán)強(qiáng)化訓(xùn)練的方法，以及通過(guò)近似窮舉法設(shè)計(jì)神經(jīng)網(wǎng)絡(luò)隱層單元數(shù)，實(shí)踐表明這些處理方法的方便、可行性;在確定訓(xùn)練函數(shù)過(guò)程中，發(fā)現(xiàn)Bayesian正則化方法具有較高的網(wǎng)絡(luò)預(yù)報(bào)精度和網(wǎng)絡(luò)泛化能力，適合于樣本數(shù)據(jù)集較小的潛艇大機(jī)動(dòng)水動(dòng)力擬合;

3）本文采用BP神經(jīng)網(wǎng)絡(luò)技術(shù)表達(dá)潛艇大攻角機(jī)動(dòng)運(yùn)動(dòng)水動(dòng)力的研究工作，為合理應(yīng)用水動(dòng)力試驗(yàn)數(shù)據(jù)并建立其數(shù)學(xué)仿真模型提供了一種新的思路，具有廣闊的工程應(yīng)用前景。

［1］張華，吳寶山．扁平型潛水器空間大機(jī)動(dòng)數(shù)學(xué)模型研究［A］．2004年船舶水動(dòng)力學(xué)學(xué)術(shù)會(huì)議論文集［C］．2004．122－129．

［2］田雨波．混合神經(jīng)網(wǎng)絡(luò)技術(shù)［M］．北京:科學(xué)出版社，2009．38－48．

［3］從爽．面向MATLAB工具箱的神經(jīng)網(wǎng)絡(luò)理論與研究［M］．合肥:中國(guó)科學(xué)技術(shù)大學(xué)出版社，2009．63－100．

［4］汪清，蔡金獅．飛機(jī)大攻角氣動(dòng)力建模研究進(jìn)展［J］．氣動(dòng)實(shí)驗(yàn)與測(cè)量控制，1994，（9）:7－16．

［5］李紅．?dāng)?shù)值分析［M］．武漢:華中科技大學(xué)出版社，2006．40－49．

［6］飛思科技產(chǎn)品研發(fā)中心．神經(jīng)網(wǎng)絡(luò)理論與MATLAB7實(shí)現(xiàn)［M］．北京:電子工業(yè)出版社，2005．99－104．