海冰工況下海洋運輸系統(tǒng)的仿真建模分析

張 明

（福建船政交通職業(yè)學院，福建 福州 350007）

中國的渤海和北黃海海域地區(qū)是我國重要的經(jīng)濟開發(fā)區(qū)，其中渤海是我國的海上交通要道，是國內(nèi)與國外溝通的海上門戶之一，且目前已發(fā)現(xiàn)豐富的油氣等資源，油氣勘探、航海運輸?shù)认到y(tǒng)正在快速的發(fā)展［1］。但是，每年的冬季，這些海域都會發(fā)生海冰，尤其是遼東海灣，海冰尤為嚴重。這些海冰的堆積和漂移對我國冬季的油氣勘探尤其是海洋運輸系統(tǒng)都造成了巨大的威脅。歷史上曾經(jīng)發(fā)生過海冰堵塞航道、撞毀海上運輸船只以及推倒石油平臺的事故［2］，為保證海洋運輸系統(tǒng)的安全、適應(yīng)冬季海洋油氣的開發(fā)需求、減少并預(yù)防海冰可能會造成的災(zāi)害，亟需建立海冰工況綜合預(yù)報系統(tǒng)。

目前，我國對于海冰的監(jiān)測方式主要包括以下三種方式：一是日常的岸站監(jiān)測，即海冰的實測數(shù)據(jù)以及氣象數(shù)據(jù)進行預(yù)測。但是這種方式由于數(shù)據(jù)較為稀疏，往往以點窺面，誤差較大；二是利用海洋、氣象與海冰之間的統(tǒng)計學關(guān)系，按照以往的經(jīng)驗進行海冰的厚度、邊界等的預(yù)測；三是利用渤?；螯S海海域的熱力學和動力學研究，預(yù)測海冰的增長、消融和漂移趨勢［3］。這兩種方式受制于數(shù)值計算本身的精度和效率，且具有一定的時效性，無法進行長期的預(yù)報。

近年來出現(xiàn)的遙測法，包括采用衛(wèi)星遙感、無人機航測等方式搜集海冰相關(guān)資料，以及后來增加的船舶和野外觀察也開拓了海冰的資料收集，這些可以作為海冰量值定標的依據(jù)。將這些數(shù)據(jù)與日常海冰監(jiān)測方法結(jié)合，對海冰進行仿真建模分析，可以有效地為海洋運輸系統(tǒng)提供依據(jù)。因此，本文將在海冰工況下對海洋運輸系統(tǒng)的仿真建模進行分析。

1 海冰監(jiān)測系統(tǒng)設(shè)計

1.1 海冰資料的收集

對海冰工況進行仿真建模分析，其中的一個重要環(huán)節(jié)是收集海冰相關(guān)資料。而海冰的研究與多個學科相關(guān)，包括海洋環(huán)境、大氣等，所需要的資料包括衛(wèi)星遙感資料、雷達監(jiān)測資料、航測資料以及固定的平臺資料等。該海冰工況仿真建模過程的設(shè)計簡圖如圖1 所示。

圖1 海冰工況仿真建模過程的設(shè)計簡圖

在進行海冰資料收集時，一般通過中分辨率成像光譜儀（MODIS）拍攝衛(wèi)星遙感數(shù)據(jù)，以獲取大面積的海冰信息。對遙感數(shù)據(jù)進行處理時，其反演流程如圖2 所示。

圖2 遙感圖像的處理流程圖

對小尺度范圍的海冰圖像或者需要海冰的細節(jié)信息，則采用在海冰監(jiān)測站進行觀察，監(jiān)測信息主要包括會影響到海冰的氣象、水文和海冰信息。氣象信息包括風速、濕度、太陽輻射、氣壓和溫度，氣象信息每一秒測一次，可在監(jiān)測站實時顯示；水文信息包括水流速、流向、水文和鹽分等，每隔15分鐘采集一次，并在監(jiān)測站顯示；海冰信息包括海冰的類別、密集度、移動速度和移動軌跡等，這些采用雷達觀測，海冰厚度采用CCD 攝像機進行觀測。

2 海冰熱力學和動力學分析

