李林斌，李洛東，劉日明

（中國船級社，北京 100007）

0 引 言

準(zhǔn)確預(yù)報船舶及海工結(jié)構(gòu)在水動力作用下的長期極值響應(yīng)是結(jié)構(gòu)設(shè)計的基礎(chǔ)，如船舶25年重現(xiàn)期的波浪彎矩或半潛式平臺100 年重現(xiàn)期的氣隙等。長期響應(yīng)預(yù)報需要考慮兩個因素，一是結(jié)構(gòu)在穩(wěn)定海況條件下響應(yīng)的短期分布，二是海況條件的長期分布，進而將每個海況條件對結(jié)構(gòu)長期極值響應(yīng)的貢獻率統(tǒng)計在內(nèi)。對于非線性響應(yīng)的情況（如波浪砰擊顫振誘導(dǎo)的船舶垂向彎矩），需要針對每個海況條件，采用非線性時域模擬和統(tǒng)計分析方法，來獲得結(jié)構(gòu)的響應(yīng)分布，因此非常耗費計算資源，不便于結(jié)構(gòu)的設(shè)計優(yōu)化和敏感性分析。

采用合理的替代方法，避免在整個海況域內(nèi)進行非線性響應(yīng)分析，是船舶和海洋工程領(lǐng)域持續(xù)研究的課題。其中，環(huán)境包絡(luò)線法在海洋工程領(lǐng)域有較為普遍的應(yīng)用[1]，這是一種基于短期極值來推算長期極值的方法，即在特定重現(xiàn)期環(huán)境包絡(luò)線設(shè)計點海況的短期最大響應(yīng)基礎(chǔ)上，乘以一系數(shù)（一般為1.1～1.3），或取短期最大響應(yīng)分布的更高分位值（一般為75%至90%分位），作為相同重現(xiàn)期的長期極值響應(yīng)，但該系數(shù)或分位值取決于結(jié)構(gòu)響應(yīng)的短期變異程度以及海況條件的長期變異程度，需要預(yù)先研究確定；另外一種方法是根據(jù)各海況對線性響應(yīng)長期極值的貢獻度，篩選出一定數(shù)量的海況，再針對這些特定海況進行更為復(fù)雜的非線性數(shù)值模擬和統(tǒng)計分析[2]，該方法成立的前提是海況的非線性和線性響應(yīng)貢獻度排序大致相當(dāng)，但這種假設(shè)并非普遍成立，同時并非所有的結(jié)構(gòu)線性響應(yīng)都存在且容易獲得。近年來，采用機器學(xué)習(xí)技術(shù)的回歸算法得到越來越多的重視和應(yīng)用，該方法一般采用平均取樣[3]或拉丁超立方體抽樣技術(shù)[4]，選取一定數(shù)量的海況進行水動力分析，獲得結(jié)構(gòu)響應(yīng)的分布參數(shù)，用以訓(xùn)練代理模型，并基于代理模型來預(yù)測其他海況對應(yīng)的響應(yīng)分布參數(shù)。由于需要保證代理模型在整個輸入域內(nèi)預(yù)測結(jié)果的準(zhǔn)確度，訓(xùn)練集樣本數(shù)往往很大，導(dǎo)致非線性響應(yīng)分析的計算量依然巨大?？傊?，這些替代性方法的應(yīng)用都存在一定的局限性。

本文提出一種新方法進行非線性水動力響應(yīng)的長期預(yù)報。該方法采用機器主動學(xué)習(xí)技術(shù)，通過機器學(xué)習(xí)模塊和水動力計算模塊的交互，以循環(huán)迭代的方式，貫序篩選對極值響應(yīng)貢獻最大的海況進行非線性時域分析，并將分析結(jié)果追加到訓(xùn)練集中，持續(xù)訓(xùn)練機器學(xué)習(xí)模型，進而實現(xiàn)預(yù)測值向真實值的快速收斂。

在后續(xù)的章節(jié)中，作者將簡要介紹采用傳統(tǒng)方法進行非線性水動力響應(yīng)長期預(yù)報的過程；闡述如何利用高斯過程回歸（GPR，Gaussian Process Regression）和貫序取樣（Sequential Sampling）技術(shù)來建立機器主動學(xué)習(xí)模型；并以兩型大型集裝箱船為算例，采用機器學(xué)習(xí)模型對彈性體砰擊顫振彎矩和剛體非線性波浪彎矩進行預(yù)測；通過與真實值的比較，驗證模型的高效性和準(zhǔn)確性，展示其在工程上的應(yīng)用價值；最后給出幾點研究結(jié)論和建議。

1 非線性水動力響應(yīng)的長期預(yù)報

