基于局部體樣條函數(shù)的海洋鹽度場(chǎng)建模

楊振發(fā)，萬(wàn) 剛，李 鋒，張宗佩

（信息工程大學(xué)，河南 鄭州450001）

隨著時(shí)代的發(fā)展，海洋對(duì)各國(guó)的價(jià)值日益凸顯，戰(zhàn)略地位不斷提升，特別是進(jìn)入21世紀(jì)后，各國(guó)正加快步伐，建立“數(shù)字海洋”。海水鹽度是海洋水文要素的重要組成部分，其在垂直方向和水平方向的變化，直接影響艦艇聲納的工作距離和艦艇的航行速度和日數(shù)。

將海洋鹽度的數(shù)據(jù)場(chǎng)進(jìn)行可視化表達(dá)既可以形象逼真的展現(xiàn)數(shù)據(jù)本身的結(jié)構(gòu)特征，又可以通過(guò)交互等操作輔助人員做出合理的判斷，而數(shù)據(jù)場(chǎng)科學(xué)可視化的關(guān)鍵在于數(shù)據(jù)的建模。針對(duì)海洋鹽度場(chǎng)的散亂大數(shù)據(jù)場(chǎng)建模，國(guó)內(nèi)外文獻(xiàn)歸納的建模方法大致可以分為兩種：一種是層次化方法，分為基于多層B樣條插值方法［1-2］和多分辨率方法［3］，前者正在解決雙變量數(shù)據(jù)建模的問(wèn)題，而后者作為提高大數(shù)據(jù)可視化的交互手段，主要用于復(fù)雜模型的簡(jiǎn)化；另一種是局部化方法，包括基于體數(shù)據(jù)的最小范數(shù)網(wǎng)絡(luò)法［4］和局部體樣條函數(shù)法［5］，前者建模精度高，但實(shí)現(xiàn)需要一個(gè)四面體化的算法和一個(gè)復(fù)雜的迭代方法，而后者將大數(shù)據(jù)場(chǎng)局部化成小區(qū)域，采用適用于中小數(shù)據(jù)量的體樣條函數(shù)建模方法，該方法精度高，且易于實(shí)現(xiàn)，但關(guān)鍵問(wèn)題在于解決如何將區(qū)域適當(dāng)?shù)胤指睢?/p>

針對(duì)以上適用于海洋鹽度散亂數(shù)據(jù)場(chǎng)建模方法的分析，本文在介紹海洋水文要素中的鹽度數(shù)據(jù)的特點(diǎn)后，采用局部體樣條函數(shù)法對(duì)海洋鹽度場(chǎng)進(jìn)行建模，并利用非線性規(guī)劃的遺傳算法，解決局部體樣條函數(shù)法在散亂鹽度場(chǎng)建模中的區(qū)域劃分優(yōu)化組合問(wèn)題，最后通過(guò)最大密度投影算法實(shí)現(xiàn)建模數(shù)據(jù)的可視化表達(dá)，為海洋戰(zhàn)場(chǎng)環(huán)境表達(dá)與分析提供較為形象、直觀的三維鹽度場(chǎng)仿真手段。

1 海水鹽度數(shù)據(jù)特點(diǎn)與來(lái)源

在海洋化學(xué)中，鹽度的確切定義為：“1kg海水中，所有碳酸鹽轉(zhuǎn)化為氧化物，溴、碘以氯置換，所有的有機(jī)物被氧化之后所含全部固體物質(zhì)的總克數(shù)”，單位為g／kg［6］。海洋表面鹽度分布的總規(guī)律是：從亞熱帶海區(qū)向高、低緯遞減，并形成馬鞍形曲線，赤道附近鹽度較低，南北緯20°附近的鹽度最高，然后隨緯度的增加而降低。海水鹽度垂直分布規(guī)律是：在北緯40°到南緯50°之間，是鹽度垂直變化最大、最復(fù)雜的地區(qū)（見(jiàn)圖1）。從海面到150m深度上鹽度高而均勻，最大鹽度值一般出現(xiàn)在100～300m之間，深層水的鹽度分布最均勻，鹽度值比表層水低、比中層水高。

