基于機(jī)器學(xué)習(xí)的脈動源格林函數(shù)預(yù)報初探

常宏宇，朱仁傳，黃 山

(上海交通大學(xué) 船舶海洋與建筑工程學(xué)院 海洋工程國家重點(diǎn)實驗室，上海 200240)

格林函數(shù)法是計算波浪與結(jié)構(gòu)物相互作用問題的一個重要方法，快速準(zhǔn)確地計算格林函數(shù)及其偏導(dǎo)數(shù)是準(zhǔn)確獲得浮體作用力的關(guān)鍵和難點(diǎn)。三維勢流水動力理論中格林函數(shù)法分為簡單格林函數(shù)法和自由面格林函數(shù)法(亦稱為復(fù)雜格林函數(shù)法)。數(shù)值計算時前者需要的離散自由面，流場截斷，遠(yuǎn)場輻射邊界的處理都是需要克服的難點(diǎn)。后者則特指采用無限流域中滿足線性自由面條件的波動格林函數(shù)，數(shù)值實現(xiàn)時，不需要離散自由面，網(wǎng)格量少，其困難在于復(fù)雜的格林函數(shù)本身的計算[1]。為此以無因次表達(dá)的脈動源格林函數(shù)為例，進(jìn)行機(jī)器學(xué)習(xí)研究，以期獲得高效高精度的計算方法。

精度和效率是格林函數(shù)法有效實施的兩個基本要素。頻域零航速自由面格林函數(shù)，按照表達(dá)形式的不同，可分為Havelock型、Haskind-Havelock型以及Kim型格林函數(shù)，其數(shù)值計算有偏微分方程法、數(shù)值積分法、級數(shù)展開法以及多項式逼近法。解析逼近計算高效但有時難以控制精度，為了避免精度損失往往需要對區(qū)域進(jìn)行有效劃分。Noblesse[2-3]采用指數(shù)積分形式提出了兩種互補(bǔ)的近場和遠(yuǎn)場單積分表示法，用漸近展開級數(shù)和收斂升序級數(shù)分別計算源點(diǎn)和場點(diǎn)之間不同距離的格林函數(shù)，當(dāng)無量綱坐標(biāo)趨近零時，給出了兩個互補(bǔ)的泰勒級數(shù)展開式。Telste和Noblesse[4]基于文獻(xiàn)[3]的工作，將計算域劃分為5個子域，發(fā)表了一個數(shù)值計算代碼來評估格林函數(shù)及偏導(dǎo)數(shù)。Newman[5]將速度勢表示為Bessel函數(shù)、Struve函數(shù)與一個有限積分之和，使用龍貝格積分求解得到雙精度的格林函數(shù)值。之后他還提出了一種將解析展開式和切比雪夫多項式逼近結(jié)合的可達(dá)10-6精度的級數(shù)展開法[6]和一種全域切比雪夫多項式逼近法[7]，其中文獻(xiàn)[6]中的級數(shù)展開方法是商業(yè)軟件WAMIT的計算核心基礎(chǔ)。Chen[8-9]等提出了一種基于多項式級數(shù)逼近的算法，采用雙重切比雪夫逼近和特殊函數(shù)法對無限水深和有限水深下的格林函數(shù)進(jìn)行計算，算法應(yīng)用在水動力計算軟件HYDROSTAR中[10]。Clement[11]基于格林函數(shù)的二階微分方程提出了一種完全不同的方法，其難點(diǎn)在于初始化。近年來，Shen等[12]提出了一種二階常微分方程算法，算法達(dá)到10-6精度，相比級數(shù)展開法更高效。Duan等[13]和Shan等[14]分別給出了一種基于全域劃分的級數(shù)展開法，全域可以達(dá)到10-9精度。Wu等[15]提出了一種全域多項式擬合方法，將格林函數(shù)及其梯度表示為基本空間奇點(diǎn)、非振蕩的近場項和波動項，針對近場項給出簡單的近似計算式，只涉及實變量和基本初等函數(shù)(代數(shù)、指數(shù)、對數(shù))，方法可達(dá)10-3精度，且極大地提高了計算效率。Shan等[16]基于Newman[7]的工作，優(yōu)化了剩余函數(shù)的選擇區(qū)域，對子域重新劃分，改進(jìn)了切比雪夫逼近法。

