基于曲面细分的地震数据三维显示方法
【技术领域】
[0001] 本发明属于计算机图形学和数据可视化领域，涉及一种基于曲面细分的地震数据 三维显示方法。
【背景技术】
[0002] 随着地震勘探技术的应用与发展，勘探家对地震解释成果的要求也越来越高。伴 随着计算机硬件的发展，尤其是GPU的快速发展，地震数据三维显示技术逐渐被人重视。传 统的地震数据显示方法:波形显示、灰度变密度显示以及彩色变密度显示都属于二维显示 方法，不能提供多个角度观察地震数据，往往该类地震数据解释系统会丢失很多细节。 [0003] 2008年Steven Lynch博士在他的博士论文《More than Meets the Eye_A Study in Seismic Visualization》提出了一种地震数据剖面三维显示方法，Steven Lynch博士 把该方法分为三个步骤：网格化、光照以及颜色显示。第一步的网格化是关键部分，由于GPU 植染单元为三角形，因此在地震数据网格化时需要指定相应的顶点索引。Steven Lynch博 士在确定地震数据的顶点索引时出现了如下问题，如图1所示，对于四个数据点组成的四边 形，有两种方式产生两个渲染三角形单元，但是对于含有位置信息的顶点来说，这两种连接 方式分别产生槽和脊，因此有时会对走势明显的地形产生破坏。如图2所示，由于槽的存在， 使得地形走势显示错误。
[0004] 1978年Edwin Catmull和Jim Clark提出经典的Catmull-Clark曲面细分算法，由 于该算法的递归特性不适用于GPU编程，基于三角形的Loop细分算法也不适用GPU编程。近 些年一大批适用于GPU编程的曲面细分算法被提出，比如：2008年CHARLES LOOP和SCOTT SCHAEFER两人于2008年提出了一种近似Catmull-Clark细分面的算法、同年Tamy Boubekeur提出的Phong细分算法。2009年微软发布Directxll，Directxll带来了细分曲面 管线，并且GPU在读取顶点缓存时可以以面片的方式读取，方便散乱点的曲面拟合以及曲面 细分算法的实现。
[0005] 地震数据采集时，是在地面上布置一条条侧线，按照一定的时间间隔进行采样，数 据是离散化保存，在进行三维解释的时候，近距离观察地质结构时会发现严重的网格化，如 图4所示。
【发明内容】
[0006] 发明目的：以解决地震数据网格化时的索引连接问题以及近距离观察时的网格化 问题。
[0007] 本发明包括:至少支持Directxll的GPU。
[0008] 本发明的特征在于，在计算机中是依次按以下步骤实现的：
[0009] 步骤(1)，SEG-Y格式地震数据的解析：
[0010] SEG-Y格式是地震数据最为普遍的格式之一。一般包括三部分，分别为EB⑶1C文件 头、二进制文件头以及地震道。用N表示数据体的第N道数据道，则第N道道头和道数据读取 位置分别为：
[0011] 3600+(N_l) (240+sample_num*data_size)
[0012] 3600+(Ν_1) (2400+sample_num*data_size)+240
[0013] SEG-Y结构是一样的，由于微机与和工作站的数据的存储格式不同，工作站高字节 在前，低字节在后，而微机则是低字节在前，高字节在后。在数据读取时需要进行高、低位数 据交换。
[0014] 步骤(2)，初始化顶点缓存：
[0015]自定义顶点结构，保存数据体的三维信息，从SEG-Y读取数据到二维数组中，记为 data[trace_num] [sample_num]。采用三阶贝塞尔曲面进行数据点的拟合，三阶贝塞尔曲面 的表达式为：
[0016]
[0017] 一个三阶贝塞尔曲面片的确定需要16个控制点，对于得到的二维数组数据源data
[trace_num] [sample_num]，GPU在进行读取时，是按照一维数组进行读取的，因此需要进行 二维数组到一维数组的转换，面片N与面片Μ左右相邻，控制点分别记为：
[0018] Ν[0][0]，Ν[1][0]，Ν[2][0]，Ν[3][0] Μ[0][0]，Μ[1][0]，Μ[2][0]，Μ[3][0]
[0019] Ν[0][1]，Ν[1][1]，Ν[2][1]，Ν[3][1] Μ[0][1]，Μ[1][1]，Μ[2][1]，Μ[3][1]
[0020] Ν[0][2]，Ν[1][2]，Ν[2][2]，Ν[3][2] Μ[0][2]，Μ[1][2]，Μ[2][2]，Μ[3][2]
[0021] Ν[0][3]，Ν[1][3]，Ν[2][3]，Ν[3][3] Μ[0][3]，Μ[1][3]，Μ[2][3]，Μ[3][3]
[0022] 为了消除Ν与Μ之间的裂缝，则需要Ν与Μ面片都包含相邻边的数据，因此Ν面片与Μ 面片消除裂缝的条件如下：
[0023] Μ[0][0]=Ν[3][0]，Μ[0][1]=Ν[3][1]，Μ[0][2]=Ν[3][2]，Μ[0][3]=Ν[3][3]
[0024] 面片的空间关系可以分为上下左右四个方位，消除裂缝算法如下：
[0025] For: i = l ,4,7---trace_num-3
[0026] For: j = l ,4,7---sample_num-3
[0027] Rowl data:data[i_l][j_l]，data[i+0][j_l]，data[i+l][j_l]，data[i+2][j_l]
[0028] Row2 data:data[i-l][j],data[i+0][j],data[i+l][j],data[i+2][j]
[0029] Row3 data:data[i-l][j+1],data[i+0][j+1],data[i+l][j+1],data[i+2][j+1]
[0030] Row4 data:data[i-l][j+2],data[i+0][j+2],data[i+l][j+2],data[i+2][j+2]
[0031] 步骤(3)，创建顶点布局：
[0032] 在自定义顶点结构中定义顶点的相关属性，此处只需要一个三维坐标信息，保存 地震数据点的位置信息。
[0033]步骤(4)，设置图元拓扑类型：
[0034] 用 1厶36七卩1'；[111；[1:;^61'0卩01087设置
[0035] D3Dll_PRmiTIVE_T0P0L0GY_16_C0NTR0L_P0INT_PATCHLIST，使得 GPU 以面片的方 式读取顶点缓存；
[0036] 为了能启用Directxll的细分阶段，设置Input Assembler的图元拓扑类型，使其 接收的是控制面片的结构，用IASetPrimitiveTopology设置
[0037] D3Dll_PRmiTIVE_T0P0L0GY_16_C0NTR0L_P0INT_PATCHLIST，使得 GPU 以面片的方 式读取顶点缓存。
[0038]步骤(5)，曲面细分：
[0039] 准备好顶点缓冲后，在GPU中实现曲面细分算法，GPU以
[0040] D3Dll_PRmiTIVE_T0P0L0GY_16_C0NTR0L_P0INT_PATCHLIST 组装方式对输入的顶 点缓存进行读取，对于给定的u，v，使用公式：
[0041]
[0042]求出细分得到的顶点的三维坐标。具体步骤是：
[0043] 第一步:设置一个静态缓冲来存储每一帧绘制时需要用到的数据，其中数据包括： 透视矩阵、摄像机位置、细分因子、地震数据最大幅值和地震数据最小幅值。
[0044] 第二步:Patch静态函数，程序对由16个控制点组成的面片执行一次，在该函数中 对每个面片的细分因子进行设置，为第一步中确定的细分因子。
[0045]第三步:hull着色程序，对每个输出控制点执行一次，设置细分domain为四边形， 输出控制点的数量为16，并指定pa t ch函数为第二步中的pa t ch函数。
[0046] 第四步：domain着色程序，在该程序中对tessellator阶段生成的顶点进行定位， 使用公式
[0047]
[0048] 侍到东疋Ui，v)生你厄相胆tf」」儿9、生你。
[0049] 步骤(6)，顶点颜色插值：
[0050] 在GPU的像素着色阶进行顶点颜色的确定，记步骤(5)第一步中得到的幅值最大值 最小值分别为MAX，MIN。对于给定顶点V的地震幅值Z，若其幅值为正，求Z与MAX比例，记为a， 则该顶点颜色为：（l_a，l_a，1)。若其幅值为负，求Z与MIN比例，记为b，则该顶点颜色为：（1， 1 -b, 1-b) 〇
[0051] 本发明使用Directxll编程工具，基于三阶贝塞尔曲面细分算法，实现了地震数据 的三维显示方法。对于准备好的顶点缓存，设置GPU以面片的方式组装控制点，因此不需要 索引缓存，从而省去了因确定索引缓冲而产生的槽和脊的问题。此外，基于GPU编程的曲面 细分技术的运用，生成了更多的控制点，使得网格的显示更加细腻平滑并具有很好的实时 性。
【附图说明】
[0052] 图1为四个采样点，采用不同的方式会产生槽和脊 [0053 ]图2为由于槽的存在，破坏了原有地形结构 [0054]图3为进行曲面细分后的效果图
[0055]图4为近距离观察时的效果图 [0056]图5为曲面细分过后显示效果图 [0057]图6为原始数据绘制得到的网格
[0058]图7为细分因子为3时绘制得到的网格
[0059]图8为基于曲面细分的地震数据三维显示方法流程图
[0060]图9为SEG-Y文件数据格式示意图
[0061]图10为左右相邻面片控制点不意图
【具体实施方式】
[0062]下面结合附图对本发明做进一步的说明。
[0063]图8是本发明的流程图，在计算机中按如下步骤进行：
[0064]步骤(1)，SEG-Y格式地震数据的解析：
[0065] SEG-Y格式是地震数据一般包括三部分，分别为EB⑶1C文件头、二进制文件头以及 地震道，如图9所示。用N表示数据体的第N道数据道，则第N道道头和道数据读取位置分别 为：
[0066] 3600+(N_l) (240+sample_num*data_size)
[0067] 3600+(Ν_1) (2400+sample_num*data_size)+240
[0068] SEG-Y结构是一样的，由于微机与和工作站的数据的存储格式不同，工作站高字节 在前，低字节在后，而微机则是低字节在前，高字节在后。在数据读取时需要进行高、低位数 据交换。
[0069]步骤(2)，初始化顶点缓存：
[0070]自定义顶点结构，保存数据体的三维信息，从SEG-Y读取数据到二维数组中。从二 维数组中初始化顶点缓冲。采用三阶贝塞尔曲面进行数据点的拟合。
[0071：
[0072] 对贝塞尔曲面求偏导数，得到偏导数可以很方便的求解顶点的法向量。得到法向 量后可以进行地震数据网格的光照显示。偏导数公式如下：
[0073]
[0074]
[0075] 一个三阶贝塞尔曲面片的确定需要16个控制点，对于得到的二维数组数据