低密集度浮碎冰數(shù)值生成方法研究

張忠宇，谷家揚，王志東，陶延武，渠基順，李文娟

（江蘇科技大學(xué) a.海洋裝備研究院；b.船舶與海洋工程學(xué)院，江蘇 鎮(zhèn)江 212003）

0 引 言

全球氣候變暖導(dǎo)致極地海區(qū)海冰發(fā)生融化，且波浪等因素推動海冰向外漂移，使得冰區(qū)密集度沿著緯度逐步遞減，至冰緣區(qū)可降低至20%左右。低密集度冰況使得極地商業(yè)航行及油氣資源開發(fā)作業(yè)變?yōu)榭尚?。冰緣區(qū)海冰主要以浮碎冰的形式存在[1]，其中直徑小于2 m的海冰稱為碎冰，大于2 m的海冰稱為浮冰。浮碎冰的形成機理較為復(fù)雜，浮冰的演變是動力學(xué)和熱力學(xué)的共同作用。首先夏季溫度場會使冰的內(nèi)力變小，冰排發(fā)生破裂形成大尺度浮冰，浮冰在溫度場的作用下會繼續(xù)發(fā)生融化。Perovich 等[2]通過一系列直升機航拍照片研究了夏季過程中浮碎冰密集度演變及浮碎冰的融化特性。Steele 等[3-4]指出大直徑浮冰的融化速率明顯大于直徑小于30 m 的浮冰。因此，對于冰緣區(qū)尤其是低密集度冰況的海域，海冰的形式主要以小尺度浮碎冰為主。若在此海域航行作業(yè)，浮碎冰區(qū)的船舶阻力性能也是極地船舶設(shè)計的重要考量因素之一。

為了研究浮碎冰區(qū)的極地船舶阻力性能，首先需確定相關(guān)海域的海冰大小及分布。直接觀測法是目前最主要的技術(shù)手段，通過航拍得到目標(biāo)海域的浮碎冰圖像，對圖像進(jìn)行處理即可?；Ｓ虻母∷楸植?。為了便于統(tǒng)計分析規(guī)律，Rothrock和Thorndike[5]提出了浮碎冰等效直徑分布的概念并討論了衡量浮碎冰直徑分布的方法。之后，一系列的研究工作[6-9]采用衛(wèi)星或航拍等手段研究了浮碎冰等效直徑的分布，并發(fā)現(xiàn)浮碎冰的數(shù)量隨其等效直徑呈指數(shù)分布，冪指數(shù)為-β。在鄂霍次克海南部，Toyota等[10-11]通過航拍和船載攝像系統(tǒng)研究了鄂霍次克海南部海域的浮碎冰分布，當(dāng)浮冰直徑在25 m到50 m 之間時，β的取值減小到1.15 左右。Toyata 等認(rèn)為這種現(xiàn)象是由于浮碎冰的膨脹導(dǎo)致的。而Steer和Lu等[12-13]通過研究航拍照片發(fā)現(xiàn)，浮碎冰數(shù)量分布在幾十米直徑范圍內(nèi)會發(fā)生急劇的變化，但對于小尺度浮碎冰區(qū)域，β的取值在1.9左右，這與Toyota等之前的結(jié)論相差較大。為了分析原因以掌握浮碎冰的分布，Toyota 等[14]又通過航拍研究了Weddell海域浮碎冰的分布情況，浮碎冰的直徑在1 m到100 m 之間。冬季由于浮碎冰受膨脹因素的影響較小，所以β的取值較小。對于鄂霍次克海，β∈[1.0,3.0]。但由于極地惡劣的環(huán)境，難以到達(dá)目標(biāo)海域進(jìn)行數(shù)據(jù)采集，且航拍的成本較高，相關(guān)資料相對較少。