在對海冰進行數(shù)值計算之前，首先需要對決定海冰凍結(jié)和融化的熱力學過程、決定海冰移動和變形的動力學過程、動力學和熱力學共同決定的海冰厚度變化、質(zhì)量和密集度的變化等進行研究［4-5］，才可對海冰進行后續(xù)的數(shù)值計算。

2.1 熱力學過程

大氣參數(shù)包括風速、相對濕度、氣溫、氣壓以及太陽輻射等，是海冰-海洋耦合的驅(qū)動條件?，F(xiàn)在研究的海冰熱力學模型都是在Stefan 公式的基礎(chǔ)上進行研究得到，該公式認為海面結(jié)冰后，海冰表面的溫度與氣溫近似，海冰底部的生長速度由結(jié)冰潛熱、冰內(nèi)部的熱傳導(dǎo)和冰內(nèi)部的溫度梯度決定，這種方法可以稱為冰凍融冰度日法。但是，不同海域的海冰狀況不同，但是基本原理相同，可以在Stefan 公式基礎(chǔ)上改動得到渤海區(qū)域的冰厚計算公式，渤海區(qū)域的平整冰厚度hs［6-7］為：

式中α 為海冰的增長系數(shù)；FD 為每日平均氣溫低于-2℃的累積的海冰冰凍日期；TD 為每日平均氣溫高于0℃的累積的海冰解凍日期；KD為海冰剛出現(xiàn)時所需要的平均海冰冰凍日期。不同海域的海冰增長系數(shù)以及所需要的平均冰凍日期不同，在進行其他海域的海冰厚度計算時，將這兩個系數(shù)更新即可。該公式忽略了太陽輻射以及海洋的熱通量等對于海冰產(chǎn)生熱影響的要素，因此在進行計算時，還需要將這些要素考慮進去。

在對海冰進行數(shù)值模擬仿真時，其熱力學模型的基礎(chǔ)是海冰內(nèi)部的熱力學平衡，這些熱力學因素主要為冰凍融冰度日法中沒有考慮到的因素，包括冰面上的太陽輻射、長波輻射、感熱和潛熱、海冰內(nèi)部的熱傳導(dǎo)以及底部的海洋熱能，海冰內(nèi)部的熱傳導(dǎo)可以通過下式計算［8］：

式中ps，cs和ks分別為海冰的平均密度、比熱容和熱傳導(dǎo)系數(shù)，Ts為海冰的平均溫度；λ 和β 為根據(jù)經(jīng)驗確定的常數(shù)；Ss（z）為海冰內(nèi)部距離海平面為z 深度的鹽度；QSe為太陽到海冰上表面的有效輻射量；γs為海冰對太陽輻射的衰減系數(shù)；τs為海冰對太陽光的透射率；以上各公式中的下標s 表示海冰。

海冰的厚度變化在海冰的上下表面同時進行。對于海冰上表面來說，若海冰表面的溫度大于海冰的融點，且表面熱量的吸收和散發(fā)的總和大于0，則海冰表面融化；對于海冰的下表面來說，其變化主要由海冰的內(nèi)部熱傳導(dǎo)、海洋熱通量以及太陽輻射透射量來決定。將海冰的上表面和下表面同時考慮，可以得到海冰厚度的增長率為：

式中，下標u和下標d分別表示海冰上表面和下表面厚度的增長率，則海冰上表面和下表面厚度的增長率分別為：

式中hs為海冰的總厚度，與平均厚度不同，也可叫做物理冰厚；Ls為海冰融化時的潛熱；Fw為海洋內(nèi)部的熱通量；Sd為到達海冰下表面的輻射透射量。

2.2 動力學過程

海冰在海洋上的漂移，其動力方程可通過牛頓定律表示，則單位面積海冰的動力學方程為：

式中M為海冰單位面積的質(zhì)量，按照下式表示：

式中Ps和A分別為海冰的密度和密集度。為海冰在海上運動時的速度矢量；f為科氏參數(shù)，通過下式計算：

