现代地球物理探测技术的高速发展使岩石圈及更深尺度精细三维结构数据(地震波速度、密度、电、磁等)的获取成为可能^[1~5]。然而纵观宏观固体地球物理现有出版物,固体地球物理三维数据展示基本限于三维模型外表面或平面和剖面等值云图显示^[3,6],三维数据较少以整体方式显示其分布特征。另一方面,断层线及浅部地质资料多为后期制图时添加以解释数据,深部地球物理探测数据多数情况下并未真正与浅部地质资料融合。因此本文给出一套岩石圈尺度三维结构模型综合及可视化方案,对深部地球物理数据和浅部地质资料进行融合,并展示基于体绘制技术(Volume rendering)的三维可视化效果。

科学数据可视化是近些年发展迅猛的一个新的学科分支。可视化的目的在于以图形方式显示数据从而帮助我们展示、理解、剖析数据的内涵。以体绘制和等值面(Isosurface)技术为代表的三维数据可视化由于包含了大量计算机新技术、新算法和特定设备而形成了一门学科。科学数据可视化已经广泛应用于物理学(天体物理、引力波、分子动力学)、地理学、生态学、数学、大气科学(https://en.wikipedia.org/wiki/Scientific_visualization)以及医学与岩石(微观)CT (Computed Tomography)图像的展示等方面^[7~10]。

本文尝试将三维数据可视化技术应用于地质和地球物理综合数据,以青藏高原东缘(100°~105°E,26°~32°N)地区为例实现该方案。用于展示的地球物理数据包括P波、S波速度以及密度;浅部地质资料侧重断层数据(断层几何形态的三维表征)及地表高程数据。具体实施方案为利用Paraview软件综合上述数据,并主要采用体绘制技术显示特定参数的异常特征。

建立青藏高原东缘(川西地区)岩石圈三维结构主要用到以下3种数据。

(1) 地球物理数据。中国地震局地质研究所地震动力学国家重点实验室在川西地区(100°~105°E,26°~32° N)布设了由297个宽频带地震台组成的密集流动观测台阵^[11],该台阵布设于2006—2009年,期间获得了大量远震观测数据。通过远震接收函数与环境噪声数据联合反演得到该地区速度结构数据^[3],包括横波波速(Vs)、纵波波速(Vp)及密度(D)3个参数,在深度上均由地表延伸至地下100 km处。台阵密度非常之高,所获得的速度结构数据可能是迄今为止分辨率最高的三维地壳速度场——地震波速度和密度数据的横向分辨率依赖于台站距离,介于10~30 km;垂向分辨率均为2 km。

(2) 浅部断层数据。从已发表论文中收集相关断层资料并提取断层的走向、倾向、倾角数据^[12~14],确定研究区主要断裂带的三维空间分布。目前断裂的信息还在不断完善中,尤其地表以下信息(倾角变化及延伸深度)不确定性较大。本文主要目的在于展示新方法新技术的应用方案及成果,断裂数据不完整并不影响主要信息的传递。

(3) 高程数据。来源于“地球系统科学数据共享平台”分辨率为1 km的网格数据,选取其中位于研究区范围内(100°~105°E,26°~32°N)的高程数据。

研究包含2个主体内容:①三维结构模型的综合;②三维模型可视化。选用开源软件Paraview作为本项研究的软件平台。

Paraview(http://www.paraview.org/)是一个开源的、跨平台的高性能数据分析与可视化应用软件,由Kitware技术公司开发。支持超过100种文件格式的读取,数据处理和显示可以通过交互式操作或通过程序批处理完成。Paraview采用分布式内存计算资源,可以分析非常大的数据集,在超级计算机和笔记本电脑上皆可运行。Paraview可以对任何点、线、面、体数据进行建模,标量、矢量及张量数据皆可表达。通过过滤器转换数据或对体模型进行体绘制以提取重要信息,使体模型的内部结构也能够直观地表现出来。Paraview软件可以在运行中同时创建多个视图以同步显示不同模型或图表。对一个文件设置多个时间步能够自动生成动画,方便用户导出。

常规的三维地质建模软件有GOCAD^®,Geomodeller^®,MaoInfo Discover^®和国产的MapGIS^®等。此类软件普遍具有三维建模、可视化、地质解译和分析的功能,以实现对表面或实体的建模与观察,在地质工程、地球物理勘探、矿业开发、遥感后期处理等方向上均有广泛应用。然而,此类软件均不具有体绘制功能,对体模型的直观观察稍有欠缺。

目前可采用的包含体绘制功能的软件较多,典型的有Paraview,ImageJ,VisIt,Drishti,Mayavi,Unidata提供的可视化工具(IDV等)以及商业软件Avizo^®和Amira^®。总体而言,非商业软件中,Paraview和VisIt的功能、用户数量和发展速度都名列前茅;而Unidata作为较新的平台,目前发展十分迅速。基于易用性、数据集成功能和开源特征,本文使用Paraview实现数据综合与三维可视化。方法与思路同样可用于其他类似软件平台。

三维模型一般由点、线、面或体的实体及其组合构成,有多种表达格式。研究中用到VTK(Visualization Toolkit)和RAW 2种文件格式。

VTK是以点要素为基本单位的数据存储格式,通过导入点的位置、点的排列方式及每个点的参数信息建立模型。VTK数据模型可以有效地表示直线、曲线网格、体数据、表面网格和三维非结构网格等不同形状和大小的网格。地表高程数据和断层形态建模采用了三维非结构网格(UNSTRUCTED GRID)格式。三维非结构网格的数据存储分为4个部分^[15]:①点数据;②单元数据(单元数据可以是线、面、体等实体,每个单元数据由特定排列方式的点数据构成);③单元数据的类型(即点的排列方式,包括线性单元类型及非线性单元类型);④对应点或单元的标量、矢量或张量参数(如密度、速度、应变等)。VTK格式对数据的描述较为灵活,但数据存储所需空间相对较大。

RAW格式是一种局限于沿坐标方向尺度确定的规则空间内(如二维空间的矩形或三维空间的长方体)的数据储存格式,空间内部被划分为均匀小格。三维RAW数据只需给出长方体内均匀小立方体(称为体像素)上的变量值(如P波速度、S波速度、密度等)。体像素是相对于二维图像数据中像素概念的扩展,一般理解为小立方体,实际应用中也可以为小长方体。不同于被广泛接受的BMP和JPG等基于像素概念的二维(图像)数据,RAW格式可以不包含文件头且没有任何压缩,并完全按照x轴、y轴、z轴正方向的顺序从原点开始,依次存储每个体像素中相应的参数数据。RAW格式完全没有存储冗余信息,占用存储空间小、读取速度快。岩石圈速度和密度结构数据选择RAW格式存储。

地表模型中,高程原始数据为DEM(Digital Elevation Model)格式。首先提取原始数据中1 km分辨率的点要素,点要素的x,y,z坐标分别对应经纬度及其高程值。 以北西—南东对角线将1 km边长的正方形划分为2个三角形,再使用VTK格式中三维非结构网格的类型6(三角形网格)形式将高程数据点连接为面,以面要素描述地形数据。

断层数据以数值化所收集的断层地表出露线作为基础线要素,根据已知断层倾向和倾角数据向深部延伸,通过三维非结构网格之类型6(三角形网格)生成断层面。部分文献^[12,14]中给出了一些断层的复杂模型,如龙门山断裂带的叠瓦状断层,但总体上,断层在三维空间的数据十分缺乏,断层深部的倾角变化鲜有准确资料,因此多数断层模型的建立依据地表出露倾角数值直线延伸到6~10 km深度上。仅龙门山断裂带依据已有信息,将断裂以铲形形态延伸至20 km深度。断层数据尚需不断补充、修改和完善,本文模型主要给出处理方案和示意效果。

岩石圈地震波速度和密度的原始数据为反演得到的各个台站下方的数据,深度方向分辨率为2 km。由于各台站观测数据以地表作为0 km,垂直向下延伸至100 km处;而不同台站之间高程差异可达3 km,因此不同观测点的相同深度值并不在同一水平面上。故首先对数据进行深度校正。校正方案为:

(1) $D=D_{0}-H$

式中:D为校正后深度值,D₀为校正前深度值,H为该台站地形数据中的高程值。校正后将每个台站下方数据重新以海平面为基准进行2 km为取样间隔的线性插值。在此基础上,对每一深度的数据逐层进行以1 km×1 km为取样间隔的双线性插值,获得模型内每个1 km×1 km×2 km单位块体的数值。矩形模型中未被流动台站覆盖的范围内的数据赋予空值。对P波速度、S波速度和密度均进行上述处理,并以RAW格式保存。

数据处理过程独立于后期可视化过程,一般需要根据初始数据状态确定具体处理方案。

模型的观察重点应用了体绘制技术(Volume rendering)。体绘制是一套将三维数据在二维屏幕上投影显示的技术。对于三维空间中一个由均匀网格组成的规则数据场,对象的属性信息(如密度,速度等)由体像素上相应的值描述。体绘制以体像素为基本操作单位计算每个体像素对投影图像的影响。技术需求包括:①对空间模型定义一个假想的光源及视角;②定义每个像素的颜色及透明度,通常用RGBA(Red,Green,Blue,Alpha)来定义每一个体像素的显示方式,这个定义也称为传递函数(Transfer Function)。传递函数可以是一个简单的斜坡函数,分段线性函数或是任意表格。体绘制技术通过传递函数确定每个体像素的色彩及透明度并投影到二维屏幕上,产生逼真的立体视觉效果。

体绘制技术已经广泛应用于医学CT或MRI (Magnetic Resonance Imaging)扫描图像数据处理(https://en.wikipedia.org/wiki/Scientific_visualization),形成的图像一般具有透明度,颜色普遍为人工指定的伪彩色。通过调节色彩及透明度,探索物体内部结构,突出显示感兴趣的区域,是体绘制技术的最大优越性。相应的代价是计算量大,对计算机硬件要求较高。

本文首次将这种技术应用于深部地球物理数据的可视化。

首先对经过高程校正并插值后数据的正确性进行验证。 图1a是通过Paraview显示的三维结构模型,与参考文献[3]中图4a完全对应,同样将模型从深度约24 km处分割为上下2个部分。图1a中断层数据为区内主要活动断裂在地表和24 km深度面上的投影。图1a模型外表面上S波速度分布与参考文献[3]几乎完全一致,证明数据处理和读取都正确。由于没有使用体绘制功能,模型内部的速度结构难以直接观察。

图1b为S波速度结构叠加断层数据。为了更直观地突出模型内部信息,图1b采用了与图1a不同的色标,并调节了透明度;同时附加断层信息。较浅深度上黄色低速体分布及延伸深度清晰可见。但深部信息被累加的半透明效果遮挡。这一问题在图像以二维打印方式呈现时尤为明显;而模型在屏幕上显示时,可以通过图形的旋转和翻滚,从各个角度看到结构体内部一定深度的信息,从而获得完整的三维结构认识。因此,科学数据可视化多以动图的方式展示,结合高性能可视化设备和技术(如超大屏幕、环形屏幕、立体视图等),能够获得绝佳的视觉效果。

图1c是在图1b基础上叠加了地表地形,给出综合了深部结构数据、断层数据和地表地形的三维模型。图1b,c中六面体所围的研究区内有淡蓝色雾,这是台站覆盖区域以外我们给定空值的区间。空值区在Paraview中以最低值显示(图1b图中速度比例尺最低值)。没有对地形面进行透明度调节,所以此时模型内部结构需要从其他角度观察。图1c的动画文件可以通过手机扫描右下角二维码下载。

图2单独对速度结构模型进行体绘制以突出显示模型西部20~30 km深部S波低速体以及深部岩石圈由东到西的速度结构变化。图2a显示研究区表层为大面积Vs≈3.1 km/s的速度层,这主要是浅部沉积层。将模型近地表10 km以上的部分切除可以更好地观察深部速度分布(图2b)。红色部分显示S波低速体主要分布于模型西部20~30 km深部范围,对应康定块体,同时在模型东南角也存在一个小区域范围S波低速异常。红色低速体整体被黄色速度区(Vs≈3.39 km/s)包围。东南部与西部2个主要黄色速度区之间的空隙大致对应龙门山断裂带和鲜水河—小江断裂带。图2c为模型北侧底部的仰视图,显示出100 km深度上,青藏高原下部岩石圈地幔速度值明显偏低。图2a~c采用相同色标。三维动态图可通过扫描色标旁所附二维码获取。

等值面(isosurface)是可视化中另一个重要手段。Paraview显示等值面可以通过输入阈值来提取出感兴趣的等值面,加以调整透明度将不需要的部分隐藏,可实现最佳显示效果。本文研究区内,四川盆地下莫霍面较清晰,西部高原区莫霍面形态复杂^[3]。理论上,莫霍面对应S波速度为3.85~4.48 km/s的跳跃层。通过提取研究区内Vs=4.0~4.2 km/s的等值面,可显示出大致与莫霍界面对应的波速跳跃层面(图3)。观察这一速度跳跃层面,显然模型东部该层面薄且分明,深度约30~40 km;而模型西部该层面较厚,深度基本上在50 km以下,且局部出现重复。注意所显示的等值面作为速度跳跃界面,仅是一定程度上代表莫霍面,并不合适由此确定莫霍面。准确确定莫霍面深度和形态需要更多资料和手段,超出了本文研究范围,不作进一步讨论。图3中下部附该模型的动态图二维码。动图可以更清楚地显示该层面的结构。

良好的体绘制和等值面视觉效果主要通过调节传递函数获得。寻找合适的传递函数也称为分类,通过恰当的色彩和透明度分类使图像突出显示感兴趣的部分并隐藏不感兴趣部分。下面以图2模型所采用的传递函数为例,给出调节方案及其效果的解释(图4)。图4窗口上部用于调整透明度,下部用于调整色标。色标中以3.3,3.5和4.3 km/s 3个值作为界线将模型分为红色、黄色、绿色与蓝紫色4个主要部分,对应图4中所标a~d区。红色区域a对应速度最低值区,适当调高不透明度使之便于观察;黄色区域b在模型中的分布基本包裹住红色区域,因此将其不透明度调至低于红色区域,使其被突出显示的同时不会完全挡住红色区域。绿色区域c对应地壳速度正常范围,本案例中以低的不透明度(高透明度)显示。蓝紫色区域d显示岩石圈地幔分布,将速度最高值区不透明度值调大,突出深部岩石圈地幔的S波速度高值区。水平色标条中可以根据需要任意增加颜色点并修改颜色,随后可以在上部区域调节透明度。色标的调整应考虑数据特征、欲突出显示的结构特征、色彩选择和匹配的美观等因素。一般占据较大空间区域的数据可认为属正常值,不建议采用高的不透明度;反之占据较小范围的数据(即异常值)可采用高的不透明度,使得异常特征得以突出显示。一定程度上,传递函数的选取和调节是图像美观或艺术性的决定因素。

通过不同体绘制处理后的速度结构模型,与浅部地质资料数据模型融合后能够更加直观地表征速度结构与浅部断层及地形之间的关系(图5)。图5a是10~100 km的模型。红色指示S波低速体,其边界与鲜水河断裂带及小金河断裂带分布基本一致,推测鲜水河断裂带及小金河断裂延伸至地下10 km以下;黄色指示S波速度较低部分,其北部边界线与龙门山断裂带分布基本一致,推测龙门山断裂带也延伸至10 km以下。图5b是24~100 km的模型。为了突出模型较深处Vs大于 4 km/s的部分, 对该模型提高了较浅部S波速度的透明度,可见在80~100 km深度上蓝色指示的S波速度分布存在明显的梯度变化界线,它与龙门山断裂带以及小金河断裂带走向基本一致。推断龙门山断裂带与小金河断裂带在深部可能是一个相连的超壳断裂,已有的地球物理数据也表明,龙门山断裂带的影响范围可能已达到地幔^[16]。

图4

传递函数实例
窗口上部的折线对应不透明度的值,下部颜色条对应色标调节

图5

通过体绘制处理并叠加断层的模型截图
(a) 切除了上部10km的S波速度结构模型俯视图;(b) 切除了上部24km的S波速度结构模型俯视图;这2种模型是由不同的颜色分布以及不同程度的体绘制处理而成

利用高性能可视化软件综合地质、地表及地球物理资料是一个新的尝试。综合多种信息资料并通过数字化方案进行表达是三维建模过程,而三维建模又是三维动力学模拟计算的基础。由于本项研究涉及的地球物理资料和以此为基础的动力学模拟计算的侧重点都在深部,且地球物理数据以地震波速度为主,更适合以体数据表达,因此选择了具有体绘制功能、更有利于观察深部结构的可视化软件。尝试结果表明,Paraview软件可以实现多种数据的综合并获得很好的可视化效果。

科学数据的三维可视化更有利于对深部结构的理解与更多信息的挖掘。通过体绘制和等值面功能及传递函数的调节,使三维图像能够充分展示速度结构数据中不同层位低速体的具体分布状态,显示地下空间高密度体与低密度体的分布位置,以及突出以特定速度值作为速度差异的分界线。在本文尝试中,也进一步推测龙门山断裂和小金河断裂的深部延伸状态。

由于地球物理探测技术的不断发展和人类探索地球的欲望,三维地下信息的获取变得越来越常规化,三维信息的有效显示成为必须考虑的问题。本文以青藏高原东缘高分辨率三维速度结构信息为例,给出了三维地球物理数据的可视化方案。该方案可用于各类地球物理数据,包括速度、密度和电、磁、热等参数。鉴于其他软件也具有相似的数据集成和可视化功能,相信诸如VisiIt和IDV等软件也可以实现类似数据的综合及可视化。希望本文提供的三维可视化方案有助于地球物理同行在未来学术交流过程中提供更美观、更直观、更易于理解并有助于分析的三维图像及动图展示。

致 谢:感谢中国地震局地质研究所刘启元老师向我们提供的数据以及在数据处理上的建议。同时感谢审稿人对论文提出的修改建议及意见。