圖1 海水鹽度場(chǎng)的分布規(guī)律

海水鹽度場(chǎng)是海洋環(huán)境中不可見(jiàn)的水文要素，與海底地貌等數(shù)據(jù)場(chǎng)存在很大的區(qū)別。鹽度數(shù)據(jù)主要來(lái)源于艦船的調(diào)查數(shù)據(jù)，其特點(diǎn)是數(shù)據(jù)分布相對(duì)離散、不規(guī)則，屬于典型的大數(shù)據(jù)場(chǎng)的散亂體數(shù)據(jù)。而僅有這些離散的、不規(guī)則的數(shù)據(jù)難以進(jìn)行鹽度場(chǎng)的三維可視化連續(xù)表達(dá)。

2 海水鹽度場(chǎng)的建模

散亂體數(shù)據(jù)的插值問(wèn)題［7］在20世紀(jì)60年代開(kāi)始引起人們的注意，經(jīng)過(guò)近半個(gè)世紀(jì)的研究，已有適用于不同需求和數(shù)據(jù)類型的建模方法，且大都只是適用于中小規(guī)模的數(shù)據(jù)場(chǎng)建模算法。Nielson在文獻(xiàn)［7］中通過(guò)實(shí)驗(yàn)證明當(dāng)數(shù)據(jù)量在200～500以內(nèi)時(shí)，體樣條函數(shù)法對(duì)于解決數(shù)據(jù)場(chǎng)的建模問(wèn)題相當(dāng)有效，而規(guī)模超過(guò)500后，求解得到的線性方程組系數(shù)矩陣的條件數(shù)將存在問(wèn)題，使得解的魯棒性和建模精度大大降低。鑒于體樣條函數(shù)法［8］在中小規(guī)模散亂數(shù)據(jù)場(chǎng)的建模中具有擬合精度高、易于實(shí)現(xiàn)等優(yōu)點(diǎn)，因此可以將大規(guī)模數(shù)據(jù)場(chǎng)適當(dāng)分區(qū)，分別構(gòu)造建模函數(shù)，再賦予其權(quán)重并將各區(qū)“拼”成一個(gè)整體。而這種局部的體樣條函數(shù)法的建模精度關(guān)鍵在于區(qū)域分割的優(yōu)化程度，本文采用非線性規(guī)劃的遺傳算法來(lái)解決這個(gè)關(guān)鍵問(wèn)題。圖2為鹽度場(chǎng)數(shù)據(jù)的建模流程。

圖2 鹽度場(chǎng)建模流程

2.1 局部的體樣條函數(shù)法

為了便于建模操作，假設(shè)鹽度數(shù)據(jù)的采樣空間為一規(guī)則的長(zhǎng)方體Ω＝｛P＝（x，y，z）｜a≤x≤b，c≤y ≤d，e≤z≤f｝。采樣點(diǎn) Pi＝ （xi，yi，zi）（1≤i≤n）的鹽度數(shù)據(jù)為fi。故對(duì)鹽度數(shù)據(jù)場(chǎng)的建模即構(gòu)造一個(gè)在Ω上的連續(xù)函數(shù)F（P）＝F（x，y，z），使Pi處采樣值滿足F（Pi）＝F（xi，yi，zi）＝fi（1≤i≤n）。實(shí)現(xiàn)步驟如下：

1）將長(zhǎng)方體區(qū)域分塊。將Ω的3個(gè)方向分別劃分為nx，ny，nz個(gè)分塊，分塊的支集為 Ωijk（0≤i≤nx，0≤j≤ny，0≤k≤nz），保證支集中的數(shù)據(jù)點(diǎn)數(shù)大致相同且不超過(guò)500，這將在下一節(jié)中采用非線性規(guī)劃的遺傳算法予以解決。