機(jī)器學(xué)習(xí)是人工智能的核心，包括許多可以通過經(jīng)驗和數(shù)據(jù)來自動改進(jìn)的計算機(jī)算法，如支持向量機(jī)和神經(jīng)網(wǎng)絡(luò)等。當(dāng)前機(jī)器學(xué)習(xí)在海洋工程領(lǐng)域上也有所應(yīng)用：董寶玉[17]研究了一種基于遺傳算法優(yōu)化支持向量機(jī)多分類的方法，擴(kuò)展到高維空間應(yīng)用于船舶柴油機(jī)故障診斷中，成功對故障模型進(jìn)行預(yù)判；遲岑[18]利用支持向量機(jī)進(jìn)行無人艇自主避障中靜止障礙物的聚類，利用ELAMN預(yù)測模型中動態(tài)障礙物的運(yùn)動狀態(tài)，順利完成避障。深度學(xué)習(xí)是機(jī)器學(xué)習(xí)的第二次浪潮，Hinton等[19]曾指出，多隱層神經(jīng)網(wǎng)絡(luò)具有優(yōu)異的特征學(xué)習(xí)能力，其實質(zhì)是通過構(gòu)建多隱層的神經(jīng)元模型，使用海量訓(xùn)練數(shù)據(jù)來學(xué)習(xí)有代表性的特征，最終提升分類或預(yù)測的準(zhǔn)確性。

前人對復(fù)雜格林函數(shù)的計算多是數(shù)值積分或以解析函數(shù)為主的多項式逼近，前者精度高但計算耗時，后者計算高效但精度難以控制。三維頻域深水繞輻射問題的脈動點(diǎn)源格林函數(shù)的數(shù)值計算分析研究較多，有直接積分計算、多項式逼近等，為便于研究和比較，以脈動點(diǎn)源格林函數(shù)為例，先參考Newman方法[5]，使用龍貝格積分得到高精度的格林函數(shù)數(shù)據(jù)庫。引入機(jī)器學(xué)習(xí)方法，采用深度學(xué)習(xí)函數(shù)庫Keras對數(shù)據(jù)庫進(jìn)行學(xué)習(xí)，通過反向傳播算法迭代更新模型參數(shù)，并調(diào)節(jié)模型深度和激活函數(shù)等超參數(shù)，建立格林函數(shù)的神經(jīng)網(wǎng)絡(luò)預(yù)報模型。分別對全局和局部的預(yù)報精度和效率進(jìn)行研究，并與相應(yīng)的數(shù)值積分和多項式逼近方法[15]比較，驗證了該方法的可行性和高效性。研究表明，基于神經(jīng)網(wǎng)絡(luò)的格林函數(shù)高精度預(yù)報模型，能夠保證較高的精度和效率，為水動力問題求解和數(shù)值計算難題提供了新的思路。

1 脈動點(diǎn)源格林函數(shù)解及其無因次表達(dá)

無限水深頻域脈動點(diǎn)源格林函數(shù)多見于教材，見文獻(xiàn)[20]。若場點(diǎn)P坐標(biāo)為(x,y,z)，源點(diǎn)Q坐標(biāo)為(ξ,η,ζ)，波數(shù)k0=ω2/g，脈動點(diǎn)源格林函數(shù)的表達(dá)式為：

(1)

(2)

將式(2)中的廣義波數(shù)m對k0做無因次化處理，得G*關(guān)于變量θ的單重積分表達(dá)式為：

(3)

記Z=z+ζ,X=k0R,Y=-k0Z=-k0(z+ζ),α=-Y+iXcosθ。由于(z+ζ)≤0，故X≥0,Y≥0，交換積分次序，格林函數(shù)歸納為：

(4)

進(jìn)一步可將格林函數(shù)表示成下面的θ型單重積分形式：

(5)

同時將式(2)按照Newman的方法[5]，將波動部分G*表達(dá)為以下無因次化形式：

(6)

其中，Y0是第二類貝塞爾函數(shù)，H0是Struve函數(shù)。采用龍貝格積分計算式(6)，得到10-15精度的格林函數(shù)G*的數(shù)據(jù)庫。

2 機(jī)器學(xué)習(xí)方法介紹