可見，數(shù)據(jù)的缺乏導(dǎo)致浮碎冰分布的理論還處在逐步完善階段。另外一種方法是基于浮碎冰隨機生成理論，對目標(biāo)海域進(jìn)行數(shù)值劃分得到浮碎冰分布。朱紅日等[15]基于二維Voronoi 圖方法對碎冰域進(jìn)行隨機切割以生成碎冰區(qū)中冰塊的幾何形態(tài)，從而確定碎冰區(qū)的初始分布場。但這種方法生成的碎冰直徑呈正態(tài)分布，與實際觀測的規(guī)律不符。Yulmetov 等[16]開發(fā)了一套生成指數(shù)分布的碎冰程序，但其只是簡要介紹了基本理論，并未詳細(xì)闡述浮碎冰的生成過程?？傊?，低密集度碎冰無論從物理生成機理還是數(shù)值生成方法均處于起步階段，相關(guān)資料較少。為了得到合理的浮碎冰場分布，本文在Yulmetov 等[16]工作的基礎(chǔ)上，將Voronoi 圖方法、概率理論等進(jìn)行有機結(jié)合，形成給定密集度下浮碎冰的生成方法，為浮碎冰場的數(shù)值建模提供合理的輸入。

1 低密集度浮碎冰隨機生成理論

1.1 目標(biāo)海域數(shù)值劃分的Voronoi圖方法

通過與實際航拍的浮碎冰場進(jìn)行對比，浮碎冰的形狀與Voronoi圖類似。Voronoi圖[17]是1908年由俄羅斯數(shù)學(xué)家Voronoi 提出，可對目標(biāo)域完成幾何隨機劃分。其在目標(biāo)域生成一組隨機點，按照最鄰近原則劃分平面，每個點與它的最近鄰區(qū)域相關(guān)聯(lián)，采用Delaunay 三角剖分法生成二維Voronoi圖，將目標(biāo)域分割為一組泰森多邊形集合?？梢姡瑸榱说玫絍oronoi圖，首先需生成二維隨機點，傳統(tǒng)的做法是通過隨機生成器生成偽隨機數(shù)，這些隨機數(shù)滿足泊松概率分布。但對于二維隨機數(shù)，在生成的過程中，新生成的隨機點與前面的點是獨立的，兩個隨機點間距離可能極大或極小，生成尺度極大或極小的浮碎冰，不符合實際海況。為了避免這種情況發(fā)生，本文采用Simple Sequential Inhibition(SSI)方法[18]生成二維隨機點，該方法的思想是首先生成一系列隨機點，給定隨機點之間的最大距離為s，若隨機點間的距離r與最大距離s的比值小于某一數(shù)值δ，則將該隨機點剔除。其中，最大距離的定義為

式中，A為碎冰域的面積，N為隨機點的數(shù)量。為了定性分析最大距離的影響，定義正則參數(shù)δ=r/s，該參數(shù)的影響將在后文中進(jìn)行討論。

1.2 量化碎冰尺寸分布概率函數(shù)

在采用Voronoi圖對碎冰域進(jìn)行劃分后，浮碎冰的數(shù)量隨其等效直徑D分布函數(shù)呈正態(tài)分布。而依據(jù)航拍的試驗數(shù)據(jù)可知，浮碎冰的數(shù)量隨其等效直徑D呈指數(shù)分布，如下式所示：

式中，NI為浮碎冰的數(shù)量，α和β為常數(shù)。

由式（2）可知，其累積分布函數(shù)F(D)為

式中，Dmin為浮碎冰的最小直徑。

由累計分布函數(shù)的性質(zhì)可知：

式中，Dmax為浮碎冰的最大直徑。

將式（4）代入式（3）可得碎冰場的概率密度分布函數(shù)f(D)為

為了使得浮碎冰分布服從指數(shù)分布，需對Voronoi圖中的泰森多邊形進(jìn)行量化處理。由泰森多邊形的性質(zhì)可知：

由隨機變量商的性質(zhì)，可得量化因子的概率密度等效函數(shù)為