2）構(gòu)造分塊區(qū)域Ωijk的權(quán)函數(shù)Wijk（P）。權(quán)函數(shù)滿足的條件如下：

a.當(dāng)P?Ωijk時(shí)，Wijk（P）＝0；

文獻(xiàn)（3）采用Hermite的方法，給出具體求解權(quán)函數(shù)的方法。

3）構(gòu)造分塊的建模函數(shù)Fijk（P）。當(dāng)數(shù)據(jù)場(chǎng)的數(shù)據(jù)個(gè)數(shù)小于500時(shí)，研究表明多重二次型法［9］具有擬合精度高等優(yōu)點(diǎn)，這里采用此法構(gòu)造函數(shù)。建模函數(shù)為

式中：‖P-Pi‖為插值點(diǎn)到采樣點(diǎn)的距離函數(shù)；K2為實(shí)驗(yàn)參數(shù)，這里設(shè)為0.1；系數(shù)ai可由下面的線性方程組的解確定

2.2 非線性規(guī)劃的遺傳算法實(shí)現(xiàn)區(qū)域的優(yōu)化分割

采用局部體樣條函數(shù)法的前提和關(guān)鍵是解決大區(qū)域的分塊問(wèn)題，即如何將分塊后的支集包含的采樣點(diǎn)數(shù)大致相同且限制在500以下。為了解決這樣一個(gè)帶約束的組合優(yōu)化問(wèn)題，本文采用非線性規(guī)劃的遺傳算法［10］。

非線性規(guī)劃問(wèn)題是研究一個(gè)n元實(shí)函數(shù)在一組等式或不等式的約束條件下的極值問(wèn)題，它具有局部搜索能力強(qiáng)等優(yōu)點(diǎn)，但易于陷入局部次優(yōu)解。遺傳算法是一類借鑒生物界自然選擇和遺傳機(jī)制的隨機(jī)搜索算法，它從隨機(jī)產(chǎn)生的初始解開(kāi)始搜索，通過(guò)一定的選擇、交叉、變異等操作逐步迭代產(chǎn)生新解。群體中的每個(gè)個(gè)體代表問(wèn)題的一個(gè)解，稱為染色體，染色體的好壞用適應(yīng)度值來(lái)衡量，根據(jù)適應(yīng)度的好壞從上一代中選擇一定數(shù)量的優(yōu)秀個(gè)體，通過(guò)交叉、變異形成下一代群體。經(jīng)過(guò)若干代進(jìn)化后，算法收斂于最好的染色體，即為最優(yōu)解或次優(yōu)解。遺傳算法在搜索時(shí)從問(wèn)題解的一個(gè)集合開(kāi)始，而不是單個(gè)個(gè)體，可大大減少陷入局部最小的可能性。結(jié)合以上兩種算法的優(yōu)缺點(diǎn)，非線性規(guī)劃的遺傳算法解決局部體樣條函數(shù)法區(qū)域劃分問(wèn)題的算法流程如圖3所示。實(shí)現(xiàn)的具體步驟為：

1）種群初始化。要使每個(gè)劃分的數(shù)據(jù)量滿足r∈（r1，r2），以X 軸向?yàn)槔齞x＝b-a，m＝ ［n／r1］，其劃分Hx編碼為dx｝，xi表示第i段的區(qū)間長(zhǎng)度，令H＝Hx×Hy×Hz，則J＝J（H）為劃分H 的一個(gè)可能解。對(duì)劃分分別隨機(jī)產(chǎn)生3個(gè)劃分?jǐn)?shù)量u，v，w∈［1，m-1］，使［u×v×w／m］＝1，令xi＝［dx／u］（0≤i≤u-1），其余碼元值為0。