使用神經(jīng)網(wǎng)絡(luò)算法對模型進(jìn)行學(xué)習(xí)。神經(jīng)網(wǎng)絡(luò)是一種經(jīng)典的機(jī)器學(xué)習(xí)算法，其中的誤差逆?zhèn)鞑?BP)神經(jīng)網(wǎng)絡(luò)是一種監(jiān)督式學(xué)習(xí)算法，具有很強(qiáng)的非線性映射能力。Hornik[21]證明，只需一個包含足夠多神經(jīng)元的隱層，多層前饋網(wǎng)絡(luò)就能以任意精度逼近任意復(fù)雜度的連續(xù)函數(shù)。神經(jīng)網(wǎng)絡(luò)中最基本的成分是神經(jīng)元模型，如圖1所示。

在每一層的連接中，輸入層向量用x表示，經(jīng)激活函數(shù)f處理，輸出為f(W1x-θ)，激活函數(shù)需使用非線性函數(shù)，如sigmoid、tanh、relu等。給定訓(xùn)練數(shù)據(jù)，神經(jīng)網(wǎng)絡(luò)的所有參數(shù)即各層之間的連接權(quán)重ωi與偏置項θ，其求解是一個最優(yōu)化問題。

損失函數(shù)求解的優(yōu)化問題中，常沿著梯度方向不斷下降來找到極小值。常用的梯度下降法有批梯度下降、隨機(jī)梯度下降以及小批量梯度下降法。訓(xùn)練過程中一個epoch表示使用訓(xùn)練集中的全部樣本訓(xùn)練一次，batchsize表示批大小，即每次訓(xùn)練使用的訓(xùn)練集樣本個數(shù)。梯度下降法在梯度平緩的維度下降非常緩慢，在梯度陡峭的維度易抖動，陷入局部極小值或鞍點(diǎn)。而且選取合適的學(xué)習(xí)率比較困難，只有在原問題是凸問題的情況下，才能保證以任意精度取得最優(yōu)解。非凸情況下，需要對梯度下降進(jìn)行優(yōu)化，比如Momentum，Adagrad，Adam等算法，可以避免陷入極大值極小值，設(shè)置得當(dāng)可以得到全局最優(yōu)解。為了防止模型過擬合，文中還將探索正則化和隨機(jī)失活的影響。

圖1 神經(jīng)元模型Fig. 1 Neuron model

圖2 神經(jīng)網(wǎng)絡(luò)中的變量符號Fig. 2 Variable symbols in neural networks

Keras由python編寫，作為Tensorflow和Theano等的高階應(yīng)用程序接口，可以方便地進(jìn)行深度學(xué)習(xí)模型的調(diào)試、評估、應(yīng)用和可視化。關(guān)于格林函數(shù)的預(yù)報，可以直接使用Keras訓(xùn)練好的h5格式模型文件，也可以導(dǎo)出模型權(quán)重等參數(shù)后重構(gòu)模型以便靈活計算。

3 格林函數(shù)的預(yù)報驗證

3.1 格林函數(shù)數(shù)據(jù)庫樣本密度

使用龍貝格積分計算式(6)得到高精度的無因次格林函數(shù)數(shù)據(jù)庫。格林函數(shù)數(shù)值計算部分結(jié)果如表1所示。計算結(jié)果與文獻(xiàn)[5]結(jié)果完全一致，且應(yīng)用于求解無航速S175實船的水動力系數(shù)和運(yùn)動響應(yīng)的結(jié)果與WAMIT一致，表明了該格林函數(shù)數(shù)值計算結(jié)果的準(zhǔn)確性。

表1 數(shù)值計算結(jié)果Tab. 1 Numerical results

由于格林函數(shù)對任何頻率通用，應(yīng)用時只需生成關(guān)于無因次水平距離X和水深Y的統(tǒng)一數(shù)據(jù)庫。考慮船長、船寬和吃水，結(jié)合性能評估的無因次頻率范圍發(fā)現(xiàn)，X不小于20，Y不小于8時即可基本滿足浮體的性能評估要求。首先探究數(shù)據(jù)庫采樣的密度，X范圍為0～20，Y范圍為0～8，分別設(shè)置采樣密度為0.05×0.01、0.05×0.05和0.10×0.10的數(shù)據(jù)為訓(xùn)練集，采樣密度為0.01×0.01的數(shù)據(jù)作為驗證集，進(jìn)行精度分析。全局訓(xùn)練設(shè)置兩個隱層，各128個隱層節(jié)點(diǎn)，不加正則化和隨機(jī)失活。模型參數(shù)設(shè)置如表2所示。

