三維地質(zhì)模型海量數(shù)據(jù)組織和可視化技術(shù)研究

魏 嘉，唐 杰，武港山，張 楊

（1.中國石油化工股份有限公司石油物探技術(shù)研究院，江蘇南京211103；2.南京大學(xué)，江蘇南京210083）

復(fù)雜油氣藏三維地質(zhì)建模技術(shù)已成功應(yīng)用于油氣藏開發(fā)等方面，然而，隨著油氣藏勘探開發(fā)的不斷深入，對油氣藏地質(zhì)建模工作提出了更高的要求［1］。油氣藏地質(zhì)建模的工作流程分為構(gòu)造建模和屬性建模兩個(gè)階段，構(gòu)造建模是復(fù)雜油氣藏地質(zhì)建模的基礎(chǔ)，其技術(shù)的形成比屬性建模技術(shù)晚，但發(fā)展速度明顯快于屬性建模技術(shù)，尤其是隨著計(jì)算幾何和三維可視化等計(jì)算機(jī)技術(shù)的發(fā)展，促進(jìn)了三維構(gòu)造建模相關(guān)技術(shù)方法的進(jìn)步和發(fā)展［2］。早在1984年，Haldorson提出了油田尺度下用于油藏動(dòng)態(tài)模擬的隨機(jī)模擬儲(chǔ)層建模方法［3］；Mallet于1989年和1992年先后發(fā)表了兩篇關(guān)于離散光滑插值［4－5］建模方法的文章，標(biāo)志著三維構(gòu)造建模技術(shù)中的地質(zhì)曲面技術(shù)獲得了突破；為了不斷提高三維地質(zhì)模型的描述精度，層位面和斷層面三角網(wǎng)格重建、地質(zhì)曲面空間拓?fù)浞治龊偷刭|(zhì)塊體劃分等技術(shù)得到了長足的發(fā)展［6－11］。復(fù)雜油氣藏三維地質(zhì)建模精度的提高，必然會(huì)產(chǎn)生巨大的計(jì)算量和數(shù)據(jù)量，這對三維地質(zhì)模型的數(shù)據(jù)管理和三維可視化顯示提出了挑戰(zhàn)；有效地組織和管理精細(xì)地質(zhì)模型的海量數(shù)據(jù)、實(shí)時(shí)地進(jìn)行三維地質(zhì)模型海量數(shù)據(jù)的三維可視化顯示成為急需解決的技術(shù)難題。

三維可視化技術(shù)改變了利用三維地震資料進(jìn)行油氣藏表征和建模的工作流程，并成為降低油氣勘探風(fēng)險(xiǎn)和不確定性的有效技術(shù)手段［12］。由于受到客觀技術(shù)條件的制約（計(jì)算機(jī)內(nèi)存和圖形卡顯存總是不能滿足三維地震數(shù)據(jù)體數(shù)據(jù)量的要求），為了能夠?qū)崿F(xiàn)海量地震規(guī)則數(shù)據(jù)體的實(shí)時(shí)顯示，許多學(xué)者對在三維地震數(shù)據(jù)體繪制中如何減少I／O操作的問題進(jìn)行了研究［13－16］，主要方法包括減少體繪制過程中一次讀入的數(shù)據(jù)量和使用壓縮方法減少需要讀入的總數(shù)據(jù)量；另外，基于數(shù)據(jù)分塊和層次細(xì)節(jié)（LOD）模型技術(shù)進(jìn)行三維地震等規(guī)則數(shù)據(jù)體的多分辨率數(shù)據(jù)組織和管理技術(shù)也得到了廣泛的應(yīng)用［17－20］，通過基于擴(kuò)展四叉樹或八叉樹等數(shù)據(jù)結(jié)構(gòu)，建立分塊多層的多分辨率模型進(jìn)行LOD數(shù)據(jù)建模［21］，可以顯著提高海量地震數(shù)據(jù)體存取和顯示的效率。

三維地質(zhì)模型數(shù)據(jù)主要包括以空間多邊形（三角網(wǎng)格為主）表達(dá)的地質(zhì)層位面和斷層面（三維構(gòu)造框架模型）以及以空間多面體（任意六面體和四面體等）表達(dá)的地層網(wǎng)格（三維地層實(shí)體模型的基本單元）；三維地震數(shù)據(jù)體主要采用規(guī)則化立方體進(jìn)行數(shù)據(jù)的描述。對于非規(guī)則化的地質(zhì)曲面數(shù)據(jù)和地層網(wǎng)格數(shù)據(jù)的多分辨率數(shù)據(jù)組織和管理，有兩點(diǎn)基本要求：①不同分辨率的數(shù)據(jù)不能改變層位面、斷層面和地質(zhì)體的空間形態(tài)；②能夠?yàn)槿S地質(zhì)模型三維可視化實(shí)時(shí)顯示提供快速的數(shù)據(jù)存取。

1 基于外存的地質(zhì)曲面海量數(shù)據(jù)多分辨率組織和繪制

目前國內(nèi)外已提出的解決大規(guī)模三角網(wǎng)格實(shí)時(shí)顯示問題的方法可以分為層次細(xì)節(jié)模型和多分辨率模型兩類。層次細(xì)節(jié)模型方法是指預(yù)先在原模型數(shù)據(jù)的基礎(chǔ)上，生成若干不同分辨率的近似模型，在實(shí)時(shí)繪制中根據(jù)當(dāng)前的視點(diǎn)參數(shù)選用相應(yīng)的近似模型進(jìn)行繪制。這類算法的優(yōu)點(diǎn)是簡單易用、速度快；缺點(diǎn)是近似模型之間不連續(xù)，顯示時(shí)會(huì)有跳躍感，而層次的增多必然導(dǎo)致整個(gè)模型占用的存儲(chǔ)空間增大。多分辨率模型方法是將網(wǎng)格組織成一個(gè)基網(wǎng)格加上一系列細(xì)化操作?？梢暬@示時(shí)，隨著視點(diǎn)的移動(dòng)，多分辨率網(wǎng)格將實(shí)時(shí)自適應(yīng)更新。相對于層次細(xì)節(jié)模型來說，這種方式建立起來的場景更加符合人的視覺特性，模型的細(xì)節(jié)是連續(xù)變化的。但是，多分辨率網(wǎng)格的建立和自適應(yīng)更新要耗費(fèi)額外的時(shí)間，其技術(shù)難點(diǎn)是要在精度和速度之間選擇一個(gè)合適的平衡點(diǎn)。