4）交叉操作算子。隨機(jī)選取兩個(gè)個(gè)體，通過(guò)染色體的交換組合，把父輩的優(yōu)秀特征傳遞給子輩，從而產(chǎn)生優(yōu)秀個(gè)體，第m個(gè)染色體和第n個(gè)染色體在l位的交叉操作為

其中，b是區(qū)間［0，1］的隨機(jī)數(shù)。

5）變異操作算子。其目的是維護(hù)種群多樣性，糾正遺傳算法過(guò)早收斂的缺點(diǎn)。從種群中隨機(jī)選擇第i個(gè)體，第j個(gè)編碼編譯操作為

其中amax是染色體aij的上界，amin是染色體的下界；f（g）＝r（1-g／Gmax）2，r是隨機(jī)數(shù)，g是當(dāng)前的迭代次數(shù)，Gmax是最大進(jìn)化次數(shù)。

6）非線性規(guī)劃尋優(yōu)操作。對(duì)進(jìn)化代數(shù)進(jìn)行判斷，這里選擇代數(shù)為5，當(dāng)滿足所定的代數(shù)時(shí)，為了加快進(jìn)化，利用Matlab的規(guī)劃工具箱函數(shù)進(jìn)行局部搜索最優(yōu)值，把所得到的局部值作為新個(gè)體繼續(xù)進(jìn)行進(jìn)化操作。

7）終止判斷條件。對(duì)所得的各代染色體進(jìn)行適應(yīng)度函數(shù)計(jì)算并比較，若為最大則終止進(jìn)化；否則返回步驟（2）。

圖3 非線性規(guī)劃的遺傳算法實(shí)現(xiàn)流程

3 實(shí)驗(yàn)及結(jié)果

從國(guó)家海洋局網(wǎng)站上獲取大規(guī)模散亂的海洋鹽度數(shù)據(jù)，利用Matlab軟件的線性規(guī)劃工具箱和遺傳工具箱，編寫(xiě)M文件程序?qū)崿F(xiàn)對(duì)數(shù)據(jù)場(chǎng)的基于非線性規(guī)劃遺傳算法的區(qū)域分塊。如圖4所示，圖4（a）（b）（c）為散亂數(shù)據(jù)空間的三個(gè)軸向數(shù)據(jù)分塊結(jié)果，整個(gè)區(qū)域被分成共48個(gè)小區(qū)域，每個(gè)區(qū)域的數(shù)據(jù)量在400～460之間。圖4（d）為某一鹽度層的插值函數(shù)三維可視化結(jié)果。

圖4 局部體樣條函數(shù)建模的實(shí)驗(yàn)結(jié)果

體繪制［11］研究的是含有物體內(nèi)部信息的體數(shù)據(jù)的表示、操作和顯示的問(wèn)題，較之于傳統(tǒng)的面繪制方法，體繪制方法能表達(dá)空間數(shù)據(jù)場(chǎng)內(nèi)部的數(shù)據(jù)特征和關(guān)系，因此特別適用于鹽度場(chǎng)數(shù)據(jù)的可視化表達(dá)。最大密度投影［12］（Maximum Intensity Projection）算法是一種被廣泛應(yīng)用于醫(yī)學(xué)CT和MRT圖像處理的體繪制算法。該方法使用光線跟蹤算法［13］跟蹤二維屏幕上的每一個(gè)像素所發(fā)射出的每一條射線與三維空間數(shù)據(jù)相交的每個(gè)體素，逐個(gè)進(jìn)行遍歷，找出每條射線中的最大數(shù)據(jù)值，無(wú)需進(jìn)行光照和累加這兩個(gè)過(guò)程，因此具有思路簡(jiǎn)單，便于編程實(shí)現(xiàn)等優(yōu)點(diǎn)。圖5給出MIP算法實(shí)現(xiàn)海洋鹽度標(biāo)量數(shù)據(jù)場(chǎng)體繪制的流程圖，圖6為MIP算法實(shí)現(xiàn)的某海域鹽度場(chǎng)體視化結(jié)果。