表2 網(wǎng)絡(luò)結(jié)構(gòu)Tab. 2 Network structure

預(yù)報結(jié)果的誤差分析(對誤差絕對值取負(fù)對數(shù))如圖3所示。預(yù)報數(shù)據(jù)的誤差分布統(tǒng)計如圖4所示。誤差分布統(tǒng)計圖中，橫軸負(fù)對數(shù)坐標(biāo)誤差為5時，表示預(yù)報誤差在10-4～10-5之間。如圖4所示，模型1預(yù)報數(shù)據(jù)的誤差分布更加均勻，平均誤差更小，預(yù)報誤差有99.81%在10-3至10-6之間，預(yù)報誤差小于10-4的比例為66.21%，明顯高于模型2的9.34%和模型3的32.76%。3種訓(xùn)練集取樣密度均可保證10-4以上的平均誤差。將全局訓(xùn)練時格林函數(shù)數(shù)據(jù)庫的取樣密度設(shè)置為0.05×0.01，分區(qū)域訓(xùn)練時的取樣密度可以適當(dāng)加大，得到原始的無因次格林函數(shù)數(shù)據(jù)庫如圖5所示。

圖3 不同訓(xùn)練樣本密度下誤差分布Fig. 3 Error distribution for different training sample densities

圖4 誤差分布統(tǒng)計Fig. 4 Statistics of error distribution

圖5 全局格林函數(shù)庫Fig. 5 Global Green function library

3.2 全局訓(xùn)練優(yōu)化和精度分析

定義全連接神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu)，利用3.1節(jié)的格林函數(shù)數(shù)據(jù)庫，在Keras下進(jìn)行學(xué)習(xí)，考慮正則化、隨機(jī)失活及Adam等優(yōu)化算法。訓(xùn)練過程顯示，模型在10 000個epoch后收斂速度減緩，根據(jù)調(diào)參結(jié)果，設(shè)置如表3所示的8種網(wǎng)絡(luò)結(jié)構(gòu)。訓(xùn)練后導(dǎo)出模型，使用X范圍為0～20，Y范圍為0～8,間隔均為0.01的密集樣本點(diǎn)作為測試集，對預(yù)報精度進(jìn)行分析。

表3 全局訓(xùn)練網(wǎng)絡(luò)結(jié)構(gòu)Tab. 3 Network structure for global training

模型5和7中加入正則化以及Dropout后，模型的訓(xùn)練和測試誤差明顯偏大，表明在密集訓(xùn)練數(shù)據(jù)下，不需要考慮過擬合。其余各模型預(yù)報的格林函數(shù)結(jié)果如圖6所示。

圖6 全局預(yù)報結(jié)果Fig. 6 Global forecast result

預(yù)報結(jié)果的誤差分析(誤差絕對值取負(fù)對數(shù))如圖7所示。由圖7可知在區(qū)域邊緣附近，誤差擾動往往比較大，實際應(yīng)用過程中這一區(qū)域要做特殊處理。結(jié)果表明sigmoid激活函數(shù)和Adam優(yōu)化算法可以加快模型訓(xùn)練速度和預(yù)報精度。其他條件一定時，增加神經(jīng)網(wǎng)絡(luò)模型深度會增加訓(xùn)練時間，但擬合精度大致不變。模型1預(yù)報誤差有99.81%在10-4至10-6之間，預(yù)報誤差小于10-5的比例為66.21%，這個精度可以滿足水動力計算需求。模型1部分預(yù)報結(jié)果展示如表4所示。

圖7 全局預(yù)報誤差分布Fig. 7 Error distribution for global forecast

表4 預(yù)報結(jié)果展示Tab. 4 Presentation of forecast results

3.3 分區(qū)訓(xùn)練優(yōu)化和精度分析

