基于局部體樣條函數的海洋鹽度場建模

楊振發，萬 剛，李 鋒，張宗佩

（信息工程大學，河南 鄭州450001）

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

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

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

1 海水鹽度數據特點與來源

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

圖1 海水鹽度場的分布規律

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

2 海水鹽度場的建模

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

圖2 鹽度場建模流程

2.1 局部的體樣條函數法

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

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

2）構造分塊區域Ωijk的權函數Wijk（P）。權函數滿足的條件如下：

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

文獻（3）采用Hermite的方法，給出具體求解權函數的方法。

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

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

2.2 非線性規劃的遺傳算法實現區域的優化分割

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

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

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

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

其中，b是區間［0，1］的隨機數。

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

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

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

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

圖3 非線性規劃的遺傳算法實現流程

3 實驗及結果

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

圖4 局部體樣條函數建模的實驗結果

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

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

圖5 最大密度投影算法實現流程

圖6 MIP算法實現的某海域鹽度場體視化結果

4 結 論

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

［1］朱振華.雙三次B樣條曲面法在聲吶測量數據空白填補中的作用［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］楊冠杰，耿震，孫菁.散亂體數據可視化：海洋水團分析的新途徑［J］.海洋通報，2000，19（2）：66-74.

［6］韓慶，李本江，高峰，等.海洋小百科全書［M］.中山：中山大學出版社，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］唐加福，汪定偉，高振，等.面向非線性規劃問題的混合式遺傳算法［J］.自動化學報，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，FISHMAN，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.