利用上述介紹的散亂數(shù)據(jù)場(chǎng)建模得到的規(guī)則鹽度場(chǎng)，利用OpenGL圖形接口的可編程著色器功能，在Intel（R）Core（TM）i5CPU 2.53GHz；內(nèi)存2.00GB；HD5730獨(dú)立顯卡的計(jì)算機(jī)上實(shí)現(xiàn)基于最大密度投影算法的鹽度場(chǎng)數(shù)據(jù)體繪制。

圖5 最大密度投影算法實(shí)現(xiàn)流程

圖6 MIP算法實(shí)現(xiàn)的某海域鹽度場(chǎng)體視化結(jié)果

4 結(jié) 論

針對(duì)海洋水文要素?cái)?shù)據(jù)量大且散亂等特點(diǎn)，本文采用局部的體樣條函數(shù)插值算法對(duì)標(biāo)量數(shù)據(jù)場(chǎng)進(jìn)行插值建模，為了解決散亂大數(shù)據(jù)場(chǎng)局部化建模時(shí)需滿足的組合優(yōu)化條件才能保證插值精度的問(wèn)題，本文利用具有較強(qiáng)全局搜索能力的非線性規(guī)劃的遺傳算法，較好地解決大數(shù)據(jù)場(chǎng)區(qū)域適當(dāng)分塊的問(wèn)題，最后通過(guò)Matlab工具箱編程實(shí)現(xiàn)區(qū)域的組合優(yōu)化分割，并將得到的規(guī)則數(shù)據(jù)場(chǎng)通過(guò)最大密度投影算法予以可視化表達(dá)，為海洋鹽度場(chǎng)的數(shù)據(jù)仿真提供較為可靠的仿真方法與手段。

［1］朱振華.雙三次B樣條曲面法在聲吶測(cè)量數(shù)據(jù)空白填補(bǔ)中的作用［J］.測(cè)繪工程，2013，22（1）：75-77.

［2］周玉娟，岳桂昌.二次B樣條修正VTEC多項(xiàng)式模型研究［J］.測(cè)繪工程，2013，22（3）：41-43.

［3］STAADT O G，GROSS M.Progressive tetrahedralizations［A］.In：Proceedings of IEEE Visualization98，North Carolina，1998.397-402.

［4］NIELSON G M，TVEDT J.Comparing methods of interpolation for scattered volumetric data［J］.Advanced Techniques for Scientific Visualization，ACM SIGGRAPH’94，II-99.

［5］楊冠杰，耿震，孫菁.散亂體數(shù)據(jù)可視化：海洋水團(tuán)分析的新途徑［J］.海洋通報(bào)，2000，19（2）：66-74.

［6］韓慶，李本江，高峰，等.海洋小百科全書(shū)［M］.中山：中山大學(xué)出版社，2012.

［7］SHEPARD D.A two dimensional interpolation function for irregularly space data［J］.Proceedings of ACM 23“National Conference，1968，517-524.

［8］NIELSON G M.Scattered data modeling［J］.IEEE CG＆A，1993，13（1）：60-70.

［9］KREGLOS O，HAMANN B.Data structures for optimizing linear spline approximations［J］.Computer ＆Graphics，2000，24（3）：353-361.

［10］唐加福，汪定偉，高振，等.面向非線性規(guī)劃問(wèn)題的混合式遺傳算法［J］.自動(dòng)化學(xué)報(bào)，2000，26（3）：401-404.

［11］KAJIYA J，VON H B.Ray tracing volume densities［J］.Computer Graphics，1984，18（3）：165-174.

［12］ELLIOT K，F(xiàn)ISHMAN，DEREK M D，NEY R.Volume Rendering versus Maximum Intensity Projection in CT Angiography：What Works Best，When，and Why［J］.RadioGraphics，2006，26（3）：905-922.

［13］LEVOY M.Efficient ray tracing of volume data［J］.ACM Trans.On Graphics，1990，9（3）：245-261.