如圖6所示，Y接近0時原始格林函數(shù)圖形波動較大。圖7表明這一區(qū)域的誤差較其他區(qū)域也偏大。為提高預(yù)報精度，考慮分區(qū)學(xué)習(xí)，將全局劃分為4區(qū)域：區(qū)域1中，X取0～3，Y取0～3；區(qū)域2中，X取3～20，Y取3～8；區(qū)域3中，X取0～3，Y取3～8；區(qū)域4中X取3～20，Y取0～3。4個區(qū)域分別進(jìn)行模型學(xué)習(xí)。

使用的分區(qū)訓(xùn)練數(shù)據(jù)完全互補(bǔ)，但并不重合，在預(yù)報時，區(qū)域交界點(diǎn)只屬于4個分區(qū)中的某一個區(qū)域。各區(qū)域的原始格林函數(shù)圖像與圖5一致。

3.3.1 區(qū)域1精度分析

區(qū)域1上，X和Y每隔0.01取樣，共90 000個訓(xùn)練數(shù)據(jù)。設(shè)置表5所示的4種網(wǎng)絡(luò)結(jié)構(gòu)。

表5 區(qū)域1網(wǎng)絡(luò)結(jié)構(gòu)Tab. 5 Network structure for zone 1

結(jié)果表明，模型2和4加入正則化后，預(yù)報誤差明顯增大，說明在訓(xùn)練樣本密集的情況下，模型不會過擬合，后續(xù)將不再考慮正則化。模型1和3預(yù)報結(jié)果的誤差分析(誤差絕對值取負(fù)對數(shù))如圖8所示。

圖8 區(qū)域1誤差分布Fig. 8 Error distribution of zone 1

平均誤差最小的模型1預(yù)報誤差有99.51%在10-3至10-6之間，誤差小于10-4的比例為94.56%。

3.3.2 區(qū)域2精度分析

在區(qū)域2上，X每隔0.02，Y每隔0.01取樣，共426 351個訓(xùn)練數(shù)據(jù)。本節(jié)及之后的分區(qū)訓(xùn)練都設(shè)置如表6所示的4種網(wǎng)絡(luò)結(jié)構(gòu)。

表6 區(qū)域2/3/4網(wǎng)絡(luò)結(jié)構(gòu)Tab. 6 Network structure for zones 2/3/4

預(yù)報結(jié)果的誤差分析(誤差絕對值取負(fù)對數(shù))如圖9所示。

圖9 區(qū)域2誤差分布Fig. 9 Error distribution of zone 2

平均誤差最小的模型1的預(yù)報誤差有96.53%在10-4至10-7之間，預(yù)報誤差小于10-5的比例為96.55%。

3.3.3 區(qū)域3精度分析

在區(qū)域3上，X和Y均每隔0.01取樣，共150 300個訓(xùn)練數(shù)據(jù)。預(yù)報結(jié)果的誤差分析(誤差絕對值取對數(shù))如圖10所示。

圖10 區(qū)域3誤差分布Fig. 10 Error distribution of zone 3

平均誤差最小的模型2的預(yù)報誤差有97.17%在10-4至10-7之間，預(yù)報誤差小于10-5的比例為99.80%。

3.3.4 區(qū)域4精度分析

在區(qū)域4上，X和Y每隔0.01取樣，共510 000個訓(xùn)練數(shù)據(jù)。預(yù)報結(jié)果的誤差分析(誤差絕對值取對數(shù))如圖11所示。

圖11 區(qū)域4誤差分布Fig. 11 Error distribution of zone 4

平均誤差最小的模型2的預(yù)報誤差有96.86%在10-4至10-7之間，預(yù)報誤差小于10-5的比例為21.94%。

結(jié)果表明，分區(qū)預(yù)報時，各區(qū)域的精度指標(biāo)基本都優(yōu)于全局預(yù)報，絕大部分區(qū)域的精度可達(dá)10-4以上。文獻(xiàn)[22]的研究表明對于水動力計算這個精度基本滿足。

分區(qū)訓(xùn)練時，分區(qū)訓(xùn)練數(shù)據(jù)完全互補(bǔ)，但是并不重合。在預(yù)報時，區(qū)域交界點(diǎn)只屬于4個分區(qū)中的一個區(qū)域，在后續(xù)的水動力系數(shù)和運(yùn)動響應(yīng)的計算中，替代的子程序只調(diào)用4個分區(qū)代碼中的某一個即可。

3.4 預(yù)報效率對比