Lindstrom等［22］在1996年首先提出了針對數(shù)字高程（DEM）模型的實(shí)時(shí)動(dòng)態(tài)多分辨率地形繪制算法，并不斷改進(jìn)了相關(guān)算法［23－24］。該算法基于一種四叉樹的數(shù)據(jù)結(jié)構(gòu)，將大范圍的DEM數(shù)據(jù)用四叉樹思想遞歸地分割成小的地形塊，再根據(jù)視覺范圍來控制地形塊的分裂和合并。該算法處理的三角網(wǎng)格是規(guī)則的DEM模型，需要在矩形參數(shù)域上等間距規(guī)則采樣。另外，Hoppe等［25－26］提出了針對非規(guī)則三角形網(wǎng)格（TIN）模型的漸進(jìn)網(wǎng)格算法，該算法首先將網(wǎng)格簡化，并記錄其簡化過程，進(jìn)而將簡化的數(shù)據(jù)組織成一個(gè)層次結(jié)構(gòu)，并定義了分裂或合并操作的有效性規(guī)則，從而將三角網(wǎng)格組織成一個(gè)基網(wǎng)格和一系列頂點(diǎn)關(guān)聯(lián)操作，方便了顯示不同的分辨率。

1.1 地質(zhì)曲面海量三角網(wǎng)格的分塊

由于對三維地質(zhì)模型中地質(zhì)曲面建模精度要求的不斷提高，在三維層位面和斷層面建模過程中，需要處理上百萬、甚至于上千萬的三角網(wǎng)格，而計(jì)算機(jī)的內(nèi)存總是有限的。我們對Hoppe等［25－26］提出的算法進(jìn)行了改進(jìn)，在創(chuàng)建地質(zhì)曲面海量三角網(wǎng)格多分辨率數(shù)據(jù)時(shí)，提出并實(shí)現(xiàn)了基于外存的海量三角網(wǎng)格多分辨率數(shù)據(jù)創(chuàng)建算法，該算法采用二叉樹的數(shù)據(jù)結(jié)構(gòu)對三角網(wǎng)格模型進(jìn)行劃分，建立層次結(jié)構(gòu)，依次對每塊數(shù)據(jù)進(jìn)行多分辨率模型創(chuàng)建，并將中間結(jié)果輸出到臨時(shí)性的數(shù)據(jù)文件，最后再反向創(chuàng)建多分辨率網(wǎng)格的數(shù)據(jù)文件。

數(shù)據(jù)分塊的方法有八叉樹劃分、四叉樹劃分和二叉樹劃分等。八叉樹劃分和四叉樹劃分等在用于劃分規(guī)則數(shù)據(jù)體（例如圖像和規(guī)則體數(shù)據(jù)）時(shí)非常有效，因?yàn)檫@些數(shù)據(jù)的拓?fù)湫畔?，即樣本之間的相互位置關(guān)系是隱含在樣本的索引值上的。這樣相鄰兩塊數(shù)據(jù)之間分辨率容易統(tǒng)一。而三角網(wǎng)格的拓?fù)湫畔⑹秋@式表達(dá)的，即三角網(wǎng)格不僅有其幾何信息（各頂點(diǎn)的空間位置），還要顯式地記錄拓?fù)湫畔ⅲǜ黜旤c(diǎn)之間的連接性信息）。各數(shù)據(jù)塊分別處理，相鄰兩塊數(shù)據(jù)之間的拓?fù)湫畔⒑茈y做到一致，輕則造成邊界處的數(shù)據(jù)裂縫，重則導(dǎo)致多分辨率模型崩潰。因此，為減少相鄰數(shù)據(jù)塊之間的同步壓力，在分析了三角網(wǎng)格二維流形的特征以后，我們選擇二叉樹作為分塊的數(shù)據(jù)結(jié)構(gòu)（圖1）。其基本算法流程為：①掃描三角網(wǎng)格文件，計(jì)算三角網(wǎng)格的頂點(diǎn)數(shù)、三角網(wǎng)格的個(gè)數(shù)以及三維構(gòu)造框架模型的包圍盒大小；②根據(jù)包圍盒的大小以及設(shè)定的每塊數(shù)據(jù)的最大頂點(diǎn)數(shù)，沿包圍盒最長方向計(jì)算二叉樹劃分順序和劃分層數(shù)；③在二叉樹劃分順序和劃分層數(shù)計(jì)算的基礎(chǔ)上，構(gòu)建二叉樹的層次父子結(jié)構(gòu)；④根據(jù)三角網(wǎng)格頂點(diǎn)的空間位置和三角網(wǎng)格邊簡化的相關(guān)性，將其歸屬到相應(yīng)層次的葉節(jié)點(diǎn)上；⑤建立三角網(wǎng)格頂點(diǎn)與二叉樹不同層次葉節(jié)點(diǎn)的對應(yīng)關(guān)系，并設(shè)置三角網(wǎng)格3個(gè)頂點(diǎn)的合并關(guān)系；⑥輸出三角網(wǎng)格數(shù)據(jù)分塊的劃分結(jié)果。

圖1 三角網(wǎng)格頂點(diǎn)的層次結(jié)構(gòu)

