基于樣本船的興波阻力智能預(yù)報方法

陳 帥，周廣利，王 超，蔣彩霞

（1.中國船舶科學(xué)研究中心，江蘇 無錫 214082；2.哈爾濱工程大學(xué) 船舶工程學(xué)院，哈爾濱 150001）

0 引 言

為了提高船型優(yōu)化設(shè)計系統(tǒng)的計算效率，近似模型的引入至關(guān)重要。目前以減小興波阻力性能為優(yōu)化目標的船型設(shè)計問題較多，而船體曲面變形技術(shù)和興波阻力預(yù)報技術(shù)為近似模型訓(xùn)練提供樣本數(shù)據(jù)，隨著計算機和人工智能技術(shù)的迅速發(fā)展，構(gòu)建基于樣本船型的興波阻力高精度近似計算模型來解決優(yōu)化過程中的海量計算具有重要研究意義。

船體曲面變形技術(shù)需要解決的是如何以較少的設(shè)計變量對船體曲面進行變換，本文根據(jù)自主開發(fā)的NURBS技術(shù)[1]，采用改進平移法和自由變形技術(shù)對船體曲面進行修改生成海量的船型樣本方案。船體興波阻力快速預(yù)報可以通過求解定常興波問題得到，在上兩個世紀中，出現(xiàn)了一些經(jīng)典的興波勢流理論，如Michell 理論、Kelvin 源（或Havelock 源）、Neumann-Kelvin（NK）理論和新細長體理論等，這些理論為后續(xù)的興波問題研究奠定了基礎(chǔ)。Dawson[2]在自由面和船體表面布置Rankine 源來求解定常興波問題，使得Rankine 源法在船舶興波阻力預(yù)報上取得了巨大的成功。本文開發(fā)了基于貼體網(wǎng)格[3]的Dawson方法計算程序，只需數(shù)十秒就可以計算一次興波，適用于樣本船型預(yù)報。

在近似模型方面，馮佰威[4-5]綜述了近似技術(shù)在船型自動優(yōu)化中的意義并進行應(yīng)用；萬德成等[6]討論了近似模型在船型優(yōu)化中的應(yīng)用現(xiàn)狀以及研宄難點；劉志強[7]以O(shè)PTShip-SJTU[12]求解器為基礎(chǔ)，分析了不同近似模型的計算精度；Huang[8]提出了一種基于徑向基函數(shù)的近似模型降低了貨船優(yōu)化的計算成本；Vesting[9]采用近似模型與局部搜索方法相結(jié)合大幅度提高了船用螺旋槳的設(shè)計效率；Priftis[10]采用近似模型研究了不確定條件下船舶的參數(shù)化設(shè)計和多目標優(yōu)化問題，在設(shè)計初期對總阻力、運費率和成本進行優(yōu)化；Kim[11]采用響應(yīng)面Kriging模型優(yōu)化Series-60船型，結(jié)果表明響應(yīng)面法大大降低計算成本。為了實現(xiàn)興波阻力的智能預(yù)報，本文分別對樣本船型生成、興波阻力預(yù)報和近似模型構(gòu)建展開研究，首先提出曲面變形技術(shù)的理論方法，以ITM 和FFD 兩種方法對船體變形進行應(yīng)用。其次，給出興波理論和數(shù)值計算方法，驗證本文興波阻力計算程序的可靠性。最后，建立基于SLE-BP 神經(jīng)網(wǎng)絡(luò)的興波阻力預(yù)報近似模型，利用KCS船型進行智能預(yù)報方法的應(yīng)用和驗證。

1 曲面變形技術(shù)

曲面變形技術(shù)是實現(xiàn)船型優(yōu)化設(shè)計的基礎(chǔ)，同時為近似模型訓(xùn)練提供樣本船型。本文根據(jù)NURBS 技術(shù)[1]，通過修改曲面控制網(wǎng)格可以實現(xiàn)曲面變形，該方法使得變形后的曲面仍為NURBS 曲面。如圖1所示，在船體曲面生成的過程中，采用ITM和FFD法對曲線控制頂點進行修改，通過變形后的控制網(wǎng)格和基函數(shù)計算出變形后的船體曲面。