1.1 長期極值響應(yīng)分析

船舶和海工結(jié)構(gòu)在水動力作用下的長期極值響應(yīng)，無論是線性的還是非線性的，都需要考慮結(jié)構(gòu)在穩(wěn)定海況條件下的響應(yīng)短期分布，以及海況條件的長期分布，即

式中：Q(y)為響應(yīng)Y的長期超越概率，目標(biāo)超越概率取決于重現(xiàn)期的要求；f(x)為海況x的長期分布，如波高-周期(Hs-Tz)聯(lián)合概率分布函數(shù)或波浪散布圖；g(y|x)為響應(yīng)Y在海況x下的短期分布（超越概率）。其中，響應(yīng)的短期分布g(y|x)是整個過程的核心。對于線性響應(yīng)，問題相對簡單，可以通過傳遞函數(shù)和譜分析的方法得到；對于非線性響應(yīng)，就需要針對每個海況，通過非線性時域水動力計算和統(tǒng)計分析獲得，需要在整個海況域內(nèi)求解，因此是一個非常耗時的過程。下面以船舶波浪砰擊顫振為例進一步說明。

1.2 船舶砰擊顫振彎矩的長期預(yù)報

大型船舶，特別是大型集裝箱船，由于首尾部船體大幅外飄，當(dāng)以較高航速在中高海況航行時，大幅的船波相對運動會使船體受到波浪的強烈沖擊，船體結(jié)構(gòu)由于砰擊載荷作用而發(fā)生瞬時高頻振動，產(chǎn)生所謂的“顫振”現(xiàn)象，具有強非線性特征，對總體波浪彎矩影響顯著，需要在船體梁極限強度分析時充分考慮。根據(jù)中國船級社《船體結(jié)構(gòu)波激振動和砰擊顫振直接計算評估指南》，應(yīng)采用時域水彈計算和統(tǒng)計分析方法進行預(yù)報，包括彈性體砰擊顫振誘導(dǎo)彎矩和剛體非線性波浪誘導(dǎo)彎矩。

對無限航區(qū)船舶，一般采用IACS REC No.34 推薦的北大西洋海浪散布圖[5]，其海況數(shù)為197 個（即：197個非零Hs-Tz參數(shù)對）。浪向角Hd的選取原則是在0°至360°之間，取間隔不大于30°，各浪向角的發(fā)生概率相等?？紤]到對稱性的因素，浪向角合計為7 個，如圖1 所示。這樣在考慮浪向后的海況總數(shù)為1379個，對如此多的海況進行逐一的非線性時域分析，對計算資源的耗費是可想而知的。

圖1 浪向角Hd及編號（括號內(nèi)的數(shù)字為Hd編號）Fig.1 Wave heading Hd and its numbering(the digit in brackets is the Hd numbering)

在短期海況xji（第j海況第i浪向角）條件下進行非線性時域計算，對得到的垂向波浪彎矩峰值y進行Weibull分布擬合，其超越概率為

式中，βji為形狀參數(shù)，ηji為尺度參數(shù)。

垂向波浪彎矩長期極值超越概率為

式中，pj為第j海況的長期概率，pi為第i浪向的長期概率。Q( )y= 10-8所對應(yīng)的響應(yīng)值yex為目標(biāo)極值響應(yīng)，即25年重現(xiàn)期的垂向波浪彎矩，本文中即為彈性體的砰擊顫振彎矩或剛體的非線性波浪彎矩。

海況xji對長期極值響應(yīng)的貢獻度Cji是一個很重要的概念，在后續(xù)的分析中會反復(fù)用到，其定義為

圖2 顯示了分析的總體過程。對Cji的進一步分析，不難發(fā)現(xiàn)長期極值響應(yīng)來自于極少數(shù)海況的貢獻。以13 500 TEU顫振中垂彎矩為例（具體見3.1節(jié)），貢獻度前30的海況，其累計貢獻度為98.3%，而這些海況僅占海況總數(shù)（1379）的0.008%，均為有義波高Hs超過12.5 m 的高海況，浪向都是迎浪或艏斜浪（1號和2號浪向角），總體情況如圖3所示。

圖2 時域非線性波浪載荷計算和統(tǒng)計分析流程Fig.2 Time domain nonlinear wave load calculation and statistical analysis process

圖3 海況對長期極值響應(yīng)的貢獻度（圖中圓的直徑與貢獻度成正比，顏色對應(yīng)浪向角）Fig.3 Contributions of sea states to the long-term extreme response(the diameter of the circles is proportional to the contributions,and the colors correspond to the wave headings)