1.2 基于拓?fù)湟恢碌亩喾直媛誓Ｐ偷膭?chuàng)建

創(chuàng)建層位面或斷層面多分辨率網(wǎng)格模型時(shí)，首先由葉子節(jié)點(diǎn)開始，依次對每個(gè)節(jié)點(diǎn)的網(wǎng)格頂點(diǎn)數(shù)量進(jìn)行簡化。在簡化過程中，不僅需要保證層位面和斷層面的空間形態(tài)基本保持不變，同時(shí)也要保證相鄰節(jié)點(diǎn)中網(wǎng)格模型的拓?fù)湫畔⒁恢?。為此，我們提出了邊界三角網(wǎng)格整形處理方法，其中關(guān)鍵環(huán)節(jié)包括三角網(wǎng)格頂點(diǎn)簡化的誤差計(jì)算、三角網(wǎng)格形狀的評估因子計(jì)算和數(shù)據(jù)分塊的邊界處理。

1.2.1 三角網(wǎng)格頂點(diǎn)簡化的誤差計(jì)算

在簡化三角形頂點(diǎn)的過程中，首先在每個(gè)頂點(diǎn)計(jì)算一個(gè)誤差系數(shù)矩陣，當(dāng)該頂點(diǎn)被合并到另一個(gè)頂點(diǎn)，也就是簡化三角網(wǎng)格的一條邊時(shí)，該系數(shù)矩陣也被加到那個(gè)頂點(diǎn)的系數(shù)矩陣上。當(dāng)三角網(wǎng)格的一條邊簡化后，三角網(wǎng)格的形狀在該條邊的鄰接范圍內(nèi)發(fā)生了變化，其新頂點(diǎn)到原來的鄰接三角網(wǎng)格的距離的平方和即為簡化誤差，可以用公式（1）計(jì)算：

其中，v表示空間一點(diǎn)；EQ（v）表示空間一點(diǎn)v到某平面的距離誤差；D（v）表示空間一點(diǎn)到某個(gè)平面的距離，

A是一個(gè)3×3的系數(shù)矩陣；b為系數(shù)向量；c為系數(shù)標(biāo)量。

1.2.2 三角網(wǎng)格形狀的評估因子計(jì)算

在三角網(wǎng)格剖分計(jì)算過程中，出現(xiàn)內(nèi)角非常小的狹長三角形，對于后期的數(shù)據(jù)處理和圖形繪制都是不利的，現(xiàn)有的大部分三角網(wǎng)格計(jì)算方法要求三角網(wǎng)格應(yīng)盡量趨于正三角形。為了避免計(jì)算結(jié)果出現(xiàn)不可接受的誤差，我們引入Gueziec［27］提出的三角網(wǎng)格形狀的評估因子：

式中：li為邊長；ω為三角片的面積。

γ值為0～1。γ越小，說明三角網(wǎng)格的狹長度越大，越不利于三角網(wǎng)格的處理；γ越大，說明三角片越接近于等邊三角形。為防止狹長三角片的出現(xiàn)，再計(jì)算刪除邊后所有目標(biāo)點(diǎn)的一階鄰接三角片形狀評估因子，如果存在小于指定閾值的值，表示有三角網(wǎng)格狹長，則在其計(jì)算出的簡化誤差上加上一個(gè)值，推后其在網(wǎng)格簡化隊(duì)列中的位置。

1.2.3 三角網(wǎng)格數(shù)據(jù)分塊的邊界處理

對葉節(jié)點(diǎn)處理完畢后，根據(jù)二叉樹層次結(jié)構(gòu)，依次將相鄰葉節(jié)點(diǎn)剩余的頂點(diǎn)和三角網(wǎng)格以及相鄰處的三角網(wǎng)格合并起來（圖2，紅圈處表示分塊處理的塊內(nèi)三角網(wǎng)格，需要進(jìn)行合并處理），再進(jìn)行簡化處理。

圖2 在分塊邊界上合并相鄰的節(jié)點(diǎn)

通過計(jì)算誤差系數(shù)矩陣，保證了地質(zhì)曲面空間形態(tài)的相似性；通過計(jì)算三角網(wǎng)格形態(tài)評價(jià)因子，避免了狹長三角形的出現(xiàn)，可以提高數(shù)據(jù)處理的穩(wěn)定性；通過綜合應(yīng)用誤差系數(shù)矩陣計(jì)算、三角網(wǎng)格形態(tài)評價(jià)因子計(jì)算和邊界三角網(wǎng)格的特殊處理方法，實(shí)現(xiàn)了地質(zhì)曲面多分辨率數(shù)據(jù)組織的幾何一致和拓?fù)湟恢拢▓D3：圖2是分塊處理的地質(zhì)曲面形態(tài)，與實(shí)際地質(zhì)曲面的拓?fù)潢P(guān)系存在一定的差異；進(jìn)行三角網(wǎng)格的整形優(yōu)化處理后得到圖3，保證了地質(zhì)曲面的幾何一致和拓?fù)湟恢拢?/p>

1.3 大規(guī)模三角網(wǎng)格多分辨率模型的實(shí)時(shí)繪制

由于地質(zhì)曲面多分辨率模型中的點(diǎn)、三角網(wǎng)格以及網(wǎng)格頂點(diǎn)的分裂和合并操作都是以連續(xù)的方式存放的，因此可以利用內(nèi)存映射機(jī)制將存在磁盤上的數(shù)據(jù)映射進(jìn)內(nèi)存，借助操作系統(tǒng)的頁面換入／換出機(jī)制實(shí)現(xiàn)數(shù)據(jù)調(diào)度。當(dāng)多分辨率數(shù)據(jù)量較大時(shí)，不能一次性映射進(jìn)內(nèi)存，需要分批將磁盤文件映射進(jìn)內(nèi)存。

