張茂亮 劉真 陳德峰 郭正府＊＊ 郭文峰

ZHANG MaoLiang1，LIU Zhen2，CHEN DeFeng3，GUO ZhengFu1＊＊ and GUO WenFeng1

1. 中國科學(xué)院地質(zhì)與地球物理研究所新生代地質(zhì)與環(huán)境重點(diǎn)實(shí)驗(yàn)室，北京 100029

2. 中國人民大學(xué)信息學(xué)院，北京 100872

3. 首都師范大學(xué)檢測成像實(shí)驗(yàn)室，北京 100048

1. Key Laboratory of Cenozoic Geology and Environment，Institute of Geology and Geophysics，Chinese Academy of Sciences，Beijing 100029，China

2. Information School，Renmin University of China，Beijing 100872，China

3. Computerized Tomography Lab，Capital Normal University，Beijing 100048，China

2014-02-08 收稿，2014-07-10 改回.

1 引言

新生代時(shí)期，青藏高原的隆升對亞洲內(nèi)陸干旱化、東亞季風(fēng)演化、全球氣候變冷等區(qū)域乃至全球尺度上的氣候環(huán)境演變均產(chǎn)生了重要的影響(An et al.，2001;Guo et al.，2002;Wu et al.，2012)，利用古高度計(jì)方法回溯大陸高原隆升歷史是探討其地球動(dòng)力學(xué)機(jī)制及氣候環(huán)境響應(yīng)的重要手段之一?；鹕饺蹘r流氣泡古高度計(jì)是國外近年來逐漸興起的一種新的探討古高度的研究方法(Sahagian and Proussevitch，2007)，其研究思路是:在測量火山熔巖流厚度、計(jì)算熔巖流氣泡體積的基礎(chǔ)上，定量計(jì)算熔巖流噴出時(shí)的古大氣壓強(qiáng)，再通過古大氣壓強(qiáng)換算成古高度。相比于其他古高度研究方法，火山熔巖流氣泡古高度計(jì)的優(yōu)勢包括:(1)適用范圍廣，可應(yīng)用于任何出露氣孔狀玄武巖的地區(qū);(2)由于新鮮火山巖是開展同位素測年的理想材料，因此熔巖流氣泡古高程計(jì)能夠給出令人信服的隆升事件的年代，這為研究高原隆升歷史提供了較可靠的年代學(xué)標(biāo)尺;(3)研究對象獨(dú)立于生物界之外，相對易于保存，更適用于高海拔、強(qiáng)烈侵蝕機(jī)制下的高原地區(qū)的古高度研究;(4)計(jì)算時(shí)，不需要假設(shè)古高度與溫度、緯度之間的關(guān)系，計(jì)算結(jié)果的影響因素少(Sahagian and Maus，1994)。目前，國內(nèi)科學(xué)家已成功地將火山熔巖流氣泡古高度計(jì)方法應(yīng)用到高原隆升研究中。例如，騰沖黑空山火山全新世熔巖流的古高度計(jì)算結(jié)果與目前實(shí)際海拔高度相吻合，表明火山熔巖流氣泡古高度計(jì)將成為研究青藏高原隆升歷史的有效手段之一(郭正府等，2011)。

熔巖流氣泡體積的計(jì)算精度直接影響著古高度計(jì)算結(jié)果的誤差，這是火山熔巖流氣泡古高度計(jì)實(shí)際應(yīng)用中的關(guān)鍵和難點(diǎn)。熔巖流氣泡體積的定量計(jì)算方法包括注膠法、巖石拋光-掃描法、體視學(xué)轉(zhuǎn)換法和三維CT (Computerized Tomography)掃描法(郭正府等，2011)。其中，測量和計(jì)算精度最高的手段是基于高分辨率X 射線層析成像技術(shù)和計(jì)算機(jī)算法的三維CT 掃描法。目前國內(nèi)已建立了相關(guān)的三維CT 檢測成像實(shí)驗(yàn)室，具備開展火山熔巖流氣泡古高度計(jì)研究的實(shí)驗(yàn)測試條件，但熔巖流氣泡的三維CT 掃描圖像解譯與氣泡體積統(tǒng)計(jì)的計(jì)算機(jī)算法實(shí)現(xiàn)仍然是亟需建立和完善的重要環(huán)節(jié)。鑒于此，本文將提出新的基于三維CT 掃描技術(shù)的熔巖流氣泡體積測量與統(tǒng)計(jì)算法，定量計(jì)算熔巖流氣泡的數(shù)量與體積，并探討該方法在火山熔巖流氣泡古高度研究中的應(yīng)用意義。