如果能夠建立一種模型，以智能和貫序方式識別出這些海況，并傳遞給水動力分析模塊進行非線性時域分析，而不是分析全部海況，將大大減少計算資源的耗費，提高結(jié)構(gòu)設(shè)計和分析效率，這正是本文研究的目的和出發(fā)點。

2 機器主動學(xué)習(xí)模型

2.1 高斯過程回歸

高斯過程回歸（GPR）是近年發(fā)展起來的一種機器學(xué)習(xí)回歸方法[6]，有著嚴(yán)格的統(tǒng)計學(xué)習(xí)理論基礎(chǔ)，對處理高維數(shù)、小樣本、非線性等復(fù)雜的問題具有很好的適應(yīng)性，且泛化能力強，有著廣泛的工程應(yīng)用[3，7]。GPR從函數(shù)空間角度出發(fā)，定義一個高斯過程來描述函數(shù)分布，直接在函數(shù)空間進行貝葉斯推理。GPR 是任意有限個隨機變量均具有聯(lián)合高斯分布的集合，其性質(zhì)完全由均值函數(shù)和協(xié)方差函數(shù)確定，即f(x)～GP(μ(x),k(x,x') )，其中

式中，x,x'∈Rd為任意隨機變量。通常會對數(shù)據(jù)作預(yù)處理，使其均值函數(shù)等于0。假設(shè)觀測過程受到加性噪聲的污染，則觀測值y可用如下模型描述：

式中，x為輸入向量，f為函數(shù)值，ε為加性噪聲污染。假設(shè)ε符合均值為0、方差為c2的正態(tài)分布，即ε～N(0,c2)，則觀測值y的先驗分布為

觀測值y和預(yù)測值f*的聯(lián)合先驗分布為

式中：K(X,X)=Kn=(kij)為n×n階協(xié)方差矩陣，kij=k(xi,xj)，為xi與xj的協(xié)方差；K(X,x*)=K(x*,X)T為測試點x*與訓(xùn)練集的輸入X之間的n× 1 階協(xié)方差矩陣；k()x*,x*為測試點x*的方差；In為n維單位矩陣。

由此可以得到預(yù)測值f*的后驗分布為

為了訓(xùn)練GPR模型，需要設(shè)定一個協(xié)方差函數(shù)，文中我們采用了高斯徑向基核函數(shù)[3]：

式中，δ2為核函數(shù)的方差參數(shù)，h為尺度參數(shù)。δ和h這兩個超參數(shù)以及c需要通過機器學(xué)習(xí)（即訓(xùn)練）獲得。

2.2 貫續(xù)取樣算法

GPR 是監(jiān)督學(xué)習(xí)的一種，是根據(jù)經(jīng)驗數(shù)據(jù)（訓(xùn)練集）來學(xué)習(xí)輸入-輸出之間的映射關(guān)系，使得給定新的輸入便可得到相應(yīng)的輸出值（預(yù)測值）。通常情況下，訓(xùn)練集需要人為預(yù)先指定。相比人工神經(jīng)網(wǎng)絡(luò)、隨機森林、支持向量機等其他機器學(xué)習(xí)方法，高斯過程回歸的獨到之處在于其預(yù)測值是一個高斯隨機變量，不但可以給出預(yù)測值的最佳估計（均值），而且可以給出變異程度（方差）。

正是GPR以上獨有的特點，為實現(xiàn)機器主動學(xué)習(xí)提供了可能，并應(yīng)用于非線性動力系統(tǒng)的長期極值響應(yīng)預(yù)報[8]，即貫序找出對長期響應(yīng)影響最大的輸入值，將其加入到訓(xùn)練集，迭代更新GPR 模型，使其預(yù)測值在對長期響應(yīng)貢獻度最大的輸入?yún)^(qū)域內(nèi)變得更為準(zhǔn)確，進而實現(xiàn)預(yù)測值向真實值的快速收斂。結(jié)合本文的研究目的，我們設(shè)計了一種算法，通過Python編程，開發(fā)了機器主動學(xué)習(xí)模塊，并通過與水動力分析模塊的交互，實現(xiàn)了對船舶結(jié)構(gòu)非線性水動力響應(yīng)的長期預(yù)報，具體步驟如下：

（1）指定初始訓(xùn)練集的輸入集，即海況x=()Hs,Tz,Hd，進行非線性時域水動力計算和統(tǒng)計分析，得到訓(xùn)練集的輸出集，即Weibull分布的形狀參數(shù)β和位置參數(shù)η；

（2）對GPR模型進行訓(xùn)練，分別獲得超參數(shù)δβ、hβ和cβ以及δη、hη和cη；

（3）利用GPR模型預(yù)測其他海況xji對應(yīng)的觀測值均值和，以及方差和;

Murphy指出，個體在實際展開工作活動的過程中，為了達到工作績效目標(biāo)所表現(xiàn)出的一系列行為，即為績效；我國在積極展開工作績效研究的過程中，指出了學(xué)習(xí)過程、創(chuàng)新行為、公民氣候、技術(shù)核心即學(xué)習(xí)績效、創(chuàng)新績效、關(guān)系績效和任務(wù)績效四個工作績效結(jié)構(gòu)[2]。

（5）計算每一個海況的“收獲函數(shù)”Aji=，找出Aji最大時對應(yīng)的海況xn；

（6）將xn傳遞給水動力分析模塊，進行非線性時域計算和統(tǒng)計分析，獲得Weibull 分布的形狀參數(shù)βn和位置參數(shù)ηn；

（7）將xn，βn和ηn追加到訓(xùn)練集中，重復(fù)（2）～（6），直到長期響應(yīng)值收斂。

以上過程中，“下一個最優(yōu)”海況的選擇應(yīng)考慮和平衡兩個因素，一是需要距離訓(xùn)練集足夠遠（也即預(yù)測值方差足夠大），二是對長期極值響應(yīng)y^ex的貢獻度足夠大，而“收獲函數(shù)”實現(xiàn)了這一平衡。

3 算例研究

3.1 誤差分析

本文以13 500 TEU 和20 000 TEU 兩型大型集裝箱船為算例進行研究，建立水動力模型并進行時域水彈計算和統(tǒng)計分析，包括彈性體砰擊顫振彎矩和剛體非線性波浪彎矩。這兩型集裝箱船的主要參數(shù)如表1所列。

表1 船舶主要參數(shù)Tab.1 Main parameters of the ships

首先采用傳統(tǒng)方法，按照1.2 節(jié)所述流程，完成了全部1379 個海況的水動力計算和統(tǒng)計分析，并獲得了長期極值響應(yīng)yex，并將其作為“真實值”。

然后采用第2 章所述的機器學(xué)習(xí)模型進行長期極值響應(yīng)預(yù)報。選擇波浪散布圖中最小和最大Hs，及其對應(yīng)Tz出現(xiàn)概率最大的海況，對這些海況進行水動力計算和統(tǒng)計分析，獲得對應(yīng)的分布參數(shù)，作為初始訓(xùn)練集。

式（13）為13 500 TEU船型顫振中垂彎矩GPR模型的初始訓(xùn)練集：

我們采用機器學(xué)習(xí)模型迭代100次，每次迭代都獲得一次長期極值響應(yīng)的預(yù)報值y^ex；同時模型將根據(jù)“收獲函數(shù)”，自動和貫序識別出“下一個”最重要的海況，傳遞給水動力模塊進行時域分析。給出的誤差e如式（14）所示，以反映機器學(xué)習(xí)模型的準(zhǔn)確性：

誤差分析的結(jié)果如圖4～7 所示。不難看出，機器學(xué)習(xí)模型的預(yù)測結(jié)果誤差隨著迭代次數(shù)的增加而快速降低和收斂。雖然由于結(jié)構(gòu)響應(yīng)的非線性、統(tǒng)計過程的不確定性，以及機器學(xué)習(xí)過程的隨機性等因素導(dǎo)致不同船型的收斂效率會有所差別，但從本文的算例看，模型迭代大約30～40 次后，其誤差均可以控制在1.5%以內(nèi)。加上初始訓(xùn)練集的14個海況，需要進行非線性水動力分析的海況總數(shù)不超過60個，相比傳統(tǒng)方法需要分析的全部1379個海況，計算效率提高了20多倍。

圖4 13 500 TEU彈性體砰擊顫振彎矩誤差Fig.4 Errors of the 13 500 TEU elastic hull whipping bending moments

圖5 13 500 TEU剛體非線性波浪彎矩誤差Fig.5 Errors of the 13 500 TEU rigid body nonlinear wave bending moments

圖6 20 000 TEU彈性體砰擊顫振彎矩誤差Fig.6 Errors of the 20 000 TEU elastic hull whipping bending moments

圖7 20 000 TEU剛體非線性波浪彎矩誤差Fig.7 Errors of the 20 000 TEU rigid body nonlinear wave bending moments

圖8 顯示了機器學(xué)習(xí)模型貫序取樣得到的前20 個海況的排序（13 500 TEU 彈性體顫振中拱彎矩）。與圖3 對比可以看出，機器學(xué)習(xí)模型很好地識別出了對長期極值響應(yīng)貢獻度最大的海況區(qū)域，顯示出機器學(xué)習(xí)技術(shù)的強大功能。

圖8 貫序取樣前20個海況的分布，0表示初始訓(xùn)練集（13 500 TEU顫振中拱彎矩）Fig.8 Scatter of the top 20 sea states by sequential sampling and 0 means the initial training dataset(the 13 500 TEU whipping hogging moment)

3.2 敏感性分析

第2章所述的機器學(xué)習(xí)模型中，初始訓(xùn)練集的輸入集（海況）需要人為預(yù)先設(shè)定，并對收斂效率產(chǎn)生影響。海況涉及到三個參數(shù)，即：有義波高Hs、跨零周期Tz和浪向角Hd。我們以13 500 TEU 顫振中垂彎矩為例，通過分別改變式（13）中的這三個海況參數(shù)，進行敏感性分析，結(jié)果匯總?cè)鐖D9所示。

