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

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

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

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

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

1 海水鹽度數(shù)據(jù)特點與來源

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

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

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

2 海水鹽度場的建模

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

圖2 鹽度場建模流程

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

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

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

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

a.當P?Ωijk時，Wijk（P）＝0；

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

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

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

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

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

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

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

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

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

5）變異操作算子。其目的是維護種群多樣性，糾正遺傳算法過早收斂的缺點。從種群中隨機選擇第i個體，第j個編碼編譯操作為

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

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

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

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

3 實驗及結(jié)果

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

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

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

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

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

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

4 結(jié) 論

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

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

［2］周玉娟，岳桂昌.二次B樣條修正VTEC多項式模型研究［J］.測繪工程，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ù)可視化：海洋水團分析的新途徑［J］.海洋通報，2000，19（2）：66-74.

［6］韓慶，李本江，高峰，等.海洋小百科全書［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ī)劃問題的混合式遺傳算法［J］.自動化學(xué)報，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.