2 熔巖流氣泡體積統(tǒng)計(jì)算法

熔巖流氣泡的形態(tài)、數(shù)量、體積及分布特征對于探討巖漿氣泡的形成與演化、火山噴發(fā)動(dòng)力學(xué)以及利用熔巖流古高度計(jì)重建高原隆升歷史等均具有重要的作用。本文提出的基于三維CT 掃描技術(shù)的熔巖流氣泡體積算法旨在從熔巖流CT 體數(shù)據(jù)中分割出獨(dú)立氣泡，并統(tǒng)計(jì)獨(dú)立氣泡的個(gè)數(shù)和體積，包括獲取CT 體數(shù)據(jù)、數(shù)據(jù)分割、標(biāo)記連通區(qū)域、骨架提取、計(jì)算候選喉面面積、篩選喉面、統(tǒng)計(jì)氣泡體積等流程(圖1)。

2.1 數(shù)據(jù)獲取

計(jì)算機(jī)斷層成像(Computerized Tomography，CT)是基于射線與物質(zhì)相互作用的無損檢測技術(shù)，其原理為:X 射線穿過物體時(shí)會發(fā)生強(qiáng)度衰減，而衰減程度與物體的密度、原子序數(shù)以及射線能量和初始強(qiáng)度有關(guān)。利用這個(gè)原理，讓X 射線經(jīng)不同角度穿過被檢測物體并被探測器接收，利用算法重建被檢測物體各部分的射線衰減系數(shù)以得到其分布圖像，以此反映被檢測物體的密度分布。目前，三維CT 掃描技術(shù)已被廣泛地應(yīng)用到工業(yè)、醫(yī)學(xué)、生態(tài)學(xué)、考古學(xué)等多個(gè)領(lǐng)域(Zubal et al.，1994;Heeraman et al.，1997;Wu et al.，2008)，并已成為獲取高分辨率三維巖石數(shù)據(jù)、探討熔巖流古高度的有效手段之一(Song et al.，2001)。

與傳統(tǒng)的二維扇束CT 相比，錐束CT 具有射線利用率高、掃描時(shí)間短、各向分辨率一致等優(yōu)點(diǎn)。實(shí)際應(yīng)用中，通常采用立式平板探測器系統(tǒng)作為錐束CT 掃描的工作平臺(圖2)。檢測時(shí)，X 射線源和二維探測器都是固定的，待測物體在樣品臺上圍繞Z 軸旋轉(zhuǎn)，從不同角度接受X 射線的照射。本文所采用的三維巖石數(shù)據(jù)采集裝置為首都師范大學(xué)檢測成像工程中心研制的微焦點(diǎn)錐束CT 實(shí)驗(yàn)設(shè)備，射線源能量范圍為10 ～225kV，系統(tǒng)最高空間分辨率為3μm。掃描的氣孔狀玄武巖巖芯樣本直徑為25mm，掃描電壓能量為90kV。待測樣品投影數(shù)據(jù)的分辨率為1024 ×1024 ×720，即轉(zhuǎn)臺每旋轉(zhuǎn)0.5°采集一幅投影數(shù)據(jù)，其中每幅投影數(shù)據(jù)尺寸為1024 ×1024。

圖1 熔巖流氣泡體積統(tǒng)計(jì)的算法流程圖Fig.1 Algorithm flowchart for vesicle volume calculation of lava flows

圖2 錐束CT 掃描示意圖Fig.2 Diagram of the cone-beam CT scanning

圖3 氣孔狀玄武巖巖芯的三維CT 體數(shù)據(jù)(a)-熔巖流體數(shù)據(jù)三維視圖;(b)-熔巖流體數(shù)據(jù)二維斷層切片F(xiàn)ig.3 CT data of vesicular basaltic lava flows(a)-three-dimensional view of lava flow data;(b)-slice of lava flow data

FDK(Feldkamp-Davis-Kress)型算法是錐束CT 設(shè)備使用的主流重建算法之一(Feldkamp et al.，1984)。目前已在該套微焦點(diǎn)CT 設(shè)備上實(shí)現(xiàn)了基于GPU(Graphics Processing Unit)硬件加速的FDK 重建算法，包括濾波和反投影全過程(趙星等，2008)。重建的氣孔狀玄武巖三維CT 體數(shù)據(jù)的分辨率為512 ×512 ×512，其中每個(gè)體素的尺寸為49μm。圖3a為重建后的氣孔狀玄武巖三維CT 體數(shù)據(jù)，圖3b 為巖芯樣品的二維斷層切片體數(shù)據(jù)，三維CT 掃描圖像清楚地展示了熔巖流在冷凝、固結(jié)成巖過程中形成的氣泡及其形態(tài)、大小、分布等特點(diǎn)。

