海冰工況下海洋運輸系統的仿真建模分析

張 明

（福建船政交通職業學院，福建 福州 350007）

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

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

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

1 海冰監測系統設計

1.1 海冰資料的收集

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

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

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

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

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

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

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

2.1 熱力學過程

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

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

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

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

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

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

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

2.2 動力學過程

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

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

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

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

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

2.3 連續方程

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

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

2.4 本構方程

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

式中σ0和σ1為海冰內部的兩個主應力；F0為海冰內部強度；e 為橢圓的主軸比。Hibler 方程中，海冰的應力和應變之間的關系如下：

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

式中F＇和A均為根據經驗確定的常數。

3 海冰數值模擬計算

3.1 海冰數值模擬計算方法的選擇

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

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

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

3.2 海冰數值的模擬計算

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

其中海冰內力這項可以寫為：

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

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

4 神經網絡對海冰面積的預報研究

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

4.1 人工神經網絡模型的設計

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

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

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

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

BP 神經網絡在進行訓練以及最終的結果預測之前，需要初始化神經網絡參數，即模型的所有神經元之間的激活函數均設定小于1 的任意值，一般取（-1，1）之間的任意值。

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

圖3 BP 神經網絡的預測步驟

圖中的正、負樣本分別指圖像中包含和未包含海冰的圖像。經過參數的調整和訓練后，可以確定BP 神經網絡的激活函數采用線性函數purelin，隱含層神經元個數為4，目標誤差為0.01，最大的訓練次數為1000。

4.2 卡爾曼濾波算法設計

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

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

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

式中Zt為描述t-1 和t 時刻圖像狀態關系的矩陣；Kt為控制圖像輸入的矩陣；μt為控制系統輸入的矩陣。其后根據先驗估計值，計算該值對應的協方差矩陣，計算方式如下：

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

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

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

4.3 海冰的結果預測和分析

為檢驗模型的對海冰面積預測的精度，采用該模型對渤海海灣的海冰面積進行預測。對2019年12 月1 日至2020 年3 月8 日的海冰面積進行預測，同時采用遙感以及固定平臺對海冰的實際面積進行觀察，將預測值與實際的測量值進行對比，計算預測誤差。結果如表1 所示。

表1 海冰預測試驗結果對比

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

5 結論

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

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

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

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

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