初始時(shí)，前端頂點(diǎn)保存基網(wǎng)格頂點(diǎn)。每一幀中，根據(jù)視點(diǎn)參數(shù)判斷活躍頂點(diǎn)是否需要合并或者分裂。當(dāng)活躍頂點(diǎn)在可視范圍內(nèi)且誤差大于指定閾值時(shí)拆分該點(diǎn)，當(dāng)前端頂點(diǎn)不在可視范圍或父節(jié)點(diǎn)誤差小于閾值時(shí)合并該點(diǎn)。最后采用映射緩沖區(qū)對象的方法來加速繪制速度。緩沖區(qū)對象是一種數(shù)據(jù)存儲(chǔ)倉庫，用于在本地視頻內(nèi)存中存儲(chǔ)數(shù)據(jù)。在頂點(diǎn)數(shù)組繪制模式下，將活躍三角形直接生成到顯存可以顯著提升繪制速率，繪制結(jié)果如圖4所示。圖4中紅框表示實(shí)際的可視區(qū)域，該區(qū)域內(nèi)采用高分辨率的三角網(wǎng)格數(shù)據(jù)進(jìn)行顯示，該區(qū)域以外的部分采用較低分辨率的三角網(wǎng)格數(shù)據(jù)。通過應(yīng)用上述技術(shù)基本達(dá)到了海量三角網(wǎng)格多分辨率實(shí)時(shí)繪制的要求（圖5）。

2 海量三維地層網(wǎng)格數(shù)據(jù)的多分辨率數(shù)據(jù)組織和可視化

地質(zhì)曲面組合構(gòu)建的是三維地層框架模型，但是，僅僅有層位面和斷層面，對于三維地質(zhì)建模是不夠的，還需要對地層的內(nèi)部實(shí)體進(jìn)行建模。地層實(shí)體建模的方法是基于已經(jīng)建立的地層框架模型，對模型中的地層按照沉積約束條件進(jìn)行地層網(wǎng)格剖分，進(jìn)而建立以三維地層網(wǎng)格為單元的地層實(shí)體模型［2］。由于復(fù)雜油氣藏的三維地質(zhì)建模的精度要求越來越高，三維地層網(wǎng)格的剖分粒度越來越細(xì)，一個(gè)三維地質(zhì)模型的地層網(wǎng)格數(shù)據(jù)量往往達(dá)到了數(shù)百GB級別，目前的三維可視化顯示技術(shù)很難滿足實(shí)時(shí)性的要求。為此，我們提出了一種基于數(shù)據(jù)分塊的層次細(xì)節(jié)模型的三維地層網(wǎng)格數(shù)據(jù)組織和管理技術(shù)。

2.1 三維地層網(wǎng)格數(shù)據(jù)的數(shù)據(jù)分塊

空間數(shù)據(jù)分塊是一種基于外存數(shù)據(jù)調(diào)度策略的技術(shù)，該技術(shù)充分利用地質(zhì)模型的空間分布性質(zhì)將網(wǎng)格數(shù)據(jù)在空間上劃分為若干子塊，每個(gè)子塊可獨(dú)立進(jìn)行處理。這樣對快速地進(jìn)行視錐判斷和圖形裁減非常有利，而且可以采用多線程的方法對多個(gè)數(shù)據(jù)子塊并行化地進(jìn)行處理，從而加快數(shù)據(jù)處理的速度，減少響應(yīng)等待時(shí)間，改善三維地質(zhì)建模過程的交互性能。

在基于數(shù)據(jù)分塊的調(diào)度系統(tǒng)中，數(shù)據(jù)塊的大小是直接影響系統(tǒng)效率的重要因素，過小的分塊導(dǎo)致過多的I／O中斷，降低磁盤的讀取效率；過大的分塊不僅會(huì)導(dǎo)致讀入過多的無用數(shù)據(jù)，而且較長的單次讀取時(shí)間會(huì)造成整個(gè)系統(tǒng)的響應(yīng)緩慢。因此，必須慎重選擇數(shù)據(jù)分塊的大小?？紤]到操作系統(tǒng)進(jìn)行頁面調(diào)度和磁盤緩存的因素，為減少缺頁中斷的次數(shù)，我們將數(shù)據(jù)分塊設(shè)置為32×32×32，這樣每個(gè)數(shù)據(jù)塊大小為8B×32×103＝256KB，并根據(jù)分塊尺度的信息建立相關(guān)模型數(shù)據(jù)的索引，從而提高數(shù)據(jù)的存取效率。

2.2 面向地層網(wǎng)格數(shù)據(jù)的LOD模型

傳統(tǒng)的三維可視化體繪制方法并不區(qū)分場景中對象的狀態(tài)，所有位于視錐范圍之內(nèi)的對象都會(huì)被以最高精度進(jìn)行繪制，一旦需要繪制的對象數(shù)據(jù)量急劇增大，這種繪制全部細(xì)節(jié)的方法便會(huì)力不從心。事實(shí)上，由于屏幕分辨率有限，無限制地提高繪制精度對于最終顯示效果的提高作用微乎其微。當(dāng)對象在場景中距離視平面較遠(yuǎn)時(shí)，可以采用較低分辨率的數(shù)據(jù)進(jìn)行繪制，而并不會(huì)明顯降低顯示的質(zhì)量。在處理大規(guī)模地層網(wǎng)格數(shù)據(jù)方面，我們采用了LOD模型進(jìn)行數(shù)據(jù)的組織和管理。