2.2 數(shù)據(jù)分割

圖4 氣孔狀玄武巖三維CT 體數(shù)據(jù)的閾值分割結(jié)果Fig.4 Results of threshold segmentation on 3D CT data of vesicular basalt

在采集并獲得氣孔狀玄武巖巖芯的三維CT 體數(shù)據(jù)后，利用不同物質(zhì)對射線吸收衰減系數(shù)的不一致性，將體數(shù)據(jù)分割為具有明顯灰度差異的巖石、氣泡和無效(背景)部分，以進(jìn)一步分析玄武巖中氣泡的形態(tài)、大小并計(jì)算氣泡體積?；谶@樣的三維CT 體數(shù)據(jù)特點(diǎn)，采用閾值分割方法處理氣孔狀玄武巖圖像，其原理為:參照不同物質(zhì)的灰度分布，選取恰當(dāng)?shù)拈撝礖0、H1。選定閾值后，遍歷體數(shù)據(jù)中的所有體素，將所有體素的灰度值與閾值H0、H1相比較。所有灰度值小于灰度閾值H0的體素，將被判定為無效(背景)部分;相反，所有灰度值大于或等于灰度閾值H1的體素，將被判定為巖石部分;其余體素則為氣泡部分(圖4)。因此，選取合適的灰度閾值是準(zhǔn)確判定氣泡形態(tài)、大小與體積的關(guān)鍵(Song et al.，2001)。

2.3 標(biāo)記連通區(qū)域

經(jīng)過閾值分割的體數(shù)據(jù)，各氣泡區(qū)域被分割成彼此分離的連通區(qū)域。所謂連通區(qū)域是指體數(shù)據(jù)在空間中的一個(gè)最大連通子集。連通性的定義如下:在一個(gè)連通集中，任意兩個(gè)元素之間都存在一條完全由這個(gè)集合中的元素構(gòu)成的路徑，可以表示為:兩個(gè)元素P、Q 是連通的，當(dāng)且僅當(dāng)存在一條路徑P1-P2-P3-…-Pn，使得P1=P，Pn=Q，?1≤i≤n-1 有Pi與Pi+1相鄰。

標(biāo)記連通區(qū)域就是尋找體數(shù)據(jù)中的所有氣泡對象，并且將屬于同一連通區(qū)域的所有體素用唯一的標(biāo)記值進(jìn)行標(biāo)記。本文采用基于深度優(yōu)先搜索原理的區(qū)域生長法對連通區(qū)域進(jìn)行標(biāo)記。首先對體數(shù)據(jù)進(jìn)行順序掃描，每遇到一個(gè)氣泡體素就分配一個(gè)標(biāo)記值，然后在該氣泡體素的鄰域進(jìn)行檢測。如果有尚未標(biāo)記的氣泡體素，則賦予相同的標(biāo)記值，反復(fù)進(jìn)行此操作，直到對體數(shù)據(jù)中所有體素掃描結(jié)束。標(biāo)記后的連通區(qū)域可能是獨(dú)立的單個(gè)氣泡，也可能是融合的多個(gè)氣泡(圖5)。為了準(zhǔn)確地提取氣泡連通區(qū)域的骨架，還需要將融合的多個(gè)氣泡分割開來，以提高氣泡數(shù)量和體積的計(jì)算精度。

2.4 骨架提取

圖5 連通區(qū)域標(biāo)記結(jié)果示意圖Fig.5 Tagging results of the connected regions

骨架是原始圖形的一種壓縮表示，與原始圖形保持了相同的拓?fù)浣Y(jié)構(gòu)，并且存在于圖形的對稱軸上，能夠反映圖形的拓?fù)渑c形狀信息。準(zhǔn)確地對標(biāo)記后的氣泡連通區(qū)域進(jìn)行三維骨架提取是生成和處理喉面的關(guān)鍵步驟。骨架提取算法通?？煞譃?拓?fù)浼?xì)化法、距離變換法、幾何分析法和廣義勢場法等(Cornea et al.，2007)。其中，拓?fù)浼?xì)化法的優(yōu)點(diǎn)在于可以較好保留原始圖形的拓?fù)涮卣鳎幢ＷC骨架連通性，同時(shí)也能保持骨架的單像素寬度。