式中的ω為地球旋轉(zhuǎn)的角速度；θ為海冰所處緯度位置。為與海面垂直的單位矢量；g為重力加速度；h0為海冰漂移處的海平面高度；為海冰內(nèi)應(yīng)力，當海冰密度較低或者與其他海冰無明顯的相互作用，則該項可忽略。和分別為海冰處的風和水流對海冰的拖拽力，計算方式為：

式中Pf和Pl分別為海洋上空風以及海水的密度；Cf和Cl分別為風和海水對海冰的拖拽系數(shù)。

2.3 連續(xù)方程

海冰的質(zhì)量變化是在熱力學與和動力學共同作將用下引起的，海冰的連續(xù)方程可以將海冰的熱力學和動力學過程耦合即可得到：

式中ΦA(chǔ)和Φh分別為海冰密集度、厚度的動力學變形系數(shù)；ΨA和Ψh分別為海冰密集度、厚度的熱力學系數(shù)；vx和vy分別為海冰運動時在x和y方向的速度分量。

2.4 本構(gòu)方程

在建立海冰的本構(gòu)方程時，將海冰看做二維的連續(xù)體，并采用Hibler 粘塑性本構(gòu)方程，該方程的建立是以符合橢圓屈服函數(shù)的二維粘塑性法則作為基礎(chǔ)。該橢圓屈服函數(shù)為：

式中σ0和σ1為海冰內(nèi)部的兩個主應(yīng)力；F0為海冰內(nèi)部強度；e 為橢圓的主軸比。Hibler 方程中，海冰的應(yīng)力和應(yīng)變之間的關(guān)系如下：

式中σmn為海冰的應(yīng)力張量；εmn為海冰的應(yīng)變速率張量；δmn為Kronecker 算子；ζ為非線性塊體的粘性系數(shù)；η為切向粘性系數(shù)。海冰內(nèi)部的靜水壓力F0為：

式中F＇和A均為根據(jù)經(jīng)驗確定的常數(shù)。

3 海冰數(shù)值模擬計算

3.1 海冰數(shù)值模擬計算方法的選擇

海冰在進行數(shù)值模擬計算時，一般將海冰看作是二維的流體或者流變體。針對流體的計算，按照采用坐標系的不同可以分為歐拉方法和拉格朗日方法。

歐拉方法適用于解決大變形的流體運動場。但對流場中包含兩種以上的流體介質(zhì)情況，介質(zhì)之間的分界面以及混合網(wǎng)格中的各類運輸量難以界定，使得人們不得不采取其他方法解決這類問題，海冰數(shù)值模擬即是這種情況。

拉格朗日方法適用于扭曲不嚴重的流體，尤其適用于多種介質(zhì)的整體或局部運動。但是，若流場變形較大，計算過程中容易出現(xiàn)網(wǎng)格畸形和交叉的情況，使計算中斷。對于海冰來說，其數(shù)值模型擾動較大且為多介質(zhì)，可以將拉格朗日和歐拉方法結(jié)合進行計算，即網(wǎng)格質(zhì)點的方法（PIC方法）。這種方法可以將網(wǎng)格的介質(zhì)劃分為若干質(zhì)點，這些質(zhì)點均有自己的質(zhì)量、動量和能量，通過對質(zhì)點進行運輸即可實現(xiàn)對流體介質(zhì)的運輸。拉格朗日坐標下的光滑質(zhì)點法（SPH）屬于網(wǎng)格質(zhì)點法（PIC 法），因此采用拉格朗日坐標下的光滑質(zhì)點法（SPH）對海冰數(shù)值進行模擬計算。

3.2 海冰數(shù)值的模擬計算

采用拉格朗日坐標系的光滑質(zhì)點法（SPH）對海冰進行數(shù)值模擬計算時，首先對海冰的動力方程進行處理可以得到：

其中海冰內(nèi)力這項可以寫為：

對于海冰中的第p個質(zhì)點，具有以下性質(zhì)：位置矢量rp，質(zhì)量mj，質(zhì)量密度Mp，海冰的平均厚度hp，質(zhì)點處的應(yīng)力張量σp，則海冰內(nèi)力可以表示為：

將上式帶入海冰的動力方程，并分別對動力方程的x和y向求解，可以得到p質(zhì)點的動量平衡方程如下式所示：