目前面向大規(guī)模數(shù)據(jù)實(shí)時(shí)繪制的LOD技術(shù)主要分為靜態(tài)LOD方法和動(dòng)態(tài)LOD方法［28］。靜態(tài)LOD方法以離線方式預(yù)先對數(shù)據(jù)模型進(jìn)行處理，生成一系列不同細(xì)節(jié)層次的網(wǎng)格模型，繪制時(shí)根據(jù)顯示參數(shù)和硬件性能選擇一個(gè)合適的層次進(jìn)行繪制。靜態(tài)LOD方法將耗費(fèi)資源的數(shù)據(jù)生成過程作為預(yù)處理，因此實(shí)時(shí)性較好，但在視點(diǎn)實(shí)時(shí)變換及細(xì)節(jié)層次突變的處理上不如動(dòng)態(tài)LOD方法。動(dòng)態(tài)LOD方法則是在繪制時(shí)才生成相應(yīng)的細(xì)節(jié)層次進(jìn)行渲染，因此需要對模型作實(shí)時(shí)簡化計(jì)算，對硬件性能要求較高［29］。考慮到大規(guī)模數(shù)據(jù)的實(shí)時(shí)處理效率，我們?yōu)楹Ａ康貙泳W(wǎng)格數(shù)據(jù)設(shè)計(jì)了一種靜態(tài)LOD模型。

八叉樹結(jié)構(gòu)是由四叉樹結(jié)構(gòu)推廣到三維空間而形成的一種三維數(shù)據(jù)結(jié)構(gòu)［30－31］，它能夠?qū)θS空間形成有效劃分和高效的分層次管理，因此適合于組織空間數(shù)據(jù)的LOD模型。空間八叉樹通常采用自頂向下的方法構(gòu)建，節(jié)點(diǎn)之間采用指針的方式表達(dá)派生關(guān)系，優(yōu)點(diǎn)是查詢效率較高，缺點(diǎn)是大量的指針數(shù)據(jù)占據(jù)了本不充裕的存儲(chǔ)空間。我們針對地層網(wǎng)格存儲(chǔ)方式具有的規(guī)則數(shù)據(jù)體的特殊性，采用自下向上的順序構(gòu)建地層網(wǎng)格的八叉樹LOD模型。不同于傳統(tǒng)的八叉樹結(jié)構(gòu)，節(jié)點(diǎn)之間的指向關(guān)系不需要用指針進(jìn)行存儲(chǔ)，而是通過節(jié)點(diǎn)索引得到，節(jié)省了存儲(chǔ)大量指針數(shù)據(jù)的空間，更有利于海量數(shù)據(jù)模型的多分辨率處理。

在構(gòu)建層次細(xì)節(jié)模型時(shí)，首先依據(jù)數(shù)據(jù)分塊的大小對原始地層的規(guī)則數(shù)據(jù)體在X，Y和Z 3個(gè)維度上進(jìn)行空間劃分，形成最高級分辨率的數(shù)據(jù)塊陣列，并作為葉節(jié)點(diǎn)存入八叉樹的最低層。然后根據(jù)八叉樹的特點(diǎn)，在X，Y和Z 3個(gè)維度上分別對緊鄰的數(shù)據(jù)塊進(jìn)行二合一操作，這樣空間關(guān)系上緊鄰的8個(gè)數(shù)據(jù)塊將被合并為一個(gè)數(shù)據(jù)塊，形成次高級分辨率的數(shù)據(jù)塊陣列，并存入八叉樹的上一層。如果在某個(gè)維度上的數(shù)據(jù)塊個(gè)數(shù)為奇數(shù)，則增補(bǔ)一個(gè)空白數(shù)據(jù)塊進(jìn)行填充，以保證八叉樹結(jié)構(gòu)的完整性。這個(gè)過程順著八叉樹以自底向上的順序迭代進(jìn)行，直到八叉樹某個(gè)層次上的數(shù)據(jù)塊個(gè)數(shù)達(dá)到閾值為止。

2.3 地層網(wǎng)格LOD模型的數(shù)據(jù)流水線

地層網(wǎng)格多分辨率LOD模型的構(gòu)建過程中，控制點(diǎn)的生成需要大量的浮點(diǎn)計(jì)算，CPU資源的開銷很大；將數(shù)據(jù)塊寫入到輸出文件則需要與慢速的磁盤I／O系統(tǒng)協(xié)同工作，經(jīng)常造成CPU等待，資源閑置。這兩個(gè)階段在整個(gè)數(shù)據(jù)流程中耗時(shí)最長，同時(shí)系統(tǒng)資源卻沒有得到充分利用。為了充分發(fā)揮硬件性能，提高程序的運(yùn)行效率，我們將LOD模型的構(gòu)建過程進(jìn)行4級流水線處理，使用4個(gè)線程分別負(fù)責(zé)流水線中的一級（圖6）。

圖6 地層網(wǎng)格數(shù)據(jù)控制流圖示

2.4 基于外存的海量地層數(shù)據(jù)調(diào)度與可視化技術(shù)

目前基于外存的數(shù)據(jù)調(diào)度和繪制方面的研究主要集中在地形數(shù)據(jù)的顯示方面，特別是數(shù)字高程模型的三維可視化顯示［32－34］，而針對三維地質(zhì)模型的地層網(wǎng)格數(shù)據(jù)的研究則相對較少。我們設(shè)計(jì)了海量地層數(shù)據(jù)的調(diào)度與可視化框架（圖7），以外部磁盤上的海量數(shù)據(jù)文件組織與管理為基礎(chǔ)，兼顧磁盤讀寫效率與用戶交互性能，實(shí)現(xiàn)了層狀實(shí)體模型的三維可視化顯示和線、道、折線方向的剖面數(shù)據(jù)快速提取與顯示。

圖7 基于外存的數(shù)據(jù)調(diào)度與可視化框架

2.4.1 多分辨率數(shù)據(jù)調(diào)度策略