本文采用Palágyi and Kuba(1999)提出的拓?fù)浼?xì)化法對氣孔狀玄武巖三維CT 圖像中的標(biāo)記連通區(qū)域進(jìn)行骨架提取。該算法的基本思想是從邊界開始向內(nèi)演化，利用反復(fù)迭代方法對三維模型進(jìn)行層層剝離。在每次迭代過程中，將邊界上的每個(gè)體素點(diǎn)及其26 鄰域的體素與3 ×3 ×3 模板進(jìn)行匹配，符合模板匹配的體素將被保留至下一次迭代過程，而不符合模板匹配的體素將被刪除(圖6)。以此類推，算法將逐步搜索到中軸骨架的位置，直至沒有任何點(diǎn)可以被刪除，最終獲得單像素寬的骨架(圖7)。

該算法是基于邊界的算法，只對邊界上的體素進(jìn)行模板匹配，極大地減少了計(jì)算量。由于算法的并行性，可在一次迭代過程中從邊界上刪除若干點(diǎn)，因而適合使用GPU 等并行處理器硬件提高運(yùn)行速率。

2.5 喉面處理

氣泡連通區(qū)域的骨架提取完畢后，計(jì)算各個(gè)骨架點(diǎn)在骨架線上的切線方向，也就是候選喉面的法平面方向。連通區(qū)域的骨架可以視為一條空間中的曲線，利用插值型的求導(dǎo)公式來計(jì)算骨架線在該點(diǎn)的切線方向。令x0，x1，x2，x3，x4為骨架上連續(xù)的5 個(gè)節(jié)點(diǎn)，則對應(yīng)每一點(diǎn)的求導(dǎo)公式為:

圖6 骨架提取算法匹配模板Fig.6 Matching template for skeleton extraction algorithm

圖7 連通區(qū)域骨架提取示意圖Fig.7 Skeleton extraction of connected regions

對于連通區(qū)域骨架線長度小于5 的情況，則利用兩點(diǎn)求導(dǎo)公式或三點(diǎn)求導(dǎo)公式(此處不再贅述)。切線方向確定后，也就能確定過該點(diǎn)且垂直于切線向量的法平面，稱為候選喉面。計(jì)算候選喉面面積是統(tǒng)計(jì)玄武巖氣泡體積的重要環(huán)節(jié)，本文提出一種新的方法——逐點(diǎn)法，用于候選喉面的面積計(jì)算。與其他計(jì)算方法相比，該方法的計(jì)算結(jié)果更為精確，且不受連通區(qū)域形狀的約束。對于計(jì)算得到的所有候選喉面面積，需從中篩選出真正的喉面用于氣泡體積計(jì)算，即確定融合氣泡的分割面。本文采用基于區(qū)域劃分的閾值計(jì)算方法，減少了傳統(tǒng)單一閾值分割方法導(dǎo)致的小連通區(qū)域過度分割的誤差，提高了氣泡數(shù)量和體積的統(tǒng)計(jì)精度。

2.6 統(tǒng)計(jì)結(jié)果

當(dāng)喉面確定后，即可對三維CT 體數(shù)據(jù)進(jìn)行重新標(biāo)記。在連通區(qū)域標(biāo)記值的基礎(chǔ)上，將獨(dú)立氣泡從連通區(qū)域中分離出來，并賦予新的標(biāo)記值。當(dāng)所有氣泡均被標(biāo)記后，則可統(tǒng)計(jì)獨(dú)立氣泡的數(shù)量與體積。不同的標(biāo)記值數(shù)量代表獨(dú)立氣泡的數(shù)量;被標(biāo)記為同一標(biāo)記值的體素則歸屬于同一個(gè)氣泡，因此可以統(tǒng)計(jì)出該氣泡的體素個(gè)數(shù)。最后，將數(shù)據(jù)獲取環(huán)節(jié)計(jì)算得到的氣泡分辨率物理尺寸換算為每個(gè)獨(dú)立氣泡對應(yīng)的實(shí)際物理體積，對氣泡體積進(jìn)行排序并分析其眾數(shù)分布特征，計(jì)算得到氣孔狀玄武巖巖芯的氣泡體積。

3 算法的主要優(yōu)勢

在本文提出的算法中，影響氣泡體積最終計(jì)算結(jié)果的關(guān)鍵步驟包括:(1)分割三維CT 體數(shù)據(jù)，以判定三維CT 掃描圖像中的巖石、氣泡和無效(背景)部分;(2)提取氣泡連通區(qū)域的骨架，以此構(gòu)建氣泡連通區(qū)域中單像素寬的連通線段;(3)計(jì)算候選喉面面積，即以骨架切線為法向量的截面面積;(4)篩選有效喉面，即從候選喉面中選取能夠代表氣孔狀玄武巖氣泡形態(tài)的分割面。其中，準(zhǔn)確選取氣泡連通區(qū)域分割面的操作難度最大，因此本文主要針對喉面面積的計(jì)算過程與有效喉面篩選方法進(jìn)行了優(yōu)化。