圖9 訓(xùn)練集的敏感性分析（13 500 TEU彈性體顫振中垂彎矩）Fig.9 Sensitivity analysis of training data(the 13 500 TEU elastic hull whipping sagging moments)

（1）改變有義波高

一般情況下，船舶結(jié)構(gòu)的水動力長期極值響應(yīng)主要來自于大波高海況（如Hs超過10 m 的海況），圖3 也清晰反映了這一趨勢。因此，在確定初始訓(xùn)練集的下限海況時，似乎應(yīng)選擇更高的海況，以便提高收斂效率和計算精度，但事實并非如此。如我們選擇(Hs,Tz)=(9.5,10.5)作為下限海況，上限海況保持不變，可以看到模型的收斂效率明顯下降，原因在于GPR 是一種適于內(nèi)插預(yù)測的機器學(xué)習(xí)方法，而外推預(yù)測時誤差往往會很大，即模型對Hs低于9.5 m海況的貢獻度預(yù)測產(chǎn)生很大誤差，進而導(dǎo)致收斂效率降低。

（2）改變波浪周期

選擇波浪散布圖中的最小波高和最小周期(Hs,Tz)=(0.5,3.5)，以及最大波高和最大周期(Hs,Tz)=(16.5,15.5)，作為初始訓(xùn)練集海況，其收斂效率會有所下降，但不是很明顯。