式中，I為關(guān)于的積分，對于給定的β和較大的浮碎冰數(shù)量N（隨機點的數(shù)量），其積分為常數(shù)。但當(dāng)浮碎冰數(shù)量較大時，影響數(shù)值的計算效率，對于其取值，本文將在后續(xù)進(jìn)行討論。

通過累積分布函數(shù)的性質(zhì)對式（7）進(jìn)行推導(dǎo)[16]可知：

式中，x為在[0,1]內(nèi)的均勻分布隨機數(shù)，本文通過Matlab的隨機生成器獲得。

通過式（8）得到量化因子，進(jìn)而得到量化后的浮碎冰直徑。通過給定的密集度，依次選取不同直徑區(qū)間的浮碎冰，使其滿足密集度的要求。

但量化浮碎冰尺寸后的浮碎冰場存在數(shù)值重疊，對于未受擾動的低密集度浮碎冰場，其不存在重疊現(xiàn)象，需對重疊區(qū)域進(jìn)行數(shù)值消除。本文主要采用分離軸定理進(jìn)行處理，遍歷泰森多邊形的每條邊，將重疊最小的軸作為分離軸，讓重疊的泰森多邊形沿分離軸向相反方向分開即可，如圖1所示，具體過程詳見文獻(xiàn)[19]，本文不再詳細(xì)贅述。

圖1 浮碎冰重疊區(qū)域分離示意圖Fig.1 Separation diagram of overlapped broken ices

2 數(shù)值結(jié)果分析與討論

本章在低密集度浮碎冰隨機生成理論的基礎(chǔ)上，開展浮碎冰場的數(shù)值模擬，在驗證方法有效性的同時，開展相關(guān)變量對浮碎冰分布特性的影響研究，得到低密集度下浮碎冰合理的空間分布。為了定性分析，數(shù)值模擬過程中目標(biāo)海域均取為100 m×100 m，最大海冰直徑為10 m、最小海冰直徑為1 m。

2.1 隨機點數(shù)N的影響研究

從低密集度浮碎冰隨機生成理論可知，隨機點數(shù)N主要影響Voronoi圖生成的泰森多邊形數(shù)量及量化過程中的式中的積分I。為了定性分析隨機點數(shù)的影響，分別取N=3 000、4 000、5 000、10 000、20 000，開展密集度為0.3的浮碎冰場數(shù)值模擬，β取為2.5，正則參數(shù)δ取為0.1。

圖2 為不同隨機點數(shù)下等效直徑的隨機概率函數(shù)。從圖中可以看到，由于數(shù)值模擬的目標(biāo)海域面積恒定，隨機點數(shù)的增加顯然會導(dǎo)致生成的泰森多邊形的尺寸較小。隨機點數(shù)為20 000 時，等效直徑的最大概率出現(xiàn)在0.78 m左右。隨著隨機點數(shù)減少，等效直徑的概率密度峰值右移，概率密度曲線也會隨之變“扁平”。對于目標(biāo)海域浮碎冰場的數(shù)值?；?，傳統(tǒng)的處理方法是將Voronoi 圖按密集度進(jìn)行等比例分布，以達(dá)到給定密集度的目的。但通過圖2可知，這種處理方法獲得的浮碎冰直徑基本上呈正態(tài)分布，與實驗測量擬合的結(jié)果不一致。

圖2 不同隨機點數(shù)下等效直徑的隨機概率函數(shù)Fig.2 Random probability function of equivalent diameter under different numbers of random points