3.1 候選喉面面積計(jì)算

作為選取氣泡連通區(qū)域分割面的基礎(chǔ)，候選喉面面積計(jì)算是整個(gè)算法中最為重要的環(huán)節(jié)。多三角形近似法是目前應(yīng)用較多的計(jì)算喉面面積的方法之一(Shin et al.，2005)，其基本思想是以氣泡連通區(qū)骨架上的每個(gè)體素為頂點(diǎn)，在以該點(diǎn)切線方向?yàn)榉ㄏ蛄康钠矫嫔希鞫鄠€(gè)固定角度的射線。以相鄰兩條射線與巖芯邊界的交匯點(diǎn)作為兩個(gè)頂點(diǎn)，與該骨架點(diǎn)構(gòu)成了一個(gè)三角形(圖8)。計(jì)算每個(gè)連通區(qū)域的單像素寬度的連通骨架線e，對于每條骨架線e 上的點(diǎn)Ve，都用多三角形近似法計(jì)算該點(diǎn)的法平面面積。根據(jù)三個(gè)頂點(diǎn)的空間坐標(biāo)計(jì)算出單個(gè)三角形的面積，所有三角形的面積之和即為該截面的面積。

在單個(gè)骨架點(diǎn)法平面的截面中(圖9)，以骨架點(diǎn)為中心、45°為固定間隔角度向巖芯作射線。當(dāng)射線與巖芯邊界相交時(shí)，就得到一個(gè)截面邊界點(diǎn)Gj。以此類推，最終得到8個(gè)截面邊界點(diǎn)。利用骨架點(diǎn)與8 個(gè)截面邊界點(diǎn)組成的8 個(gè)三角形，計(jì)算所有三角形的面積之和，并將其作為候選喉面面積的近似值。多三角形近似法的精度可以通過不斷增加三角形的數(shù)量來逐步逼近截面的面積，所選則的角度越小，則得到的三角形越多，那么三角形面積之和也就越近似于候選喉面面積的實(shí)際值。

圖8 熔巖流候選喉面面積計(jì)算過程示意圖Fig.8 Calculation process of candidate throat surface in lava flow

圖9 多三角形近似法示意圖Fig.9 Diagram of the multi-triangular approximation

多三角形近似法適用于候選喉面輪廓比較規(guī)則的情況，但實(shí)際上熔巖流氣泡連通區(qū)域的截面往往是很不規(guī)則的(圖10)，若使用多三角形近似法計(jì)算候選喉面面積，則會產(chǎn)生較大的誤差。多三角形近似法的誤差主要由兩部分組成:(1)部分三角形計(jì)算結(jié)果偏大，對于氣泡連通區(qū)域內(nèi)凹的情形，統(tǒng)計(jì)時(shí)會將該部分巖石計(jì)入候選喉面面積;(2)部分三角形計(jì)算結(jié)果偏小，對于連通區(qū)域外凸的情形，統(tǒng)計(jì)時(shí)會忽略凸出的連通區(qū)域部分，不計(jì)入候選喉面面積。例如，圖10 中左側(cè)閉合部分為氣孔狀玄武巖的氣泡連通區(qū)域，閉合區(qū)域以外的部分為玄武巖。當(dāng)采用多三角形近似法計(jì)算氣泡連通區(qū)域面積時(shí)，得到的計(jì)算結(jié)果(即斜線陰影部分)將會包括多算的部分和少算的部分，造成難以估計(jì)的誤差(圖10)。

圖10 采用多三角形近似法計(jì)算截面面積導(dǎo)致的誤差示意圖Fig.10 Errors of cross section area resulting from multitriangular approximation

圖11 逐點(diǎn)法示意圖Fig.11 Diagram of traversal-points algorithm

鑒于此，我們提出一種新的計(jì)算候選喉面面積計(jì)算的方法——逐點(diǎn)法，以解決應(yīng)用多三角形近似法計(jì)算不規(guī)則氣泡連通區(qū)域截面面積時(shí)產(chǎn)生的不確定性誤差問題。首先，我們指定一個(gè)閾值Hpointonplane作為篩選有效體素點(diǎn)的判斷標(biāo)準(zhǔn)。然后，計(jì)算體素點(diǎn)到平面的距離，并與指定的閾值比較。

點(diǎn)到平面的距離公式:

其中，(a，b，c)為截平面的法向量，(x，y，z)為三維空間中一個(gè)點(diǎn)的坐標(biāo)。