圖1 曲面變形技術(shù)的程序?qū)崿F(xiàn)流程Fig.1 Program implementation flow of surface deformation technology

1.1 改進平移法

根據(jù)修改函數(shù)式（1），以曲面v方向某曲線的變形為例，ITM 法原理示意圖如圖2 所示，f1(x)和f2(x)為變形前后的函數(shù)曲線，g(x)為修改函數(shù)曲線。針對船體曲面變形，以控制網(wǎng)格點作為輸入即可得到變形中的位移大小和方向，再通過曲面重構(gòu)就可以得到變形船體。以KCS船體曲面變形為例，從圖3可以看出，式（1）中x1、x2、α1和α2四個參數(shù)含義均有體現(xiàn)，給定參數(shù)擾動區(qū)間即可實現(xiàn)船體曲面的整體變形。ITM 法的特點是設(shè)計變量涉及較少，非常利于船體幾何（半船、整船）的修改，并且基于NURBS技術(shù)的船體曲面控制網(wǎng)格變換可以保證新船型的光順性和連續(xù)性。

圖2 ITM法原理示意圖Fig.2 Schematic diagram of ITM method

圖3 KCS船體曲面變形示意圖Fig.3 Deformation diagram of KCS hull surface

1.2 自由變形方法

自由變形（free-form deformation,FFD）方法[3]是Sederberg 和Parry 提出的一種網(wǎng)格自由變形方法。它作為常用的幾何造型技術(shù)，目前被廣泛應(yīng)用于計算機動畫、機器人導(dǎo)航和船體幾何重構(gòu)等領(lǐng)域。

在一個含有待變形物體的長方體中構(gòu)造局部坐標系O'-STU，在笛卡爾坐標系O-XYZ下的任意點X在局部坐標系中的坐標為(s,t,u)，則有

式中，B代表Bernstein多項式基函數(shù)，定義為

由式（5）可知，初始船體曲面控制網(wǎng)格點是長方體內(nèi)控制頂點的線性函數(shù)。將部分控制頂點的位置作為設(shè)計變量，通過控制長方體框架的變形來達到船體曲面控制網(wǎng)格變形的目的。長方體框架內(nèi)的船體曲面控制網(wǎng)格點的局部坐標為(s,t,u)，部分控制頂點的位置變化將導(dǎo)致X移動到X'：

如圖4 所示，以KCS 球艏變形為例，在船艏處布置控制框架確定曲面變形區(qū)域，利用等分數(shù)設(shè)置框架內(nèi)控制點位置，控制點在x、y、z方向上的位移控制球艏曲面在三個方向上的變形，如圖4（a）所示，在控制框架內(nèi)的曲面控制網(wǎng)格隨著部分控制點的位移而變形，通過曲面正算可計算出曲面節(jié)點的變形，如圖4（b）所示，船艏縱剖線變形前后對比顯示，整個球艏區(qū)域變形節(jié)點保持協(xié)調(diào)性。

圖4 KCS船體球艏變形Fig.4 KCS hull bulbous bow deformation

2 興波阻力預(yù)報

2.1 興波理論

本文采用C++語言開發(fā)，基于Dawson 法的興波阻力計算步驟[2]，首先在船體表面上用面元法計算疊模擾動流，然后在船體表面和自由面上利用物面邊界條件、自由邊界條件和輻射條件計算興波擾動流。計算模型離散如圖5所示，把求解的疊模繞流和興波擾動流問題轉(zhuǎn)化為求解船體表面和自由面分布源點強度問題，有

圖5 計算模型離散示意圖Fig.5 Schematic diagram of discrete calculation model

式中：Φ為總速度勢，滿足拉普拉斯方程；φr為疊模擾動速度勢；φw為興波擾動速度勢；為疊模表面源強；σB和σF分別為船體表面和自由面源強；SB為船體表面；SF為自由面；為場點p到船體表面上點qb和自由面上點qw的距離。

2.2 數(shù)值計算驗證

以Wigley 和KCS船型作為算例，圖6～7為自由面在各方向的加密增長率示意圖，自由面采用貼體網(wǎng)格。通過計算域大小變化對結(jié)果的影響分析，確定自由面上游擴展0.5L，下游擴展2.5L，y方向延伸1.5L。如圖所示對船體首尾部及船寬方向進行了加密。

