程榮， 錢生平， 孫添力， 周懷陽，*

(1.同濟大學海洋與地球科學學院， 上海 200092;2.同濟大學海洋地質國家重點實驗室， 上海 200092;3.同濟大學電子與信息工程學院， 上海 201804)

氣孔構造是火山巖的主要代表性特征之一。硅酸鹽巖漿在上升過程中，壓力降低，揮發(fā)份過飽和后發(fā)生氣體出溶，并經過氣泡長大、合并、變形，或者氣體逃逸的改造，形成火山巖的氣孔構造[1-3]。根據巖漿中揮發(fā)份壓力與環(huán)境壓力平衡的原理，利用火山巖中氣孔的大小分布可用來估計大陸溢流玄武巖地區(qū)的古大氣壓和古海拔[4-6]。對火山巖氣孔特征的準確評估可以提供有關巖漿上升速度和噴發(fā)方式的信息，從理論上確定火山噴發(fā)的時間和噴發(fā)的性質，提高預測火山爆發(fā)的能力[7-9]。不僅如此，對火山巖氣孔特征(氣孔的大小、空間和形狀分布，整體的連通性、定向性等)的研究，還可以幫助人們了解巖漿的原始揮發(fā)份含量，源巖(區(qū))性質，巖漿的物理(黏性、流動性等)性質，以及氣泡的形成和生長過程[2,10-12]。

人們很早就開始使用各種實驗方法研究火山巖中氣孔的物理特征。Houghton等[13]根據阿基米德原理，通過分別測量火山碎屑巖在空氣和水中的質量計算巖石的孔隙率(水浸法)。Sahagian等[14]在水浸法的基礎上，先使用低黏度聚合物(液體環(huán)氧樹脂)充填氣孔，然后用氫氟酸溶解巖石基體，最后通過對殘余填充物進行計數和測量獲得氣孔的大小分布。Whitham等[15]使用樹脂進行真空浸漬，證實浮石中所有氣孔是相互連接的，并進一步通過氮吸收技術和壓汞法分別得到氣孔的表面積和不同孔徑氣孔的孔隙率分布。Shea等[16]基于二維圖像數據，將氣孔視為規(guī)則球體進行體視學轉換來分析火山巖氣孔的三維分布特征。隨著光學和電子顯微鏡成像技術的應用和普及，人們從依據薄片的光學照片[17]和背散射電子圖像[18]統(tǒng)計氣孔大小分布和孔隙率，發(fā)展為使用立體掃描電鏡(SSEM)構建樣品表面的三維數字高程模型來觀察氣孔形態(tài)和計算氣孔體積[19]。然而，這些方法通常很耗時，采集的數據量較少，并且會破壞樣品的原有組成和結構，其推廣應用受到一定的限制。

高分辨率X射線計算機斷層掃描，也稱工業(yè)CT，是一種現(xiàn)代化的材料檢測技術。目前，該技術已經廣泛應用于醫(yī)學、工業(yè)、考古學、生態(tài)學等多個領域，并且正在逐步發(fā)展為地球科學領域的重要分析手段[20-21]。自Song等[22]的開創(chuàng)性工作以來，越來越多的研究使用工業(yè)CT討論火山學問題，這些研究大多集中在氣孔結構和大小分布的測量上，并強調相互貫通的氣孔構造對巖漿滲透率和氣體散失的影響，顯著提高了人們對火山脫氣動力學的理解[23-26]。對巖石氣孔進行分析時，工業(yè)CT具有無損、快速、實時、高分辨率(微米到亞微米級)立體成像的特點，相比于傳統(tǒng)的稱重法、浸泡法、注膠法和分層磨片光學成像法具有明顯的優(yōu)勢，可以讓研究者對巖石的內部結構(即晶體、氣孔和基質)進行三維空間的可視化和量化，并通過上述任何一種技術直接測量同一樣品而獲得可重復的結果和補充性研究[27-30]。

本文通過將工業(yè)CT掃描與計算機圖像處理技術相結合，使用了三種分析火山巖中氣孔含量及大小分布的軟件或程序(VG Studio MAX、ImageJ、Matlab代碼)，此過程能快速、便捷地提取出氣孔可靠的形態(tài)學數據；同時，通過改進火山巖中氣孔體積轉化為揮發(fā)份質量分數的方法，定量評價揮發(fā)份在巖漿中的富集程度。