4 神經(jīng)網(wǎng)絡(luò)對海冰面積的預(yù)報研究

海冰的數(shù)值計算可以實現(xiàn)海冰面積的短期預(yù)測，盡管預(yù)報結(jié)果具有一定的精度的，但仍然存在一定的誤差。為能夠及時、有效地識別海冰并進行海冰預(yù)警，以盡早采取措施降低海冰造成的損失，本文引入了預(yù)報精度更高的神經(jīng)網(wǎng)絡(luò)模型對海冰面積進行預(yù)報，同時為保證海冰面積在檢測過程中的精度，采用卡爾曼濾波算法對預(yù)測結(jié)果進行跟蹤，以進一步提高預(yù)報精度。

4.1 人工神經(jīng)網(wǎng)絡(luò)模型的設(shè)計

人工神經(jīng)網(wǎng)絡(luò)是一種數(shù)學模型，通過利用類似大腦神經(jīng)突觸的結(jié)構(gòu)進行數(shù)據(jù)和信息的處理。在這些人工神經(jīng)網(wǎng)絡(luò)模型中，BP 神經(jīng)網(wǎng)絡(luò)或者其改進的形式是最常用的網(wǎng)絡(luò)模型。根據(jù)本文的需求，適合對海冰的面積進行預(yù)測的模型只有前饋網(wǎng)絡(luò)中的BP 神經(jīng)網(wǎng)絡(luò)。另一方面，BP 神經(jīng)網(wǎng)絡(luò)具有精度高、學習效率快、泛化能力強以及快速收斂的特點，因此本文采用BP 神經(jīng)網(wǎng)絡(luò)算法完成海冰面積的預(yù)測。

BP 神經(jīng)網(wǎng)絡(luò)模型在進行設(shè)計時，主要對其網(wǎng)絡(luò)層數(shù)、每層的神經(jīng)元個數(shù)、各層之間的激活函數(shù)、初始值和學習速率進行設(shè)計。對于BP 神經(jīng)網(wǎng)絡(luò)的輸入值層和輸出層，主要根據(jù)需要解決的問題以及最終的需求決定。根據(jù)設(shè)計要求，本文設(shè)計的BP 神經(jīng)網(wǎng)絡(luò)的輸入層主要包括海冰冰凍溫度、海冰融化溫度、每日平均氣溫和離岸大風引起的溫度變化，輸出層為未來四年的預(yù)測海冰面積。

隱含層節(jié)點數(shù)量對BP 神經(jīng)網(wǎng)絡(luò)的預(yù)測結(jié)果影響較大，節(jié)點數(shù)量過多會使BP 神經(jīng)網(wǎng)絡(luò)過于復(fù)雜，延長網(wǎng)絡(luò)訓(xùn)練以及后續(xù)計算的時間，缺少泛化能力，甚至降低識預(yù)測的準確率；而節(jié)點數(shù)量過少則會無法準確提取訓(xùn)練樣本特征，無法準確識別未進行訓(xùn)練的樣本，容錯性差。在設(shè)計過程中，確定隱含層節(jié)點數(shù)量k 主要通過經(jīng)驗公式計算，如下式所示：

式中m和n分別為輸入層和輸出層的神經(jīng)元個數(shù)，a為（1，10）范圍的常數(shù)。

BP 神經(jīng)網(wǎng)絡(luò)在進行訓(xùn)練以及最終的結(jié)果預(yù)測之前，需要初始化神經(jīng)網(wǎng)絡(luò)參數(shù)，即模型的所有神經(jīng)元之間的激活函數(shù)均設(shè)定小于1 的任意值，一般取（-1，1）之間的任意值。

其后選擇BP 神經(jīng)網(wǎng)絡(luò)的訓(xùn)練樣本，樣本盡量選擇輸入層所包含的范圍，才能獲得良好的網(wǎng)絡(luò)性能。根據(jù)以上原則，基本可以確定BP 神經(jīng)網(wǎng)絡(luò)的預(yù)測模型，訓(xùn)練步驟如圖3 所示。

圖3 BP 神經(jīng)網(wǎng)絡(luò)的預(yù)測步驟