圖6 Wigley船型自由面網(wǎng)格示意圖Fig.6 Schematic diagram of Wigley ship free surface grid

圖7 KCS船體及自由面網(wǎng)格示意圖Fig.7 Schematic diagram of KCS hull and free surface grid

圖8（a）和圖9（a）為以自由面橫向網(wǎng)格數(shù)為變量的興波阻力系數(shù)收斂圖，從圖中可以看出，隨著面元數(shù)的增加興波阻力系數(shù)趨近一個值，因此Dawson方法具有較好的收斂性。圖8（b）和圖9（b）為兩種船型興波阻力系數(shù)隨不同F(xiàn)r數(shù)下航行的變化曲線，通過對比國內(nèi)外不同文獻的興波理論計算值，計算結(jié)果雖有些不同，但隨Fr數(shù)的變化較為一致，說明了本文程序計算興波阻力的可靠性。

圖8 Wigley船型興波阻力系數(shù)Fig.8 Wave-making resistance coefficient of Wigley ship form

圖9 KCS船型興波阻力系數(shù)Fig.9 Wave-making resistance coefficient of KCS ship form

3 近似模型

3.1 SLE方法

圖10 SLE方法樣本生成過程Fig.10 Sample generation process of SLE method

圖11 試驗樣本空間填充性能對比Fig.11 Comparison of filling performance of test sample space

3.2 BP神經(jīng)網(wǎng)絡(luò)

圖12 為BP 神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，包含輸入層、隱含層和輸出層。輸入層對當(dāng)前輸入的數(shù)據(jù)進行處理，隱含層把輸入層的數(shù)據(jù)處理好后傳向輸出層，輸出層給出當(dāng)此正向傳播的輸出值和訓(xùn)練誤差。步驟如下：

圖12 BP神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)Fig.12 BP neural network structure

步驟1：初始化BP 網(wǎng)絡(luò)參數(shù)，包括輸入層、隱含層和輸出層神經(jīng)元個數(shù)，初始化權(quán)值w和閾值b，設(shè)置最大訓(xùn)練次數(shù)、最小精度和學(xué)習(xí)率。

步驟2：正向過程計算，通過輸入試驗樣本值，計算各層輸出和訓(xùn)練損失函數(shù)，如果滿足終止訓(xùn)練的要求，則近似模型訓(xùn)練完成，此時神經(jīng)網(wǎng)絡(luò)達到預(yù)期目標。

步驟3：誤差計算，判斷是否滿足收斂條件。

步驟4：反向過程計算，如果不滿足終止條件，將誤差反向傳輸并更新權(quán)值和偏置，繼續(xù)訓(xùn)練。

3.3 近似模型程序設(shè)計

如圖13 所示，本文通過SLE 試驗設(shè)計方法和BP 神經(jīng)網(wǎng)絡(luò)相結(jié)合，建立的近似模型可代替Dawson方法預(yù)報船體興波阻力。

圖13 近似模型流程圖Fig.13 Flow chart of approximate model

4 KCS船型計算分析

4.1 船型參數(shù)化設(shè)計

采用ITM 和FFD 方法對KCS 前半部分變形，設(shè)計變量如圖14 所示，α1和α2控制前半船在船長和船寬方向的變形，dx、dy、dz為控制框架里部分控制點沿著三個方向的位移，表1顯示了約束范圍均轉(zhuǎn)換為相對于船長的無因次量。采用SLE 方法生成400 個含有5 個特征變量的樣本船型，取Fr=0.26 并計算興波阻力值，從圖15中可以看出，較多樣本船在興波阻力方面表現(xiàn)優(yōu)于母型船。

圖14 KCS船體變形示意圖Fig.14 Diagram of KCS hull deformation

表1 變形方法參數(shù)設(shè)置Tab.1 Parameter setting of deformation method

圖15 船型樣本集Fig.15 Sample set of ship type

4.2 網(wǎng)絡(luò)參數(shù)對性能的影響分析

為了分析不同網(wǎng)絡(luò)參數(shù)對近似模型預(yù)測精度的影響，本文選擇均方誤差Emse和相關(guān)系數(shù)R2為評價指標：