為了評估機(jī)器學(xué)習(xí)方法的效率，將神經(jīng)網(wǎng)絡(luò)模型與數(shù)值積分方法及多項式逼近方法進(jìn)行效率對比，對比的數(shù)值積分方法有Newman方法[5]和θ型單重積分，多項式逼近方法為文獻(xiàn)[15]的全局逼近。

單次計算平均用時，指在不同的預(yù)報數(shù)據(jù)的規(guī)模下，全部完成預(yù)報的用時與數(shù)據(jù)規(guī)模的比值，表示多次計算的平均用時。python中的矩陣計算依賴于numpy函數(shù)庫，對向量和矩陣的計算進(jìn)行了優(yōu)化加速。神經(jīng)網(wǎng)絡(luò)模型在進(jìn)行較大規(guī)模數(shù)據(jù)的預(yù)報時，進(jìn)行前饋網(wǎng)絡(luò)的逐層計算，等價于多個矩陣計算。算法在預(yù)報時將全部輸入視為大規(guī)模矩陣進(jìn)行加速，在數(shù)據(jù)規(guī)模更大時，神經(jīng)網(wǎng)絡(luò)算法會表現(xiàn)出更好的性能，單次計算的平均用時更短。而這也是機(jī)器學(xué)習(xí)算法的優(yōu)勢，數(shù)據(jù)庫以及訓(xùn)練算法一旦形成，可以廣泛應(yīng)用于不同對象，而不需要每次都重新計算。預(yù)報效率的對比如表7所示。

表7 不同數(shù)據(jù)規(guī)模下單次計算平均用時 Tab. 7 Average time of single calculation for different data sizes (μs)

如表7所示，機(jī)器學(xué)習(xí)預(yù)報的計算效率低于以解析函數(shù)為主的多項式逼近，高于θ型單重積分?jǐn)?shù)值方法，稍低于Newman的龍貝格積分?jǐn)?shù)值方法。

4 結(jié) 語

采用深度學(xué)習(xí)函數(shù)庫Keras，對脈動點(diǎn)源格林函數(shù)數(shù)據(jù)庫進(jìn)行學(xué)習(xí)，建立了神經(jīng)網(wǎng)絡(luò)預(yù)報模型，并展開了模型預(yù)報精度和效率的分析研究。得到如下結(jié)論：

1) 機(jī)器學(xué)習(xí)格林函數(shù)庫建立的神經(jīng)網(wǎng)絡(luò)模型預(yù)報格林函數(shù)是可行的。機(jī)器學(xué)習(xí)模型預(yù)報的格林函數(shù)在99%以上的區(qū)域可以保證10-4及以上的精度，能夠達(dá)到實用精度的要求。

2) 全局?jǐn)M合中，區(qū)域邊緣處的預(yù)報誤差偏大。論文采用的分區(qū)訓(xùn)練思路，有效地提高了精度。在分區(qū)擬合中，區(qū)域邊緣處的誤差會大于區(qū)域內(nèi)部，后續(xù)可考慮擴(kuò)展訓(xùn)練區(qū)域來提高精度。

3) 機(jī)器學(xué)習(xí)模型預(yù)報的效率優(yōu)于直接數(shù)值積分計算，低于以解析函數(shù)為主的多項式逼近，稍低于Newman的龍貝格積分?jǐn)?shù)值方法。

以上工作主要為建立在非常成熟的脈動點(diǎn)源格林函數(shù)研究基礎(chǔ)上的機(jī)器學(xué)習(xí)初步探索研究，此格林函數(shù)的計算有Newman的龍貝格數(shù)值積分方法和相關(guān)的解析逼近方法。機(jī)器學(xué)習(xí)建立的預(yù)報模型在此格林函數(shù)的計算效率上并沒有優(yōu)勢，但比直接求解數(shù)值計算的效率高很多。由此可以想象，對于那些復(fù)雜格林函數(shù)，特別是不能用解析多項式逼近表達(dá)的，如有航速三維頻域輻射格林函數(shù)，以及難以計算且效率低下的復(fù)雜函數(shù)等，可以合理地應(yīng)用機(jī)器學(xué)習(xí)方法，進(jìn)行建模替代計算。本方法為提高水動力問題求解效率，解決傳統(tǒng)計算難題提供了新的思路。