柴柏林 張雅靜 李 琦
(1.中國石油大學(華東)計算機科學與技術學院 青島 266580)
(2.中國石油大學(華東)海洋與空間信息學院 青島 266580)
流場可視化技術是一種將矢量場數據映射為圖像或動畫,從而表現出其結構與特征的技術。海洋流場可視化通過對高度概括的海洋流場動態演變特征的分析提取與可視化表達,能夠深層次理解和歸納海洋流場的時空變化規律。因此在海洋應用領域使用越來越廣泛。陳長勝[1]等提出了有限體積沿海海洋模型的3D 非結構化三角網格(稱為FVCOM)。近年來,FVCOM 模型廣泛應用于海洋可視化、模擬流場、預測沿海地區特征[2~5]等領域,并實現了較好的效果。但目前世界上并沒有針對FVCOM 海洋模式的非結構化三角網格的可視化應用。
目前,關于種子點的放置問題,提出了隨機、均勻、圖像引導、流場引導、最遠距離種子點放置等一系列算法[6~9]。例如Zhan 和Robert[10]通過拓撲種子點的優先順序提出均勻分布的流線布置算法。Turk[11]等提出了一種圖像引導流線放置算法。Bhatia[12]等提出了在單純性內通過n+1個行列式快速判斷臨界點。但這些算法都不適用于FVCOM模式數據。點定位方面,杜小甫[13]提出了通過相鄰網格單元表進行點定位,但在數據結構的建立上存在存儲冗余。
在曲線積分過程中,由于固定步長精度較差,人們提出了自適應步長的計算方法。Buning[14]提出了步長與格網大小和速度矢量成倒數的方法,但計算較為復雜。又有部分學者通過計算相鄰流線向量之間的夾角來自適應步長[15~19]。季民等[20]提出了一個新的計算流線追蹤步長的模型(AMFCA算法),但此算法只針對于規則網格,沒有針對非結構化三角網格稠密漸變的結構特性而分析。王文濤[21]提出對于網格單元差異很大時采用局部步長,但并沒有展開討論。
綜上所示,本文針對FVCOM 模式非結構化三角網格的網格特性,通過添加局部網格尺寸和視角像素參數對AMFCA算法進行優化。同時結合粒子隨機化和點定位方法對粒子運動軌跡進行計算,并設計實現了基于非結構化三角網格的流場可視化應用。
如上述流程圖1 所示,本文首先對FVCOM 海洋模式數據預處理。在數據預處理中,分析數據中網格關系,構建非結構化三角網格的數據結構;求解每個三角網格的面積,根據三角網格面積和臨界點進行粒子隨機化。然后以隨機化的點為初始點,通過二階龍格庫塔法、改進的AMFCA 算法和點定位對粒子軌跡計算。

圖1 系統總體設計流程圖
本文基于Cesium[22~23]三維引擎進行粒子繪制和可視化應用的搭建。Cesium 是一個開源的Web地理空間數據可視化框架,它提供了一個效果精確、操作簡易的渲染引擎,同時支持三維地球、地圖、數據管理和用戶界面,相比其他框架較為成熟而且功能也更為豐富,在眾多領域都有廣泛的應用。
FVCOM 模式數據的非結構化三角網格如圖2所示。網格在海陸邊界處分布密集,在海洋內部分布稀疏。如果按照網格序號隨機選取網格,會導致海陸邊界處粒子數量過多、流場遮擋問題嚴重、視覺效果混亂,而在大洋內部粒子數量過少很容易丟失流場中表征重要特征的流線。

圖2 網格效果圖
基于上述分析,我們將粒子按照網格面積大小進行隨機化,算法流程如下所示:
步驟一:遍歷所有三角網格,計算每個三角網格的面積和三角網格面積總和。
步驟二:取最大和最小網格的面積差S,將[0,S]區間分成等大的n 份,將所有三角網格根據面積大小分成n個數組。
步驟三:計算n 個子數組中所有網格的面積之和,根據n個子數組間的面積和之比,按比例在n個數組中隨機選取粒子。
在隨機化的基礎上,為了保證粒子隨機化不丟失重要的網格特征點,我們通過Bhatia 提出的臨界點快速檢測方法提取海洋中的臨界點。Bhatia 證明了一個單純形內部是否包含臨界點可以通過它的n+1 個頂點a0,a1,…,an 構成的行列式所決定,如式(1)所示:
用原點替換第n 個頂點后,通過這n+1 個行列式符號是否相同來確定當前網格內是否存在臨界點。臨界點的檢測效果圖如圖3所示。

圖3 臨界點效果圖
本文將臨界點周圍區域與根據網格面積的隨機化相結合,最終形成粒子初始化集合,結果如圖4所示。