式中，yi、y?i分別為樣本的理想值和預(yù)測值，為理想值的平均值，n為樣本數(shù)。

設(shè)網(wǎng)絡(luò)模型訓(xùn)練樣本數(shù)為400，最大訓(xùn)練次數(shù)為50 000次，訓(xùn)練精度為0.000 01。圖16顯示的是訓(xùn)練模型評價指標隨學(xué)習(xí)率的變化，可以看出相關(guān)系數(shù)越大，均方誤差越小。如圖17所示，隱含層的神經(jīng)元節(jié)點數(shù)影響近似模型的預(yù)測精度和訓(xùn)練產(chǎn)生的時間成本。根據(jù)數(shù)據(jù)對比結(jié)果可看出，本文學(xué)習(xí)率為0.99時，網(wǎng)絡(luò)模型的訓(xùn)練性能較好，當(dāng)隱含層節(jié)點數(shù)為6層時，相關(guān)系數(shù)較大，時間成本相對較低。

圖16 學(xué)習(xí)率對性能評價指標的影響Fig.16 Influence of learning rate on performance indicators

圖17 隱含層節(jié)點數(shù)對性能評價指標的影響Fig.17 Influence of number of hidden layer nodes on performance

4.3 訓(xùn)練實驗及結(jié)果分析

隱含層節(jié)點數(shù)和學(xué)習(xí)率參數(shù)確定之后，按照表2 網(wǎng)絡(luò)參數(shù)開始訓(xùn)練近似模型。為了驗證近似模型預(yù)測的準確性，在設(shè)計空間內(nèi)，隨機生成60 個測試樣本船型并計算興波阻力值。近似模型總誤差隨訓(xùn)練次數(shù)的變化曲線如圖18所示，可以看出總誤差E逐漸收斂為一個極小值。然后，使用訓(xùn)練之后的近似模型對測試樣本船型集的60組樣本進行預(yù)測，圖19為測試樣本預(yù)測值的相對誤差統(tǒng)計圖。

表2 近似模型結(jié)構(gòu)參數(shù)Tab.2 Structure parameters of approximate model

圖18 訓(xùn)練次數(shù)與總誤差的關(guān)系Fig.18 Relationship between training times and total error

圖19 測試樣本的相對誤差變化Fig.19 Relative error variation of the test sample

利用訓(xùn)練好的模型對60 組測試樣本Cws進行預(yù)測，平均相對誤差為0.282%，相關(guān)系數(shù)R2為0.9998，均方誤差為3.14E-04，運行時間為2 ms。其中，部分測試樣本的網(wǎng)絡(luò)輸出如表3 所示，表中輸入數(shù)據(jù)為歸一化后的值，通過分析表明，本文構(gòu)建的智能預(yù)報方法預(yù)測樣本船型的興波阻力的精度高且速度快。

表3 測試樣本的網(wǎng)絡(luò)輸出Tab.3 Network output of test samples

5 結(jié) 論

本文建立了基于樣本船的興波阻力智能預(yù)報方法，結(jié)合曲面變形技術(shù)和興波阻力預(yù)報技術(shù)，實現(xiàn)了在設(shè)計空間內(nèi)船體興波阻力的快速預(yù)報，并對近似模型的相關(guān)參數(shù)和訓(xùn)練展開分析，得出如下結(jié)論：

（1）提出了基于曲面控制網(wǎng)格的曲面變形方法，編寫了ITM 和FFD 方法的曲面變形程序，對KCS船型進行參數(shù)化設(shè)計并生成樣本船型集。

（2）在曲面變形技術(shù)和興波阻力預(yù)報技術(shù)基礎(chǔ)上，通過SLE試驗設(shè)計方法和BP神經(jīng)網(wǎng)絡(luò)相結(jié)合，建立的近似模型可快速預(yù)報船體興波阻力，耗時成本為毫秒量級。

（3）通過網(wǎng)絡(luò)模型參數(shù)分析和訓(xùn)練，測試樣本預(yù)測值相關(guān)系數(shù)為0.9998，均方誤差為3.14E-04，本文近似模型的預(yù)測精度較高。