當(dāng)點(diǎn)到平面的距離小于或者是等于指定的閾值時(shí)，認(rèn)為該點(diǎn)在平面上，即被判定為有效體素點(diǎn)。如圖11 所示，實(shí)心點(diǎn)代表落在平面上的有效體素點(diǎn)，而空心點(diǎn)代表非平面上的體素點(diǎn)。候選喉面面積越大，落在該平面上點(diǎn)的數(shù)目也會越多;反之亦然。

利用逐點(diǎn)法計(jì)算候選喉面面積的基本步驟為:

記Px為編號為x 的所有標(biāo)記為氣泡的體素集合;Ex為編號為x 的連通區(qū)域的骨架線;Vex，j為Ex上的第j 個(gè)骨架點(diǎn);Vpx，j為落在第x 區(qū)域第j 個(gè)點(diǎn)的截面的體素點(diǎn)集合;編號Ax，j為x 區(qū)骨架上第j 個(gè)點(diǎn)的截面面積。氣泡體素v∈Px，且Dv，Vex，j≤Hpointonplane，則v∈Vpx，j。將落在平面上的有效體素個(gè)數(shù)來作為候選喉面面積，則Ax，j= |Vpx，j|(表1)。

表1 應(yīng)用逐點(diǎn)法計(jì)算候選喉面面積的流程Table 1 Procedure of traversal-points algorithm for calculation of candidate throat surface area

雖然在應(yīng)用多三角形近似法計(jì)算候選喉面面積過程中，可以通過減小試探角度來提高精度，但如何選取試探角度以及確定射線與巖芯邊界的交點(diǎn)都是較難解決的問題，需要通過多次實(shí)驗(yàn)來調(diào)整合適參數(shù)以得到理想的計(jì)算結(jié)果。此外，由于不同成因熔巖流的氣泡形態(tài)、截面、分布特點(diǎn)等差別懸殊，若要通過多三角形近似法獲得較高精度的熔巖流氣泡體積，只能針對不同的熔巖流分別調(diào)整參數(shù)，因此多三角形近似法的適用性較差。

相較而言，逐點(diǎn)法減少了不規(guī)則氣泡連通區(qū)域形狀對候選喉面面積計(jì)算的干擾。更重要的是，影響逐點(diǎn)法計(jì)算結(jié)果的因素只有一個(gè)，即像素點(diǎn)到平面的距離閾值，而距離閾值的設(shè)定與氣泡截面的變化無關(guān)，使得逐點(diǎn)法可被方便地應(yīng)用于不同成因的熔巖流氣泡體積計(jì)算，其適用范圍比多三角形近似法更廣泛。此外，逐點(diǎn)法易于利用計(jì)算機(jī)算法實(shí)現(xiàn)，可以方便地移植到GPU 等并行計(jì)算平臺上。因此，在實(shí)際應(yīng)用中逐點(diǎn)法比多三角形近似法更具有優(yōu)勢。

3.2 喉面選取

喉面選取的原則是從骨架線上的截面中選擇最有可能成為氣泡之間分割面的平面。通用的作法是為所有的平面選擇一個(gè)全局閾值，小于閾值截面面積的平面將被視為分割面。然而，這種閾值選取方式忽略了一個(gè)重要的問題:不同連通區(qū)域的截面面積存在較大差異，而全局閾值方法可能導(dǎo)致氣泡被過度分割。例如，圖12 中標(biāo)記值為9 的連通區(qū)域的寬度明顯大于其它區(qū)域。如果對所有連通區(qū)域都采用同一個(gè)全局閾值，將會使一些體積較小的連通區(qū)域(如標(biāo)記值為1、2、6、8 的區(qū)域)分割過碎，導(dǎo)致無法準(zhǔn)確地選擇合適的喉面。解決過度分割問題的思路就是對不同寬度的氣泡連通區(qū)域采用不同的喉面閾值。因此，在本文中我們提出將閾值分段應(yīng)用于不同連通區(qū)域的方法——基于區(qū)域劃分的閾值計(jì)算方法。相對于傳統(tǒng)的僅采用單一閾值的判斷方法，基于區(qū)域劃分的閾值計(jì)算方法自適應(yīng)地調(diào)整了不同寬度連通區(qū)域的閾值，提高了分割精度。

圖12 氣泡連通區(qū)域的過度分割示意圖Fig. 12 Diagram of excessive separation on connected regions