圖3 為量化后不同隨機點下浮碎冰的數(shù)值分布，其中紅色實線為通過公式得到的浮碎冰直徑的理論概率密度函數(shù)，藍(lán)色直方圖表征數(shù)值模擬的浮碎冰概率密度函數(shù)。從圖中可以看到，不同隨機點下數(shù)值結(jié)果均與理論值吻合相對良好，從而驗證了本文方法的合理性。雖然量化前浮碎冰直徑的隨機概率函數(shù)不同，但量化后隨機點數(shù)對數(shù)值結(jié)果基本沒有影響。顯然，對量化過程中式中的積分I，假設(shè)其為常數(shù)也是成立的，隨機點數(shù)對其影響較小。表1 給出了各隨機點數(shù)下不同直徑區(qū)間內(nèi)浮碎冰的數(shù)量，五種隨機點數(shù)下生成的浮碎冰總數(shù)分別為697 塊、691 塊、701 塊、696 塊和691 塊。從表中可以看出，在大于4 m 的直徑區(qū)間生成的浮碎冰數(shù)量均相同，僅在較小的直徑區(qū)間內(nèi)浮碎冰的數(shù)量略有差異。以隨機點數(shù)3 000下浮碎冰的概率密度為參照，表1還給出了各隨機點下不同直徑區(qū)間內(nèi)浮碎冰概率密度誤差。隨著隨機點數(shù)的增加，概率密度誤差整體仍較小。可見，基于本文數(shù)值方法生成的浮碎冰場隨著隨機點數(shù)增加是收斂的。為了提高計算效率，本文在后續(xù)數(shù)值模擬過程中隨機點數(shù)均取為3 000。

表1 不同隨機點數(shù)下浮碎冰數(shù)量及概率密度誤差Tab.1 Quantityandprobabilitydensityerrorsofbrokenicesunderdifferentnumbersofrandompoints

圖3 量化后不同隨機點數(shù)下浮碎冰數(shù)值分布Fig.3 Numerical distribution of quantified broken ices under different numbers of random points

2.2 正則參數(shù)δ的影響研究

為了定性分析正則參數(shù)δ的影響，分別取δ=0、0.05、0.1、0.2，開展密集度為0.3 的浮碎冰場分布數(shù)值模擬，β取為2.5。圖4為不同正則參數(shù)δ下浮碎冰數(shù)值分布，橫坐標(biāo)為浮碎冰的等效直徑，縱坐標(biāo)為浮碎冰的數(shù)量，曲線上的標(biāo)記點表征數(shù)值模擬中小于該直徑區(qū)間的浮碎冰數(shù)量。由于隨機點在生成的過程中是隨機的，因而本文對每個正則參數(shù)開展了5 次重復(fù)性數(shù)值試驗。從圖中可以看到，當(dāng)δ=0時二維點的生成是SSI方法的一種特例，即在初始二維點生成時，通過隨機生成器生成偽隨機數(shù)，新生成的隨機點與前面的點是獨立的，兩個隨機點間距離可能為極大或極小，不存在篩選過程。由于隨機點的生成是完全隨機的，導(dǎo)致重復(fù)生成的浮碎冰數(shù)量出現(xiàn)一定的離散度。從圖4（a）的數(shù)值結(jié)果中可以看到，5 次數(shù)值結(jié)果的離散度相對較大，不利于后續(xù)對極地船舶冰阻力特性的定性研究。隨著正則參數(shù)的增加（圖4（b）），5次數(shù)值結(jié)果的離散度相對減小。當(dāng)δ=0.1時，數(shù)值結(jié)果的離散度最小，各直徑區(qū)間的浮碎冰數(shù)量基本相當(dāng)。這是因為δ表征的是SSI 生成方法中隨機點間的距離與最大距離之間的比值，通過控制δ的取值，刪掉一些不符合要求的隨機點，從而一定程度上降低了碎冰生成的隨機性，進(jìn)而間接控制生成浮碎冰的大小。但當(dāng)δ=0.2 時，數(shù)值結(jié)果的離散度尤其是小直徑區(qū)間的離散度反而進(jìn)一步增大。顯然δ的增加，會刪掉一些處于小直徑區(qū)間的隨機點，使得生成浮碎冰的離散度增大。朱紅日等[15]研究發(fā)現(xiàn)，當(dāng)δ逐漸繼續(xù)增加并趨于1時，所有隨機點之間的距離趨近于相等。因而，此時生成的浮碎冰具有相同的幾何形態(tài)，生成的浮碎冰直徑難以滿足指數(shù)分布。綜合以上分析，本文在數(shù)值模擬過程中正則參數(shù)δ均取為0.1。