圖4 隨機化結果圖
目前流線數值積分的主要方法有歐拉法、二階龍格庫塔法、四階龍格庫塔法等[24]。歐拉法雖然計算量小,但是誤差較大。四階龍格庫塔法計算精度高,但計算也較為復雜。同時由于FVCOM 全球海洋數據集較為龐大,結合四階龍格庫塔法計算會大大降低可視化效率,故本文使用二階龍格庫塔方法,其公式如式(2)、(3)、(4)所示:
其中,k1、k2代表斜率,而h是所取的時間步長。
2014 年季民等提出了一種同時顧及網格流速和流向的AMFCA算法。該算法綜合考慮了當前網格的流速和流向兩個因素,達到了同向網格積分步長的快速增長和逆向網格積分步長的迅速銳減。具體模型如式(5)所示:
式中,Dij表示規則網格第i 行、第j 列網格的積分步長;D 表示流場中的數據網格尺寸;μ表示全局步長增長控制參數,主要用于步長增長的平衡控制;αij表示當前網格流速與流場中最大流速的比值,其取值范圍為[0,1];δ表示步長增長倍數,一般取大于1的值;θij表示流入與流出當前格網兩個流向的夾角,取值范圍為[0,π]。
AMFCA 算法在規則網格下流線追蹤具有較好的特性,但只適用于結構化網格。對于非結構化三角網格,若以相同的網格尺寸進行計算,在近岸密集的小網格上步長會過大,導致計算精度低。在大網格上步長過小,造成無謂增多的遞推計算。同時,基于動態粒子的流場可視化方法,根據可視化地球視角的放大縮小而改變粒子運動的快慢,會有更好的可視化效果。故本文在AMFCA算法基礎上兼顧網格大小、視角像素兩個因素,具體模型如式(6)所示:式中,K 代表網格尺寸。將三角網格根據面積分為小于0.03、0.03~0.2、大于0.2三個區間,若三角網格面積小時我們采用小步長,這樣可以使單次步長不會跨越過多的網格,保證可視化的準確性,面積較大時采用原始步長,K 根據網格面積尺寸取1/4、1/2、1。為了保證可視化效果的相對速度,當K 取1/4時,我們將4 次小步長的值相加得到最終步長。μ表示全局步長增長控制參數;ω表示視角像素值參數;α表示速度分量,δ表示流向局部控制參數,θ表示流入與流出網格的兩個流向的夾角,取值范圍為[0,π]。
在流線追蹤的過程中,模型綜合考慮了當前網格的大小、流向、流速、視角像素四個因素,并能通過μ、δ雙自由度調整。達到同向網格步長增加、逆向網格步長銳減;視角像素值大的步長略大、視角像素值小的步長略小;根據網格尺寸大小采用局部步長。這樣使可視化結果更精確,適應多種比例尺和曲率變化情況。
通過二階龍格庫塔法和自適應步長法可以計算出粒子運動終點。但由于非結構化三角網格的特性,無法通過某一經緯度坐標直接定位到對應的三角網格,因此需要建立邊列表數據結構,通過邊列表查詢粒子運動后所對應的網格。
杜小甫提出的點定位方法中,網格邊界會存儲大量-1 值,造成數據冗余。本文以“頂點序號-頂點序號”為鍵,代表兩個頂點間的邊。以“網格序號,網格序號”為值,代表邊所對應的1 或2 個三角形網格。查詢粒子坐標對應網格的算法步驟如下:
1)判斷粒子下一時刻位置是否在當前三角形內,若在則返回當前三角形序號,若不在則進入步驟2)。
2)設連接粒子當前位置和下一時刻位置的線段(圖5 中黑線所示),查詢該線段與當前三角形有交點的邊。

圖5 粒子計算原理
3)查詢該三角形的邊是否存在所對應的另一個三角形(一條邊最多可對應兩個三角形),若不存在則證明是網格邊界(海陸邊界),粒子消亡。存在則進入步驟4)。
4)把此三角形賦值給當前三角形。循環執行步驟1)。
本文實驗的PC 配置為操作系統為64bit 的Windows10,處 理 器 為 Intel(R) Core (TM)i7-4702MQ,內存為16GB,顯卡為NVIDIA GeForce GTX 760M,開發軟件為Microsoft VS Code,瀏覽器為Google Chrome。數據來自廣州國家超級計算中心提供的FVCOM 海洋模式數據,數據分全球網格和珠江口區域網格數據,包含了流場U 分量、流場V分量、垂向深度自海洋表面至海底分50層。
本文采用的全球網格數據時間為2019 年6 月21日0時,珠江口區域網格數據時間為2019年6月1 日0 時。垂向深度為海洋表面流場。在實驗中,我們通過粒子隨機化方法生成了256*256 個粒子。ω視角像素值在小于2000、2000~9000、大于9000 下分別取值1、2、3。取μ=0.5,δ=2。我們將改進后的AMFCA算法與AMFCA算法進行對比,保證其他參數一致的情況下,在小網格區域單次步長平均穿越網格數減少了2.2 倍,提高了可視化的準確性。可視化效果如圖6、圖7所示,通過實驗結果可以清晰的展示海洋中的流場特征和動態變化規律。

圖6 我國珠江口區域流場細節圖

圖7 大西洋區域流場可視化效果圖
本文通過基于網格面積的粒子隨機化方法能簡單、快速、均勻地覆蓋全球海洋區域及珠江口區域。在AMFCA算法的基礎上根據網格面積自適應步長大小,提高了可視化的準確性;同時添加視角像素參數,使改進后的算法能隨著可視化視角更好地描繪大型的流場軌跡和局部流場細節。基于邊列表的點定位方法避免了數據冗余,使得數據存儲更小。本文在Cesium 引擎上實現了非結構化三角網格海洋流場可視化應用,準確美觀地表現出海洋流場的真實軌跡,可視化效果較好。