1 實驗部分

1.1 樣品來源與實驗設備

實驗樣品為2018年“嘉庚號”無人遙控深潛航次在南海珍貝海山通過水下機器人(ROV)抓取到的氣孔狀玄武巖(圖1a)。樣品表面呈深灰色，塊狀構造，斑狀結構，具輕微蝕變，基質主要由微晶組成，常見的斑晶礦物有橄欖石、輝石、長石等。由X射線熒光光譜儀(XRF熔片法)測試的全巖主量成分含量為：SiO244.63%，Al2O315.87%，F(xiàn)e2O311.11%，MgO 6.89%，CaO 8.69%，Na2O 3.96%，K2O 2.76%，P2O50.83%，MnO 0.28%，TiO22.94%。樣品中的氣孔構造十分發(fā)育，大多數氣孔為扁圓形或不規(guī)則，且相互獨立，與一般玄武巖中多為偏圓形氣孔有所不同。少數氣孔有合并現(xiàn)象，在空間中隨機分布，既沒有聚類也沒有規(guī)則的間距。氣孔尺寸相差較大，大的可達2～3mm，小至微米級?？紤]到樣品大小與CT掃描的分辨率密切相關，以及獲取一系列規(guī)則切片的需要，本實驗將檢測樣品切割成19mm×22mm×43mm的四方柱體。

a—“嘉庚號”無人遙控深潛航次ROV水下采樣現(xiàn)場情況，采樣水深1488m；b—微焦點X射線CT掃描成像系統(tǒng)示意圖，XYZ線代表三維空間坐標軸，探測器用于接收投影數據。圖1 氣孔狀玄武巖樣品和檢測設備Fig.1 Vesicular basalt samples and testing equipment

a—樣品三維CT重建效果； b—巖石中的氣孔提取過程，藍色為渲染的氣孔； c—去除巖石基體后的氣孔三維立體分布； d—CT切片的灰度直方圖，氣孔與巖石基體灰度差異大，容易區(qū)分； e—CT切片圖，樣品區(qū)域內黑色為氣孔，灰色為基質和低密度礦物斑晶，白色為高密度礦物斑晶； f—切片中的氣孔分布圖； g—批量處理CT切片的Matlab代碼運行框架和操作流程。圖2 基于工業(yè)CT掃描的三種氣孔含量及大小分布統(tǒng)計方法Fig.2 Types of methods used to estimate content and size distribution of vesicles based on X-ray CT

玄武巖樣品在中國石油化工股份公司石油物探技術研究院使用微焦點錐束CT實驗設備(美國通用電氣公司)進行測試，設備型號為Phenix v|tome|x s，測試電壓150kV，束流220mA，檢測樣品的系統(tǒng)分辨率為19.9μm。典型的CT布局系統(tǒng)如圖1b所示，X射線源和探測器分別置于轉臺兩側，轉臺帶動樣本360°旋轉，每轉動一個微小的角度后，探測器會接收由X射線照射樣本獲得的投影圖數據。其原理是利用X射線在穿過物體的斷面時，會因物質密度的不同而發(fā)生不同程度的衰減，計算機系統(tǒng)通過將衰減系數轉為灰度值，再經過一系列校正后就可以重構樣品的三維像素體圖像[29]。因此，CT圖像可以反映被檢測物體內部物質密度的分布關系。對于具有較大氣孔和簡單氣孔結構的氣孔狀玄武巖，由于氣孔與基體密度差異大，CT圖像上的灰度會明顯不同(圖2d)。一般氣孔呈現(xiàn)灰度值較低的黑色，圖像軟件很容易借助灰度閾值分割，邊緣檢測，形態(tài)學運算(開閉、腐蝕膨脹等)等方法提取出巖石內部的氣孔形貌，這一轉化過程所帶來的誤差也可忽略[31-32]。

1.2 氣孔含量及大小分布統(tǒng)計方法