海量的地層網(wǎng)格數(shù)據(jù)經(jīng)過多分辨率組織與存儲(chǔ)技術(shù)處理后，形成了不同精細(xì)程度的幾個(gè)層次細(xì)節(jié)模型，每一層次的模型又在空間上以數(shù)據(jù)分塊的形式進(jìn)行組織并存儲(chǔ)。繪制時(shí)可根據(jù)當(dāng)前的顯示參數(shù)和所用機(jī)器的硬件性能選擇一個(gè)合適的細(xì)節(jié)層次，同時(shí)也不需要將該層次上所有數(shù)據(jù)塊全部調(diào)入主存，只需根據(jù)當(dāng)前的顯示參數(shù)實(shí)時(shí)調(diào)入正在或?qū)⒁L制的部分?jǐn)?shù)據(jù)塊。我們采用了一種基于優(yōu)先級隊(duì)列的自頂向下的多分辨率數(shù)據(jù)調(diào)度策略：優(yōu)先級隊(duì)列采用最大堆實(shí)現(xiàn)，存放當(dāng)前進(jìn)入調(diào)度系統(tǒng)的數(shù)據(jù)塊索引，按照優(yōu)先權(quán)函數(shù)值進(jìn)行排列，越靠近隊(duì)列頭的數(shù)據(jù)塊就越應(yīng)該在調(diào)度的下次操作中被細(xì)化。每當(dāng)有新的數(shù)據(jù)塊索引從優(yōu)先級隊(duì)的尾部入隊(duì)，或具有最大優(yōu)先權(quán)的數(shù)據(jù)塊索引從優(yōu)先級隊(duì)的頭部出隊(duì)，優(yōu)先級隊(duì)列都要進(jìn)行自上向下或自下向上的堆調(diào)整以保證正確的優(yōu)先權(quán)排序。

2.4.2 緩存管理方案

基于外存的數(shù)據(jù)組織技術(shù)要求數(shù)據(jù)調(diào)度模塊不斷地訪問外存設(shè)備獲取新的數(shù)據(jù)來更新當(dāng)前顯示的舊數(shù)據(jù)。從用戶瀏覽地質(zhì)體數(shù)據(jù)的方式來看，當(dāng)前可見的數(shù)據(jù)塊在不久的將來依然可見的可能性是相當(dāng)高的（幀－幀相關(guān)性），也就是說對外存數(shù)據(jù)的訪問具有時(shí)間局部性。因此，為了減少磁盤I／O交互次數(shù)，有效利用之前已經(jīng)讀入內(nèi)存的數(shù)據(jù)，我們設(shè)計(jì)了一個(gè)針對地層LOD模型的緩存管理方案，將內(nèi)、外存數(shù)據(jù)交互這個(gè)瓶頸對整個(gè)多分辨率數(shù)據(jù)調(diào)度與可視化系統(tǒng)性能帶來的不利影響降低到最小。

2.4.3 利用頂點(diǎn)緩沖的大規(guī)模繪制

經(jīng)過調(diào)度系統(tǒng)輸出的數(shù)據(jù)塊隨后進(jìn)入裝配與繪制過程。我們將需要繪制的數(shù)據(jù)塊中的控制點(diǎn)在空間上進(jìn)行拓?fù)潢P(guān)系檢查，以四邊形網(wǎng)格的形式作為線、道、折線方向提取的剖面及其屬性。為了更有效地將網(wǎng)格數(shù)據(jù)提交到渲染引擎，我們將剖面數(shù)據(jù)組織成頂點(diǎn)、索引、屬性3個(gè)數(shù)據(jù)流的形式，利用OpenGL提供的頂點(diǎn)數(shù)組機(jī)制和頂點(diǎn)緩沖對象的擴(kuò)展技術(shù)，減少OpenGL函數(shù)調(diào)用的次數(shù)，提高網(wǎng)格數(shù)據(jù)的傳輸效率，減少了每幀渲染時(shí)間，提高了幀率，保證數(shù)據(jù)顯示的實(shí)時(shí)性（圖8）。

圖8是對一個(gè)20GB數(shù)據(jù)量的三維地層網(wǎng)格數(shù)據(jù)進(jìn)行三維可視化測試的結(jié)果，從圖8a到圖8d的變化可以看出：在不適用緩存的情況下，幀率極不穩(wěn)定，經(jīng)常出現(xiàn)嚴(yán)重的掉幀，顯示過程中頻繁出現(xiàn)畫面停滯，嚴(yán)重影響了軟件交互的性能。隨著緩沖區(qū)中緩存數(shù)據(jù)塊的增加，幀率的穩(wěn)定性有較大幅度提高，當(dāng)緩存數(shù)據(jù)塊達(dá)到1 000個(gè)以上，雖仍有部分輕微的掉幀情況出現(xiàn)，但數(shù)據(jù)繪制的過程自然流暢，用戶交互操作反映較為靈敏，基本滿足了三維地質(zhì)模型的實(shí)時(shí)可視化要求。

圖8 不同緩存方案的幀率曲線

3 實(shí)際數(shù)據(jù)的算法驗(yàn)證

利用本文提出的方法對某實(shí)際地震構(gòu)造解釋數(shù)據(jù)進(jìn)行了層位面和斷層面的建模，并在地層框架模型的基礎(chǔ)上建立了地層實(shí)體模型，通過三維構(gòu)造建模的全過程對以上所述相關(guān)算法進(jìn)行了系統(tǒng)驗(yàn)證。

試驗(yàn)工區(qū)的三維地質(zhì)建模包括了3個(gè)層位面和若干個(gè)斷層面，每個(gè)層面的離散數(shù)據(jù)點(diǎn)數(shù)為130×104～150×104。圖9給出了其中一個(gè)層位面（共有133×104個(gè)數(shù)據(jù)點(diǎn)）的解釋數(shù)據(jù)在不同分辨率時(shí)的建模效果。從圖9a到圖9d可見，不同分辨率的三角網(wǎng)格數(shù)據(jù)空間形態(tài)和拓?fù)潢P(guān)系一致，并且在硬件條件相同的情況下，計(jì)算效率優(yōu)于某商業(yè)軟件得到的結(jié)果。

圖9 相同層位不同分辨率三維顯示結(jié)果