（3）改變浪向角

只保留1/4/7三個浪向，將式（13）的海況訓(xùn)練集樣本數(shù)縮減至6個。可以看出，誤差收斂效率有所降低，需要迭代更多次才能達到預(yù)期的誤差范圍。

根據(jù)以上分析，在確定初始訓(xùn)練集海況時，我們應(yīng)選擇波浪散布圖中最小和最大Hs及其對應(yīng)Tz出現(xiàn)概率最大的海況，以及所有的浪向角，以保證迭代效率和預(yù)測精度。

4 結(jié) 語

船舶和海工結(jié)構(gòu)在水動力作用下的響應(yīng)往往具有非線性特征，對長期響應(yīng)極值的準(zhǔn)確預(yù)報通常需要在整個海況域上進行時域非線性數(shù)值模擬和統(tǒng)計分析，計算資源耗費巨大。本文采用高斯過程回歸和貫序取樣技術(shù)，設(shè)計了機器主動學(xué)習(xí)的算法，并通過編程開發(fā)出了相應(yīng)的應(yīng)用模型，實現(xiàn)了基于少數(shù)關(guān)鍵海況的分析結(jié)果來準(zhǔn)確預(yù)測長期響應(yīng)極值的目標(biāo)。

本文以兩型大型集裝箱船為算例，采用機器學(xué)習(xí)模型，對彈性體砰擊顫振彎矩和剛體非線性波浪彎矩的長期響應(yīng)極值進行了預(yù)報，結(jié)果表明其預(yù)測值隨著迭代次數(shù)增加而快速收斂，在可接受的誤差范圍內(nèi)，較傳統(tǒng)方法提高計算效率20多倍，達到了工程化應(yīng)用的程度。對貫序取樣結(jié)果的分析表明，模型可以很好地識別出對極值響應(yīng)影響最大的輸入?yún)^(qū)域，顯示了機器學(xué)習(xí)技術(shù)的強大功能。

以上機器學(xué)習(xí)模型中沒有引入特定假設(shè)，因此具有廣泛的適用性，做稍許適應(yīng)性改造就可以用于其他非線性響應(yīng)的長期預(yù)報。如果將這一機器學(xué)習(xí)模塊嵌入到水動力分析軟件中，實現(xiàn)模型迭代結(jié)果和水動力計算結(jié)果的自動交互，將大幅提高軟件的智能化水平和計算效率。