在確保工業(yè)CT的X射線束空間分辨率足以解析火山巖氣孔特征的前提下，使用配套商用軟件VG Studio MAX(https://www.volumegraphics.com/en/products/vgstudio-max.html)重構樣品的三維立體圖像，并調節(jié)灰度閾值范圍分離出氣孔，測量其結構形態(tài)等參數(圖2a～c)。然后，由該3D數據處理軟件輸出樣品主視、側視、俯視三個方向的二維CT切片，并通過開源軟件ImageJ(http://rsb.info.nih.gov/ij/)對CT切片逐個進行圖像處理和氣孔形態(tài)學運算，從而在二維平面上更為精細地展現(xiàn)氣孔細節(jié)，研究氣孔在巖石內部的空間分布差異(圖2d～f)。

盡管使用ImageJ處理切片可以獲取所需的氣孔含量、大小及形態(tài)數據，但這對于處理高精度CT掃描獲得的大批量CT切片來說：①操作較為繁雜，耗時費力；②批量處理時占用電腦內存過大，不適用一般配置的個人電腦；③軟件本身對單次處理文件也有一定的限制。為此，作者編寫了一套專門處理CT切片數據的Matlab程序。本程序在Win10 X64位操作系統(tǒng)、Intel(R) Core(TM) i5-7200U CPU@2.5GHz2.71GHz、RAM 4GB的運行環(huán)境下，使用MATLAB R2018b批量處理900張1266×2204分辨率TIF格式的主視方向切片圖僅耗時6.8min(其中CT Process為4.5min)，而且操作簡便，占用內存小，可以隨時中斷然后繼續(xù)運行(源代碼見附件)。

該程序的運行框架和操作流程如圖2g所示。首先，創(chuàng)建CT_image、binary_image、mask_image三個文件夾，并把CT_process.m、area_calculation.m、hole_stats.m三個Matlab代碼置于同級目錄下，然后將CT切片放在CT_image文件夾中，依次運行三個代碼，最終輸出每張切片中氣孔的統(tǒng)計結果。其中，CT_process算法是該計算程序主要的圖像處理部分。具體過程是：先將CT切片圖轉換為8位的二值化灰度圖，在區(qū)分邊界后，創(chuàng)建蒙版保護背景部分不被處理，再使用圖像閉操作算法獲取樣本區(qū)域的大小，最后輸出二值化氣孔圖和樣本區(qū)域圖，并依次循環(huán)處理下一張切片。在此基礎上，area_calculation和hole_stats算法通過對氣孔圖和樣本區(qū)域圖中像素點的運算、統(tǒng)計，可以計算出感興趣區(qū)的氣孔含量和氣孔大小分布情況。因為CT切片中灰度值的絕對大小還與原子序數和X射線能量相關，像元灰度值只能反映各物質密度的差異，并沒有嚴格的物理意義[33]，所以在此之前需要基于VG Studio MAX和ImageJ的氣孔含量測定結果(與之相等或相當)，采用逆分析的辦法手動設定統(tǒng)一的灰度閾值分割點，以此減少因為閾值選取不恰當導致氣孔提取偏少或部分基體也被提取出來的影響。

1.3 氣孔體積-揮發(fā)份質量分數轉換方法

大量研究表明，巖漿中的揮發(fā)份主要是二氧化碳和水蒸氣，其含量和存在形式能夠顯著影響巖漿的流變學性質，在控制火山噴發(fā)行為中起著重要作用[34-36]。為定量分析揮發(fā)份在巖漿中的富集程度，將一定體積的氣孔含量轉化為火山巖中揮發(fā)份的質量分數是必要的。在假定溶解在熔體中的氣體遵循“亨利定律”的前提下，Head等[37]將火山巖氣孔的體積分數視為巖漿中逸出CO2或H2O的理想氣體的含量，計算海底火山在臨界爆破條件下不同水深的氣體質量分數。但是，該方法采用的氣體密度是通過理想氣體狀態(tài)方程得到的，隨著水深加深壓力增大，氣體性質會逐漸偏離理想氣體，最終得到的氣體質量分數也會偏離實際情況。

已知平衡狀態(tài)下一種氣體的密度只與飽和溫壓有關，我們可以通過調用美國國家標準與技術局(NIST)熱物性數據庫中給定溫壓下的氣體密度，對火山巖氣孔體積-揮發(fā)份質量分數轉換方法加以改進。該NIST熱物性數據庫(https://www.nist.gov/srd/refprop)使用的是當前可用的、最精確的氣體狀態(tài)方程和模型，收集的數據經過了嚴格的評估和審查，理論上會更為準確。如表1所示，本文對方法作了改進之后，轉換結果與Head等[37]的方法相比，在同樣的條件下，兩者計算值的相對標準偏差(RSD)均在合理的范圍內，在水深不超過3500m時，CO2含量小于4.30%，H2O含量小于0.95%，CO2偏差更大，并且隨著水深加深壓力增加，氣體狀態(tài)逐漸偏離理想氣體的緣故，RSD有進一步擴大的趨勢。因此，改進后的方法并非基于簡單理想氣體狀態(tài)方程，而是采用數據庫氣體密度計算的方法得到的氣體質量分數更加準確，特別是對于水深較大的海底火山巖樣品。

2 結果與討論

2.1 氣孔的三維空間分布特征

在使用VG Studio MAX重構樣品的立體空間后，可以直接獲得氣孔的個數、體積、表面積、直徑等參數。通過該方法得到玄武巖樣品的氣孔體積分數為12.32%。如圖3a所示，氣孔尺寸分布呈現(xiàn)出對數正態(tài)分布的特點，相關系數R2皆大于0.96，這種分布可能反映巖漿中的氣泡是由簡單的成核和生長事件引起的，偏右分布的大氣泡是許多小氣泡發(fā)生聚集合并的結果[39]。其中，與氣孔體積相等的等效球直徑分布在100～800μm區(qū)間內，在180～200μm之間的數量最多；氣孔最大外接圓直徑(即最大直徑)分布在140～1500μm區(qū)間內，分布區(qū)域相對較廣，在340～360μm之間的數量最多。最大直徑比等效球直徑的分布明顯偏右，這應該是氣泡生長受到巖漿性質(黏度、流動性)影響導致氣泡變形的結果[11,40]。如圖3b所示，氣泡從成核到不斷長大的過程中，比表面積會不斷減小，在達到一定尺度后趨于穩(wěn)定，穩(wěn)定值約0.02μm-1。這種變化特征除了反映氣泡比表面積與直徑近似呈反比例函數的關系，也在一定程度上反映了氣泡生長向規(guī)則球形轉變的趨勢，因為氣泡在應力松弛時本質上有利于巖漿-氣泡界面上應力的均勻分布，從而減少其表面能[16,41]。

表1 本文與其他學者獲得的氣孔體積-質量分數轉換結果對照

Table 1 Comparison of conversion results of vesicle volume to mass fraction between this paper and another scholar

水深(m)氣體全部為CO2(%)氣體全部為H2O(%)本文Head等(2003)[37]RSD(%)本文Head等(2003)[37]RSD(%)5001.9131.9350.570.8010.8020.0610003.6783.7611.121.5751.575015005.3435.5201.632.3352.3370.0420006.9147.2162.143.0843.0870.0525008.4018.8522.613.8213.8250.0530009.80910.4323.084.5454.5530.09350011.14512.1474.305.2585.3590.95

注：氣體體積分數ηCO2/H2O(vol/%)=ηvesicle(vol/%)=75%，深水壓強P(MPa)=1026(kg/m3)×9.8(m/s2)×Depth(采樣水深,m)/106+0.101325(MPa,標準大氣壓)，巖漿溫度T=1255℃，熔體密度ρlava=2700(kg/m3)；氣體密度ρCO2/H2O(kg/m3)采用的數值不同，Head等[37]由理想氣體狀態(tài)方程計算得到，本文是先在系統(tǒng)中安裝好Nist Refprop9.1，并將其加載到Excel宏中，然后在單元格中輸入氣體的密度函數[38]調用NIST數據庫特定溫壓下的氣體密度得到；巖石密度ρrock(kg/m3)=ρlava(kg/m3)×[1-ηCO2/H2O(vol/%)]+ρCO2/H2O(kg/m3)×ηCO2/H2O(vol/%)，氣體質量分數ηCO2/H2O(wB/%)=ρCO2/H2O(kg/m3)×ηCO2/H2O(vol/%)/ρrock(kg/m3)×100%。RSD為兩種方法CO2質量分數計算值的相對標準偏差。

a—氣孔尺寸分布直方圖，紅線為體積等效球直徑擬合曲線，綠線為最大外接圓直徑擬合曲線；b—氣孔尺寸與比表面積關系，比表面積為氣孔表面積與體積之比。圖3 氣孔的三維空間分布情況Fig.3 Three-dimensional spatial distribution of vesicles

2.2 氣孔的二維切片面分布特征

通過ImageJ和Matlab處理主視、左視、俯視三個方向上的二維切片，獲得了氣孔含量和氣孔數密度在剖面上的變化特征(圖4)。二維切片統(tǒng)計的氣孔含量在各個剖面上圍繞三維體積分數12.32%上下波動。兩種方法在三個方向上所有切片氣孔的平均含量分別為主視12.69%/12.48%、左視12.66%/12.50%、俯視12.66%/12.51%，與三維統(tǒng)計結果的相對標準偏差(RSD)均小于1.48%，并且兩種方法統(tǒng)計的氣孔含量的變化趨勢和數值大小基本重合，同步性近乎一致。這說明基于切片圖作數據再處理的方法在統(tǒng)計精度上是穩(wěn)定可靠的。

從剖面方向來看，隨著深度加深，氣孔含量在主視、左視方向上的切片波動較大，但無明顯分段，在俯視方向上的切片大體可以分為三段：10mm以下13%～16%，10～35mm之間11%～13%，35mm以上11%～14%，氣孔含量整體呈現(xiàn)先減小后增加的趨勢。而且，三個剖面上的氣孔含量變化基本與氣孔數密度(單位面積的氣孔個數)的變化正向同步，表明氣孔數密度是決定氣孔含量大小的決定性因素。少數不同步的區(qū)間段(圖4a，8～16mm；圖4b，4～10mm)是因為局部切片中大氣孔與小氣孔占比差異較大的緣故，氣孔數密度比較一致的情況下，大氣孔多的氣孔含量會有所升高。

a—主視圖方向剖面；b—左視圖方向剖面；c—俯視圖方向剖面。

圖4 二維切片中氣孔含量和氣孔數密度的面分布特征

Fig.4 Profile distribution features of vesicle content and vesicle number density in two-dimensional slices

火山巖中氣孔的大小，形狀和空間分布記錄了巖漿脫氣過程中氣泡的形成和生長歷史，最終的氣泡尺寸分布取決于成核的早晚和次數，并受氣泡生長、合并和變形的影響[3,14,39]，這些過程與巖漿性質、氣體擴散速率、溫壓條件變化等因素密切相關[1,11,42]，所以該樣品中氣孔面分布特征(含量、數量、大小)的較大差異，表明氣泡的成核和生長過程受到了多種因素的影響。

2.3 揮發(fā)份的質量分數

根據第1.3節(jié)改進后的氣孔體積-揮發(fā)份質量分數轉換方法，計算南海氣孔狀玄武巖樣品的揮發(fā)份質量分數。將采樣水深近似等于古水深Depth=1488m，巖漿溫度和熔體密度根據全巖主量成分數據，使用Kurt Hollocher設計的CIPW norm Program計算(Fe3+/Fe=0.12)，其中T=1286℃，ρlava=2971kg/m3，氣孔體積分數采用三維統(tǒng)計值η(vol/%)=12.32%，最終計算出氣體的質量分數為ηCO2(wB/%)=0.233%，ηH2O(wB/%)=0.099%。若巖漿逸出的氣體全部為CO2，揮發(fā)性氣體過飽和點立刻成核，根據玄武質巖漿中CO2的溶解度公式，可以得到巖漿開始脫氣釋放出CO2，大致發(fā)生的最大深度為海底以下3829m(圖5)。

巖漿的脫氣過程與揮發(fā)份在巖漿中溶解度的變化密切相關，在恒定的溫度下，揮發(fā)份的溶解度受壓力和熔體成分的控制[34]。當巖漿沿著裂隙通道上升侵位或噴出地表時，由于玄武質巖漿中CO2含量高，溶解度遠比H2O小(圖5)，因此在海底高壓(>10MPa)環(huán)境下，巖漿脫氣進入氣泡中的揮發(fā)份絕大部分是CO2[43-45]。盡管CO2的擴散系數相對較小，但仍可快速實現(xiàn)熔體溶解的CO2與氣泡之間的平衡，所以對于一般情況下的熔體上升速率(未爆破)，巖漿能夠保持近似平衡脫氣的過程[44,46]。

溶解度公式：n=0.0023P(CO2)，n=0.1078P0.7(H2O)。式中：n(wB/%)代表溶解氣體的量，壓力P(MPa)=Depth(m)×9.8(m/s2)×2700(kg/m3)/106。

圖5 氣體在玄武質巖漿中的溶解度與壓力的關系(據Parfitt等[47]修改)Fig.5 Relationship between gas solubility in basaltic magma and pressure (Modified according to Parfitt,et al[47])

2.4 方法的適用性和有待改進之處

基于高分辨率工業(yè)CT掃描，本文使用的三種統(tǒng)計氣孔的軟件程序各有其利弊。其中：①商業(yè)化軟件VG Studio MAX可視化功能強大，可以進行樣品的三維重建，但無法對內部細節(jié)逐個分析，功能也相對有限，且數據處理成本較高；②ImageJ是開源軟件，免費使用，功能齊全，對圖像格式限制少，但操作繁雜，效率較低，占用內存大，不適用于大量切片的批量處理；③Matlab代碼的圖像處理功能和氣孔提取算法，可以快速獲取巖石內部不同方向上各個剖面的氣孔分布情況，批量處理，操作簡便，占用內存小，針對個性化的需求，還可再度編譯。

上述氣孔統(tǒng)計方法仍存在一些局限性，如無法區(qū)分氣孔中充填的杏仁體，在通過灰度閾值分割提取切片中的氣孔時會存在視覺誤差。這對于氣孔尺寸大且結構簡單的玄武巖樣品來說影響可以忽略不計，但如果是微氣孔含量高、氣孔結構復雜的樣品，最直接的解決辦法是使用更高精度的工業(yè)CT進行測試，或從切片重構算法、三維數據逆處理、多尺度單元體表征三個方面對CT檢測和數據處理過程進行改進[28]。對于巖石中氣孔相互連接組成的多孔網絡，則需要基于“骨骼”的概念設計算法將單個氣孔分離、標記和計數，從而了解多孔空間的拓撲結構[48]。這些信息可以與巖漿的物理、流變和化學特性相結合，為模擬和理解巖漿過程提供基礎。

火山巖氣孔體積轉化為揮發(fā)份質量分數的方法假定的是封閉巖漿體系，巖漿從洋殼深部到海底表層發(fā)生平衡脫氣，氣體全部為單一氣體的理想情況。對于實際的巖漿脫氣過程則需要考慮更多的因素；如果能選取具有代表性的多個火山巖(確保低中高氣孔含量的樣品都能采集到)，并將氣體之間的比例、氣泡定型時上覆巖層的壓力、古水深的變化代入運算，那么由火山巖氣孔得到的氣體質量分數才可能代表真正意義上的巖漿脫氣量。

3 結論

工業(yè)CT能夠對巖石進行三維成像，可視化樣品內部的空間結構，具有分辨率高和無損樣品的特點。三種氣孔統(tǒng)計方法能夠快速提取巖石中的氣孔，獲取氣孔的含量和大小分布情況。其中，商用軟件VG Studio MAX適用于氣孔的三維重構和體積測量，開源軟件ImageJ適用于少量CT切片的圖像處理和氣孔的二維形態(tài)學運算，自主開發(fā)的程序代碼使用方便快捷，適用于大量CT切片的批量處理。通過調用美國國家標準與技術局熱物性數據庫的氣體密度，改進的海底表層環(huán)境下火山巖中氣孔體積-揮發(fā)份質量分數轉換方法，在深水條件下氣體質量分數計算結果更加準確。

本文方法還可通過計算機算法分割巖石中相連的氣孔和去除充填的杏仁體，并根據熔體包裹體成分和各氣體的溶解度關系推測每種氣體所占的比例來進行改進，從而反映更加真實的氣孔構造特征和揮發(fā)份含量的變化，幫助人們了解火山噴發(fā)和巖漿脫氣的動力學過程。