圖10給出了試驗(yàn)工區(qū)地層實(shí)體模型不同分辨率數(shù)據(jù)三維可視化顯示結(jié)果。該模型地層網(wǎng)格的精度為0.29m×0.29m；地層網(wǎng)格的層間距為0.29m；數(shù)據(jù)塊的數(shù)量為84 105塊；原始數(shù)據(jù)量為22 047.6MB。為了三維可視化顯示的流暢性，建立了7級LOD多分辨率數(shù)據(jù)，多分辨率數(shù)據(jù)的實(shí)際數(shù)據(jù)量為42.02GB，大約是原始數(shù)據(jù)量的1.91倍，生成7級LOD模型共耗時(shí)1 440s。

從圖10中可以看出，本文提出的地層網(wǎng)格多分辨率組織技術(shù)，不僅保證了地層網(wǎng)格組成的地質(zhì)體在不同分辨率的情況下幾何形態(tài)和拓?fù)潢P(guān)系一致，同樣，也保證了地層網(wǎng)格的屬性（圖10中的顏色顯示）在不同分辨率數(shù)據(jù)中，空間變化規(guī)律一致。本文提出的地層網(wǎng)格多分辨率數(shù)據(jù)組織方法不僅能夠支撐三維可視化顯示，而且也為三維地質(zhì)模型的交互操作提供了便利（圖11）。

4 結(jié)束語

1）針對海量地質(zhì)三角網(wǎng)格曲面的實(shí)時(shí)顯示問題，本文提出了一種基于二次誤差和邊簡化操作的三角網(wǎng)格多分辨率模型的數(shù)據(jù)組織及其創(chuàng)建算法。該算法首先計(jì)算模型中所有邊的簡化代價(jià)，并從小到大進(jìn)行排序；然后依次簡化模型中三角網(wǎng)格的邊，進(jìn)行向下采樣；最后在恢復(fù)分辨率的同時(shí)，創(chuàng)建該三角網(wǎng)格的多分辨率模型。與其它算法相比，本文提出的算法可以處理海量三角網(wǎng)格，不受內(nèi)存等系統(tǒng)資源的限制。

2）針對海量地層網(wǎng)格數(shù)據(jù)的特點(diǎn)，本文提出了一種基于空間八叉樹結(jié)構(gòu)分塊處理非規(guī)則化地層實(shí)體模型的多分辨率數(shù)據(jù)組織方法，該方法綜合運(yùn)用數(shù)據(jù)分塊和層次細(xì)節(jié)模型等技術(shù)手段，為海量地層網(wǎng)格數(shù)據(jù)生成適合于多分辨率應(yīng)用的LOD模型，并以此建立了地層網(wǎng)格數(shù)據(jù)的多層級索引，實(shí)現(xiàn)了海量地層網(wǎng)格數(shù)據(jù)的高效存取。

3）以海量地層網(wǎng)格數(shù)據(jù)組織與管理為基礎(chǔ)，本文提出了一種針對海量地層網(wǎng)格數(shù)據(jù)的調(diào)度與可視化框架，充分運(yùn)用基于外存技術(shù)的各種手段，實(shí)現(xiàn)了包含地層屬性的層狀實(shí)體模型的線、道、折線方向的剖面數(shù)據(jù)快速提取與顯示，同時(shí)有效提升了大規(guī)模數(shù)據(jù)調(diào)度時(shí)的磁盤讀寫效率，改善了大規(guī)模數(shù)據(jù)繪制時(shí)的交互操作性能。

4）三維地質(zhì)模型海量數(shù)據(jù)的組織管理和實(shí)時(shí)三維可視化是三維地質(zhì)建模的關(guān)鍵技術(shù)，本文從用戶體驗(yàn)出發(fā)，解決了海量三角網(wǎng)格數(shù)據(jù)、海量地層網(wǎng)格數(shù)據(jù)等三維地質(zhì)模型基本數(shù)據(jù)的多分辨率數(shù)據(jù)組織與三維可視化實(shí)時(shí)顯示的關(guān)鍵問題。上述技術(shù)與其它三維建模技術(shù)（地質(zhì)曲面重建、地質(zhì)曲面拓?fù)浞治?、地層網(wǎng)格剖分和相控屬性建模等技術(shù)）是研發(fā)面向復(fù)雜油氣藏的三維地質(zhì)建模軟件的技術(shù)基礎(chǔ)，這些技術(shù)的發(fā)展將進(jìn)一步促進(jìn)地震綜合解釋和地質(zhì)建模的一體化。

［1］魏嘉.地質(zhì)建模技術(shù)［J］.勘探地球物理進(jìn)展，2007，30（1）：1－6 Wei J.Review of geologic model building techniques［J］.Progress in Exploration Geophysics，2007，30（1）：1－6

［2］魏嘉，唐杰，岳承琪，等.三維地質(zhì)構(gòu)造建模技術(shù)研究［J］.石油物探，2008，47（4）：319－327 Wei J，Tang J，Yue C Q，et al.Study of 3－D geological structure model building［J］.Geophysical Prospecting for Petroleum，2008，47（4）：319－327

［3］Haldorson H，Lake L.A new approach to shale management in filed scale simulation models［J］.The Journal of Petroleum Technology，1984，24（4）：1－4

［4］Mallet J L.Discrete smooth interpolation［J］.ACM Transactions on Graphics，1989，8（2）：121－144

［5］Mallet J L.Discrete smooth interpolation in geometric modeling［J］.Computer Aided Design，1992，24（4）：177－191

［6］Euler N.A new tool to seal 3Dearth model：a cut with constraint［J］.Expanded Abstracts of 68thAnnual Internat SEG Mtg，1998，321－324