首先，設(shè)置算法閾值Dneck(0≤Dneck≤1)。將不同連通區(qū)域骨架線上點(diǎn)的截平面劃分為不同的組，并對不同分組的候選喉面應(yīng)用閾值Dneck。通過這種設(shè)定，體積較大的連通區(qū)域所選用的喉面閾值就會相應(yīng)大一些，而體積較小的連通區(qū)域?qū)玫捷^小的閾值，從而緩解了僅采用全局單一閾值所產(chǎn)生的較小連通區(qū)域被過分分割的問題。然后，將Dax定義為標(biāo)記值x 對應(yīng)的氣泡連通區(qū)域的喉面閾值，其實(shí)質(zhì)為不同連通區(qū)域的骨架截面積之和乘以一個(gè)人工設(shè)定的閾值Dneck，計(jì)算公式如下:

如果?a∈Ax且a≤Dax(a ＞0)，則此時(shí)的截面a 即為有效的氣泡連通區(qū)域分割面。

4 程序運(yùn)算及結(jié)果

統(tǒng)計(jì)程序采用C+ +語言編寫，包括8000 余行代碼。通過程序運(yùn)算得到氣孔狀玄武巖巖芯的氣泡數(shù)量及體積，然后將氣泡體積進(jìn)行排序并分析其眾數(shù)分布特征(圖13)。

在生成眾數(shù)統(tǒng)計(jì)結(jié)果時(shí)，參照Sahagian et al. (2002)對橫坐標(biāo)取對數(shù)的方法。為了更好地分析統(tǒng)計(jì)結(jié)果，對于取對數(shù)后的橫坐標(biāo)乘以3，即拉伸3 倍以觀察氣泡體積眾數(shù)分布。從圖13 中可以得知，獨(dú)立氣泡體積Ln (volume)×3 的眾數(shù)分布區(qū)間為［6.5，7.5］，則獨(dú)立氣泡體積Volume 眾數(shù)分布區(qū)間為［9，12］(單位:體素)。在熔巖流氣泡眾數(shù)統(tǒng)計(jì)軟件計(jì)算結(jié)果的基礎(chǔ)上，利用火山熔巖流氣泡古高度計(jì)原理，郭正府等(2011)恢復(fù)了騰沖黑空山火山全新世噴發(fā)形成的氣孔狀玄武質(zhì)熔巖流的古高度，結(jié)果表明基于三維CT 掃描技術(shù)的熔巖流氣泡體積定量計(jì)算方法是有效和可行的。

圖13 氣泡眾數(shù)統(tǒng)計(jì)結(jié)果Fig.13 Results of vesicle volume mode

在算法運(yùn)算效率方面，對于512 ×512 ×512 分辨率的體數(shù)據(jù)，目前算法的計(jì)算耗時(shí)是約為6 小時(shí)(測試平臺為CPU:Intel Core2 Duo P8700@2.53GHz，內(nèi)存:4GB，硬盤320GB@5400RPM)。

5 結(jié)論

本文提出了基于計(jì)算機(jī)斷層成像技術(shù)的熔巖流氣泡眾數(shù)統(tǒng)計(jì)新方法，并針對影響熔巖流氣泡數(shù)量和體積統(tǒng)計(jì)精度的喉面面積計(jì)算與喉面選擇兩個(gè)關(guān)鍵技術(shù)環(huán)節(jié)提出了改進(jìn)方法，結(jié)果顯示:

(1)利用逐點(diǎn)法計(jì)算熔巖流體素點(diǎn)落在平面上的數(shù)量，以統(tǒng)計(jì)截面面積大小，能夠有效地克服多三角形近似法對于截面邊界形狀敏感的缺點(diǎn)，并提高喉面面積計(jì)算的精度;

(2)本文提出的基于區(qū)域劃分的閾值計(jì)算方法，解決了傳統(tǒng)喉面選擇方法中單一閾值所產(chǎn)生的小連通區(qū)域過度分割的問題，減小了氣泡數(shù)量統(tǒng)計(jì)的誤差。

根據(jù)上述新方法，我們設(shè)計(jì)并實(shí)現(xiàn)了熔巖流氣泡眾數(shù)統(tǒng)計(jì)軟件，其有效性和可行性得到了氣孔狀玄武巖巖芯樣品測試結(jié)果的驗(yàn)證。該方法有望提高熔巖流氣泡數(shù)量、體積以及火山熔巖流氣泡古高度計(jì)的計(jì)算精度。在后續(xù)工作中，需要通過進(jìn)一步的實(shí)驗(yàn)調(diào)整系統(tǒng)參數(shù)，減小算法誤差，以達(dá)到最佳的計(jì)算精度。在算法性能優(yōu)化方面，可利用并行計(jì)算硬件進(jìn)行算法加速(如多核CPU、GPU 等)，以提高統(tǒng)計(jì)運(yùn)算的效率。通過軟件與硬件兩方面的進(jìn)一步優(yōu)化，使熔巖流氣泡眾數(shù)統(tǒng)計(jì)方法成為高原隆升研究的快速分析工具。