圖4 不同正則參數(shù)下浮碎冰數(shù)值分布Fig.4 Numerical distribution of broken ices under different regular parameters

圖5 給出了五次重復(fù)性數(shù)值模擬下浮碎冰的空間分布（δ=0.1）。從圖中可以看到，雖然五次重復(fù)性數(shù)值模擬過程中各直徑區(qū)間的浮碎冰數(shù)量基本相當(dāng)，但在理論推導(dǎo)及數(shù)值計算過程中，未有變量控制浮碎冰的空間分布，因而五次重復(fù)性數(shù)值試驗中最大海冰的空間位置各不相同。這種空間位置的不確定性對極地船舶冰阻力性能的影響規(guī)律，將是未來研究的一個重點。

圖5 δ=0.1下浮碎冰的空間分布Fig.5 Spatial distribution of broken ices at δ=0.1

2.3 參數(shù)β的影響研究

由公式可知，參數(shù)β是影響浮碎冰分布概率密度函數(shù)的一個重要參數(shù)。為了定性分析參數(shù)β的影響，分別取β=2.5、3.0、3.5、4.0，開展密集度為0.3的浮碎冰場數(shù)值模擬。

圖6為不同參數(shù)β下碎冰的數(shù)值分布。從圖中可以看到，參數(shù)β是影響浮碎冰分布的一個重要因素。隨著參數(shù)β的增加，較大直徑區(qū)間的浮碎冰數(shù)量迅速下降，這從圖7 中也可以看出。為了保證目標(biāo)海域的密集度不變，較小直徑區(qū)間的浮碎冰數(shù)量迅速增加。圖7（d）中較小直徑區(qū)間的浮碎冰數(shù)量達(dá)到776塊，目標(biāo)海域被小尺度碎冰覆蓋?？梢?，為了較為合理地數(shù)值?；繕?biāo)海域，參數(shù)β的取值至關(guān)重要。從相關(guān)文獻(xiàn)中可知，通過對航拍得到的實驗數(shù)據(jù)進(jìn)行擬合發(fā)現(xiàn)：不同海域參數(shù)β的取值不一樣，低密集度區(qū)域參數(shù)β的取值通常在2.5左右。

圖6 不同參數(shù)β 下浮碎冰數(shù)值分布Fig.6 Numerical distribution of broken ices under different β

圖7 不同參數(shù)β下碎冰空間分布Fig.7 Spatial distribution of broken ices under different β

3 結(jié) 論

本文基于低密集度浮碎冰隨機生成理論，開展浮碎冰場的數(shù)值模擬，在驗證方法有效性的同時，開展隨機點數(shù)N、正則參數(shù)δ、參數(shù)β等變量對浮碎冰分布特性的影響研究，合理確定相關(guān)變量的取值。主要的結(jié)論如下：

（1）隨機點數(shù)雖然會改變Voronoi 圖中泰森多邊形等效直徑的概率密度函數(shù)，但浮碎冰直徑經(jīng)過量化后，隨機點數(shù)對數(shù)值結(jié)果基本沒有影響；

（2）通過重復(fù)性數(shù)值實驗可知，正則參數(shù)δ影響數(shù)值結(jié)果的離散度。當(dāng)δ=0.1時，數(shù)值結(jié)果的離散度最??；

（3）參數(shù)β是影響浮碎冰分布概率密度函數(shù)的一個重要參數(shù)。隨著參數(shù)β的增加，較大直徑區(qū)間的碎冰數(shù)量迅速下降，較小直徑區(qū)間的碎冰數(shù)量迅速增加。