［7］楊欽，徐永安，陳其明，等.三維約束Delaunay三角化的研究［J］.計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào)，2000，12（8）：590－594 Yang Q，Xu Y A，Chen Q M，et al.Research on 3D constrained delaunay triangulation［J］.Journal of Computer Aided Design and Computer Graphics，2000，12（8）：590－594

［8］Thore P，Shtuka A，Lecour M，et al.Structural uncertainties：determination，management，and applications［J］.Geophysics，2002，67（3）：840－852

［9］王秀閩，劉洪，孟小紅，等.用于塊體劃分的最小回路法及模型試驗(yàn)［J］.石油物探，2007，46（3）：249－253 Wang X M，Liu H，Meng X H，et al.Model test and study on minimum circuits for block partition［J］.Geophysical Prospecting for Petroleum，2007，46（3）：249－253

［10］Lorensen W E.Marching cubes：a high resolution 3D surface construction algorithm［J］.Computer Graphics，1987，21（4）：163－169

［11］袁國棟，施寅，勵(lì)曉建.一種新的抽取等值面的四面體分解方法［J］.計(jì)算機(jī)工程與科學(xué)，2002，20（6）：48－51 Yuan G D，Shi Y，Li X J.A new tetrahedral decomposition scheme for extracting isosurfaces［J］.Computer Engineering ＆ Science，2002，20（6）：48－51

［12］Chambers H，Brown A L.3－D visualization continues to advance integrated interpretation environment［J］.First Break，2003，21（5）：31－34

［13］Hodgson P，Cunnell C，Krueger A，et al.Cluster visualization rises to the seismic processing challenge［J］.First Break，2005，23（3）：33－35

［14］Gray G.Observations of azimuthal anisotropy in prestack seismic data［J］.Expanded Abstracts of CSPG ＆ CSEG Convention，2007，373－377

［15］Molnar S，Cox M，Ellsworth D，et al.A sorting classification of parallel rendering［J］.IEEE Computer Graphics and Applications，1994，14（4）：23－32

［16］Varadhan G，Manocha D.Out－of－core rendering of massive geometric environments［J］.Expanded Abstracts of Proceedings of IEEE Visualization Mtg，2002，69－76

［17］Imma B，Isabel N，Roberto S.Multiresolution volume visualization with a texture－based octree［J］.The Visual Computer，2001，17：185－197

［18］Guthe S，Wand M，Gonser J，et al.Interactive rendering of large volume data sets［J］.Expanded Abstracts of Proceedings of IEEE Visualization Mtg，2002，104－115

［19］Schneider J，Westerm R.Compression domain volume rendering［J］.Expanded Abstracts of Proceedings of the IEEE Visualization Mtg，2003，293－300

［20］魏嘉，唐杰，武港山，等.地震疊前數(shù)據(jù)三維可視化技術(shù)探討［J］.勘探地球物理進(jìn)展，2009，32（1）：12－17 Wei J，Tang J，Wu G，et al.On 3－D visualization of prestack seismic data［J］.Progress in Exploration Geophysics，2009，32（1）：12－17

［21］魏嘉，唐杰，武港山，等.三維地震數(shù)據(jù)體多分辨率數(shù)據(jù)組織與管理技術(shù)研究［J］.石油物探，2010，49（3）：240－244 Wei J，Tang J，Wu G S，et al.Multi－resolution organization and management of massive 3Dseismic data volume［J］.Progress in Exploration Geophysics，2010，49（3）：240－244

［22］Lindstrom P，Koller D，Ribarsky W，et a1.Real－time continuous level of detail rendering of height fields［J］.Computer Graphics，1996，30（1）：109－118

［23］Lindstrom P，Pascucci V.Terrain simplification sim－plified：ageneral framework for view－dependent outof－core visualization［J］.IEEE Transactions on Visualization and Computer Graphics，2002，8（3）：239－254

［24］Lindstrom P，Isenburg M.Lossless compression of hexahedral meshes［C］∥Expanded Abstracts of IEEE Data Compression Conference.New York：IEEE Publishing，2008：192－201

［25］Hoppe H.Progressive meshes［C］∥Expanded Abstracts of Computer Graphics Mtg.New York：Association for Computing Machinery，1996：99－108

［26］Hoppe H.View－dependent refinement of progressive meshes［C］∥Expanded Abstracts of Computer Graphics Mtg.New York：Association for Computing Machinery，1997：189－198

［27］Gueziec A.Surface simplification with variable tolerance［C］∥Expanded Abstracts of 2ndAnnual Internat Symposiumon Medical Robotics and Computer Assisted Surgery Mtg.Wiley：Orthopaedic Research Society，1995：132－139

［28］Duchaineauy M，Wolinsky M.Roaming terrain：realtime optimally adapting meshs［C］∥Expanded Abstracts of IEEE Visualization Mtg.New York：IEEE Publishing，1997：81－88

［29］Rottger S，Heidrich W，Slussallek P，et al.Rea1－time generation of continuous levels of detail for height fields［J］.Expanded Abstracts of Proceedings of the 6thInternational Conference in Central Europe on Computer Graphics and Visualization，1998，315－322

［30］Samet H.The design and analysis of spatial data structures［M］.New Jersey：Addison－Wesley，1990：223－235

［31］Boada I.An octree－based multiresolution hybrid framework［J］.Future Generation Computer Systems，2004，20（8）：1275－1284

［32］Vitter J S.External memory algorithms and data structures：dealing with massive data［J］.ACM Computing Surveys，2001，33（2）：209－271

［33］Ulrich T.Chunked LOD：rendering massive terrains using chunked level of detail control［C］∥Computer Graphics Proceedings Annual Conference Series Expanded Abstracts.Texas，USA：ACM，2002：202－209

［34］Vo H T，Callahan S P，Smith N，et al.iRun：Interactive rendering of large unstructured grids［C］∥Eurographics Symposium on Parallel Graphics and Visualization Expanded Abstracts.Geneva：The Eurographics Association，2007：93－100