致謝 熔巖流氣泡體積統(tǒng)計(jì)與計(jì)算方法由第二作者劉真完成;氣孔狀玄武巖三維CT 掃描測試與圖像處理由第三作者陳德峰完成;熔巖流氣泡測量與統(tǒng)計(jì)結(jié)果的地質(zhì)解譯由通訊作者和第一作者完成;樣品測試過程中得到首都師范大學(xué)張朋教授的熱心幫助;在此一并表示感謝。

An ZS，Kutzbach JE，Prell WL and Porter SC. 2001. Evolution of Asian monsoons and phased uplift of the Himalaya-Tibetan Plateau since Late Miocene times. Nature，411(6833):62 -66

Cornea ND，Silver D and Min P. 2007. Curve-skeleton properties，applications，and algorithms. IEEE Transactions on Visualization and Computer Graphics，13(3):530 -548

Feldkamp LA，Davis LC and Kress JW. 1984. Practical cone-beam algorithm. Journal of the Optical Society of America A，1(6):612-619

Guo ZF，Zhang ML，Cheng ZH，Liu JQ，Zhang LH and Li XH. 2011. A link of measurements of lava flows to Palaeoelevation estimations and its application in Tengchong volcanic eruptive field in Yunnan Province (SW China). Acta Petrologica Sinica，27(10):2863 -2872 (in Chinese with English abstract)

Guo ZT，Ruddiman WF，Hao QZ，Wu HB，Qiao YS，Zhu RX，Peng SZ，Wei JJ，Yuan BY and Liu TS. 2002. Onset of Asian desertification by 22Myr ago inferred from loess deposits in China.Nature，416(6877):159 -163

Heeraman DA， Hopmans JW and Clausnitzer V. 1997. Three dimensional imaging of plant roots in situ with X-ray computed tomography. Plant and Soil，189(2):167 -179

Palágyi K and Kuba A. 1999. A parallel 3D 12-Subiteration thinning algorithm. Graphical Models and Image Processing，61(4):199-221

Sahagian DL and Maus JE. 1994. Basalt vesicularity as a measure of atmospheric pressure and palaeoelevation. Nature，372(6505):449-451

Sahagian DL，Proussevitch AA and Carlson WD. 2002. Analysis of vesicular basalts and lava emplacement processes for application as a paleobarometer/paleoaltimeter. The Journal of Geology，110(6):671 -685

Sahagian DL and Proussevitch AA. 2007. Paleoelevation measurement on the basis of vesicular basalts. Reviews in Mineralogy ＆Geochemistry，66(1):195 -213

Shin H，Lindquist WB，Sahagian DL and Song SR. 2005. Analysis of the vesicular structure of basalts. Computers ＆ Geosciences，31(4):473 -487

Song SR，Jones KW，Lindquist WB，Dowd BA and Sahagian DL. 2001.Synchrotron X-ray computed microtomography: Studies on vesiculated basaltic rocks. Bulletin of Volcanology，63(4):252-263

Wu GX，Liu YM，He B，Bao Q，Duan AM and Jin FF. 2012. Thermal controls on the Asian summer monsoon. Scientific Reports，2:Article Number:404，doi:10.1038/srep00404

Wu XJ，Liu W，Dong W，Que JM and Wang YF. 2008. The brain morphology of Homo Liujiang cranium fossil by three-dimensional computed tomography. Chinese Science Bulletin，53(16):2513-2519

Zhao X，Zhang HT，Chen M et al. 2008. Progress in CT research at the computer tomography laboratory of Capital Normal University.Chinese Journal of Stereology and Image Analysis，13(3):158 -165 (in Chinese with English abstract)

Zubal IG，Harrell CR，Smith EO，Rattner Z，Gindi G and Hoffer PB.1994. Computerized three-dimensional segmented human anatomy.Medical Physics，21:299 -302

附中文參考文獻(xiàn)

郭正府，張茂亮，成智慧，劉嘉麒，張麗紅，李曉惠. 2011. 火山“熔巖流氣泡古高度計(jì)”及其在云南騰沖火山區(qū)的應(yīng)用. 巖石學(xué)報(bào)，27(10):2863 -2872

趙星，張慧滔，陳明等. 2008. 首都師范大學(xué)檢測成像實(shí)驗(yàn)室CT 研究進(jìn)展. 中國體視學(xué)與圖像分析，13(3):158 -165