圖中的正、負樣本分別指圖像中包含和未包含海冰的圖像。經(jīng)過參數(shù)的調(diào)整和訓(xùn)練后，可以確定BP 神經(jīng)網(wǎng)絡(luò)的激活函數(shù)采用線性函數(shù)purelin，隱含層神經(jīng)元個數(shù)為4，目標誤差為0.01，最大的訓(xùn)練次數(shù)為1000。

4.2 卡爾曼濾波算法設(shè)計

在使用航測資料進行海冰面積的預(yù)測時，需要采用卡爾曼濾波算法對圖像進行跟蹤處理。算法對圖像的求解過程如圖4 所示。

圖4 卡爾曼濾波算法對圖像的求解過程

圖像在t 時刻的先驗估計值Xt|t-1根據(jù)濾波方程計算：

式中Zt為描述t-1 和t 時刻圖像狀態(tài)關(guān)系的矩陣；Kt為控制圖像輸入的矩陣；μt為控制系統(tǒng)輸入的矩陣。其后根據(jù)先驗估計值，計算該值對應(yīng)的協(xié)方差矩陣，計算方式如下：

第三步為計算系統(tǒng)當前的卡爾曼增益Kt，計算方式如下：

式中Qt為描述t 時刻圖像狀態(tài)映射到無人機觀測區(qū)域的矩陣；Gt為按照高斯分布表示的協(xié)方差矩陣。第四步為計算系統(tǒng)更新后的后驗估計值Xt|t和協(xié)方差矩陣，分別如下式所示：

式中Zt為當前時刻觀察的圖像值。其后，系統(tǒng)將此時刻的圖像信息作為上一時刻，依次循環(huán)該求解過程。

4.3 海冰的結(jié)果預(yù)測和分析

為檢驗?zāi)Ｐ偷膶１娣e預(yù)測的精度，采用該模型對渤海海灣的海冰面積進行預(yù)測。對2019年12 月1 日至2020 年3 月8 日的海冰面積進行預(yù)測，同時采用遙感以及固定平臺對海冰的實際面積進行觀察，將預(yù)測值與實際的測量值進行對比，計算預(yù)測誤差。結(jié)果如表1 所示。

表1 海冰預(yù)測試驗結(jié)果對比

由表1 可知，對海冰面積的預(yù)測值和觀察值的結(jié)果以及海冰面積的變化趨勢基本一致，說明該模型可用于對海冰面積的預(yù)測，可以采用該模型為海洋運輸系統(tǒng)的正常運行提供理論依據(jù)。

5 結(jié)論

當我國的渤海和北黃海海域發(fā)生海冰時，會對航海運輸和油氣勘探等造成破壞，為保證海洋運輸系統(tǒng)的安全、減少并預(yù)防海冰可能會造成的災(zāi)害，本文主要進行了以下工作：

（1）建立了多個收集海冰相關(guān)資料的方法，包括衛(wèi)星遙感、雷達監(jiān)測、航測以及固定的監(jiān)測站等。

（2）對海冰的熱力學和動力學要素進行了分析。通過對冰面太陽輻射、長波輻射、感熱、冰內(nèi)熱傳導(dǎo)和冰下海洋熱通量等熱力學要素對海冰生消機理的影響進行分析，討論海冰厚度存在的條件；通過對拖曳力、海冰內(nèi)力等動力學因素進行分析，討論了海冰漂移的動力學方程。同時建立了海冰的連續(xù)方程和本構(gòu)方程。

（3）充分考慮了渤海海域海冰生消和漂移的特點，采用拉格朗日坐標系的光滑質(zhì)點法對海冰的演化過程進行了數(shù)值模擬。

（4）本文采用BP 神經(jīng)網(wǎng)絡(luò)對海冰面積進行預(yù)報，并采用卡爾曼濾波算法對預(yù)測結(jié)果進行跟蹤，以進一步提高預(yù)報精度。采用該模型對渤海海冰面積進行預(yù)測，結(jié)果表明該模型的精度較高，可用于對海冰面積進行預(yù)測，為海洋運輸系統(tǒng)的正常運行提供理論依據(jù)。