基于機器學習的脈動源格林函數預報初探

常宏宇，朱仁傳，黃 山

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

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

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

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

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

1 脈動點源格林函數解及其無因次表達

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

(1)

(2)

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

(3)

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

(4)

進一步可將格林函數表示成下面的θ型單重積分形式：

(5)

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

(6)

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

2 機器學習方法介紹

使用神經網絡算法對模型進行學習。神經網絡是一種經典的機器學習算法，其中的誤差逆傳播(BP)神經網絡是一種監督式學習算法，具有很強的非線性映射能力。Hornik[21]證明，只需一個包含足夠多神經元的隱層，多層前饋網絡就能以任意精度逼近任意復雜度的連續函數。神經網絡中最基本的成分是神經元模型，如圖1所示。

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

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

圖1 神經元模型Fig. 1 Neuron model

圖2 神經網絡中的變量符號Fig. 2 Variable symbols in neural networks

Keras由python編寫，作為Tensorflow和Theano等的高階應用程序接口，可以方便地進行深度學習模型的調試、評估、應用和可視化。關于格林函數的預報，可以直接使用Keras訓練好的h5格式模型文件，也可以導出模型權重等參數后重構模型以便靈活計算。

3 格林函數的預報驗證

3.1 格林函數數據庫樣本密度

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

表1 數值計算結果Tab. 1 Numerical results

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

表2 網絡結構Tab. 2 Network structure

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

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

圖4 誤差分布統計Fig. 4 Statistics of error distribution

圖5 全局格林函數庫Fig. 5 Global Green function library

3.2 全局訓練優化和精度分析

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

表3 全局訓練網絡結構Tab. 3 Network structure for global training

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

圖6 全局預報結果Fig. 6 Global forecast result

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

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

表4 預報結果展示Tab. 4 Presentation of forecast results

3.3 分區訓練優化和精度分析

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

使用的分區訓練數據完全互補，但并不重合，在預報時，區域交界點只屬于4個分區中的某一個區域。各區域的原始格林函數圖像與圖5一致。

3.3.1 區域1精度分析

區域1上，X和Y每隔0.01取樣，共90 000個訓練數據。設置表5所示的4種網絡結構。

表5 區域1網絡結構Tab. 5 Network structure for zone 1

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

圖8 區域1誤差分布Fig. 8 Error distribution of zone 1

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

3.3.2 區域2精度分析

在區域2上，X每隔0.02，Y每隔0.01取樣，共426 351個訓練數據。本節及之后的分區訓練都設置如表6所示的4種網絡結構。

表6 區域2/3/4網絡結構Tab. 6 Network structure for zones 2/3/4

預報結果的誤差分析(誤差絕對值取負對數)如圖9所示。

圖9 區域2誤差分布Fig. 9 Error distribution of zone 2

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

3.3.3 區域3精度分析

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

圖10 區域3誤差分布Fig. 10 Error distribution of zone 3

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

3.3.4 區域4精度分析

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

圖11 區域4誤差分布Fig. 11 Error distribution of zone 4

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

結果表明，分區預報時，各區域的精度指標基本都優于全局預報，絕大部分區域的精度可達10-4以上。文獻[22]的研究表明對于水動力計算這個精度基本滿足。

分區訓練時，分區訓練數據完全互補，但是并不重合。在預報時，區域交界點只屬于4個分區中的一個區域，在后續的水動力系數和運動響應的計算中，替代的子程序只調用4個分區代碼中的某一個即可。

3.4 預報效率對比

為了評估機器學習方法的效率，將神經網絡模型與數值積分方法及多項式逼近方法進行效率對比，對比的數值積分方法有Newman方法[5]和θ型單重積分，多項式逼近方法為文獻[15]的全局逼近。

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

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

如表7所示，機器學習預報的計算效率低于以解析函數為主的多項式逼近，高于θ型單重積分數值方法，稍低于Newman的龍貝格積分數值方法。

4 結 語

采用深度學習函數庫Keras，對脈動點源格林函數數據庫進行學習，建立了神經網絡預報模型，并展開了模型預報精度和效率的分析研究。得到如下結論：

1) 機器學習格林函數庫建立的神經網絡模型預報格林函數是可行的。機器學習模型預報的格林函數在99%以上的區域可以保證10-4及以上的精度，能夠達到實用精度的要求。

2) 全局擬合中，區域邊緣處的預報誤差偏大。論文采用的分區訓練思路，有效地提高了精度。在分區擬合中，區域邊緣處的誤差會大于區域內部，后續可考慮擴展訓練區域來提高精度。

3) 機器學習模型預報的效率優于直接數值積分計算，低于以解析函數為主的多項式逼近，稍低于Newman的龍貝格積分數值方法。

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