南海是西太平洋一个独特的边缘海，保留了大陆边缘共轭张裂到海盆扩张的丰富信息(李家彪等，2011)，反射/折射地震方法可以得到可靠的地壳结构，直接为共轭张裂研究提供依据.南海大陆边缘地震测线的分布现状(图 1)为北部陆缘多，南部陆缘非常少(丘学林等，2011)，穿越南海南部礼乐滩的广角反射/折射地震剖面(OBS973-2)是南海南部陆缘仅有的几条深地震测线之一，其地壳结构特征在南海共轭扩张等研究中有重要意义(丘学林等，2011；阮爱国等，2011).

图 1

Fig. 1

图 1(a)南海深地震测线分布图，本文研究测线标记为红色(OBS973-2)和黄色(NH973-2)，黄色方框为图b在本图中的位置;(b)本文研究地震测线的台站布设图，白色圆圈代表成功回收的OBS，红色为丢失的OBSFig. 1(a)The distribution of deep seismic profiles in the South China Sea，the profiles of this study are marked in red(OBS973-2) and yellow(NH973-2)respectively. And the yellow rectangle is the location of Fig.b. And(b)，OBS dataset of this study，the white circles represent the OBS data available， and the red circles represent the lost OBS

阮爱国等(2011)对南部陆缘礼乐滩附近的OBS973-2 广角地震剖面进行了反演研究，并将其与北部陆缘的OBS2006-1剖面进行对比，进而探讨南海南北陆缘的共轭性及特征.但由于WARRPI软件的局限性(Ditmar and Makirs， 1996)，未对速度模型进行定量化误差分析，进而缺少了对速度模型可靠性等方面的评估.另外，研究区还存在一条与OBS973-2广角地震剖面共线的(相距2~3 km)的多道反射地震剖面(NH973-2)(丁巍伟和李家彪，2011)，多道反射地震和广角折射/反射地震是都常用于构造研究，但多道反射地震方法精于浅层而疏于深部构造(宋海斌等，2001；徐华宁等，2006；丁巍伟和李家彪，2011)，而广角折射/反射地震手段探测深度大，可直接应用于对大陆边缘深部结构、构造属性及张裂变形特征的研究(Yan et al., 2001；郝天珧和游庆瑜，2011；Zhao et al., 2010；阮爱国等， 2004，2011；吴振利等，2011；丘学林等，2011；黄海波等，2011；卫小冬等， 2011a，2011b).

基于以上考虑，本文使用正演迭代方法(Zelt and Smith， 1992)和自动反演方法(Hobro et al., 2003)重新对OBS973-2剖面进行了研究，同时，为了得到更精细的地壳结构信息，本文还对多道反射地震剖面(NH973-2)(丁巍伟和李家彪，2011)进行了时深转换研究.结合广角地震剖面和多道反射地震剖面，将反射剖面得到的界面信息作为折射剖面反演模型的约束，同时将折射剖面的反演速度用于反射剖面时深转换，通过相互约束，得到折射地震剖面的正演迭 代(Zelt and Smith， 1992)和自动反演(Hobro et al., 2003)速度模型，再结合多道地震深度剖面及其地质解释模型，探讨三个模型的差异及其缘由，最终得到礼乐滩及邻区的精细地壳结构，并对模型进行了定量化误差分析，以得到更可靠的深部动力学特征.

广角地震测线OBS973-2，是2009年由“实验2 号”科学考察船完成的.测线穿越礼乐滩东北部，向西北方向延伸进入东部次海盆，呈NW—SE走向，长369 km，共布设17台OBS，成功回收15台，阮爱国等(2011)已详细介绍了数据采集过程以及数据处理方法，本文不再赘述.分析表明，15个台站的数据质量良好，深部震相清晰(图 2)，所有台站识别出直达水波Pw和结晶基底以下上地壳的折射波Pg1(Watremez et al., 2011)，大部分台站识别出莫霍面的反射波PmP和上地幔的顶部的折射波Pn，少量台站识别出经过下地壳的折射波Pg2(Watremez et al., 2011)，最远震相可以追踪到120 km以外.另外，2009年“探宝号”调查船在该区域完成了一条多道地震剖面NH973-2，与OBS973-2测线基本重合，且最大间距仅为3 km，丁巍伟和李家彪(2011)对NH973-2剖面进行了处理和解释，可为本研究进行正演模拟提供浅部约束.

图 2

Fig. 2

图 2 OBS973-2测线 OBS3站位的地震记录(折合速度6 km·s-1)Fig. 2 Seismic data from OBS3 along OBS973-2 profile with reduced velocity of 6 km·s-1

首先，分析所有OBS记录剖面中的震相，根据多道地震剖面、区域地质资料以及前人在该区域所做的工作，建立地壳结构的初始模型.然后，采用2D射线追踪模拟计算各震相的理论走时曲线，并将该理论计算的走时与实际观测的走时进行对比，用试错法不断修改模型，使理论计算结果逐步向实测曲线逼近，获得一个较理想的模型，最后逐层对射线密集区域进行循环迭代计算，使得所有震相总的均方根走时残差最小，获得各台站的理论射线路径和二维地壳速度结构.采用的软件为RayInvr(Zelt and Smith， 1992).

本文根据OBS973-2测线各台站的震相特征，参考前人的结果(阮爱国等，2011)建立初始模型，对拾取的震相进行射线追踪计算.例如，OBS8台站位于COT区，剖面右侧震相偏移距大于60 km，识别出了Pg1和PmP震相，剖面左侧为海盆区，水深比右侧深，记录环境好，记录到的震相达到80 km以外，主要有Pg1、PmP和Pn震相(图 3a)，从获得的射线(图 3b)和震相拟合结果(图 3c)，可以看出结果很好.对测线上所有台站进行相互约束，共同拟合得到理想的速度结构模型.再对该模型进行迭代计算，得到最佳二维地壳结构和射线密度分布(图 4).

图 3

Fig. 3

图 3 OBS8台站地震记录(a)、射线追踪(b)和走时拟合情况(c)Fig. 3(a)Seismic section of vertical component of OBS 8 on profile OBS973-2，(b)simulation of ray-tracing，(c)The fitting of the calculated travel time to observed

图 4

Fig. 4

图 4 OBS973-2测线正演模型地壳速度结构(a)和射线密度分布图(b)(计算网格：0.5 km × 0.25 km)V.E.=5，图中黑色圆圈为OBS在剖面上的位置，其上数字为OBS编号.图a中白线为本文中展示的1D速度曲线的位置，白色数字为模型中的位置，单位为km.Fig. 4(a)Crustal velocity structure of 973-2 profile from forward model;(b)The distribution of rays density with the size of statistic network: 0.5 km×0.25 km V.E.=5. The black circles with numbers are OBS positions and its numbers. The white lines with numbers are the positions of 1D velocity curves discussed in this paper.

检验模型的可靠程度是通过在射线追踪过程中计算的走时残差以及 ²(实际观测走时与理论计算走时的拟合程度，越接近1表示拟合越好)来进行，同时还可以通过绘制射线密度覆盖图、计算最终模型的速度和界面扰动以及检测板测试等方法对数据的恢复能力和可靠性进行判断.本研究中，震相拾取时，根据数据质量的 好坏人为设定不确定性为±50~±80 ms. 结果表明各种震相走时残差均方根(RMS)均较小，随着深度略有增加，最大不超过120 ms. ²介于1.1到1.9之间，表明震相拟合较好(表 1).整条测线速度模型的射线密度分布(图 4b)表明，射线覆盖次数普遍大于5次，主要集中在10~40次，对整个模型有较好的覆盖，保证了模型的可靠性，因而模型结果有较好的约束和分辨率.最终模型的速度和界面扰动(Zelt and Smith， 1992；牛雄伟等，2014)(表 2)表明，地壳速度的不确定性很小(<0.25 km·s^-1)，莫霍面埋深不确定性为±0.29 km.

表 1(Table 1)

表 1 OBS973-2测线的震相拾取和模拟Table 1 Picks and rays simulation of profile OBS973-2 震相 参与计算的震相数 RMS(ms) χ2 Pw 513 77 1.185 Pg1 2632 112 1.830 Pg2 310 101 1.272 PmP 1014 120 1.762 Pn 460 95 1.085 总计 4929 110 1.643

表 1 OBS973-2测线的震相拾取和模拟Table 1 Picks and rays simulation of profile OBS973-2

表 2(Table 2)

表 2 速度模型的层速度和界面不确定性分析Table 2 Uncertainty analysis of the model parameters 模型参数 不确定性 上地壳顶部速度/(km·s-1) ±0.08 上地壳底部速度/(km·s-1) ±0.18 下地壳顶部速度/(km·s-1) ±0.23 下地壳底部速度/(km·s-1) ±0.22 莫霍面埋深/km ±0.29

表 2 速度模型的层速度和界面不确定性分析Table 2 Uncertainty analysis of the model parameters

自动反演采用Jive3D 软件(Hobro，1999; Hobro et al., 2003). 该软件适用于广角地震走时数据正演建模和层析成像反演，其算法使用规则网格节点来描述地壳速度和深度，允许数据拟合误差和模型复杂度同时降到最小，得到最小结构的1D、2D或3D模型，也支持折射、广角反射和大偏移距多道地震数据的同时反演(Hobro，1999).

其反演过程即使目标函数Ψ最小化的过程.目标函数Ψ定义为：

本文使用Jive3D软件进行2D建模反演，只需把3D模型中的某一水平维度的长度设置为1即可(Hobro，1999).反演采用改进了层剥法(Paulatto et al., 2010)，使得对各层和界面的约束随着深度增加更加连续有序.在每一步反演中，当前层及其上伏的层和界面都被反演.由于自动反演方法缺少直观性，故最好先使用正演方法对震相进行识别和确认(Paulatto et al., 2010).所以将正演迭代法获得的模型作为初始模型，包括海水层、地壳层和上地幔顶部(图 5a)，并使各层内速度场连续且光滑，设地 壳内速度从上向下在1.8～7.0 km·s^-1范围内变化，上地幔速度为8.0～8.2 km·s^-1. 首先根据网格节点，将初始模型自动插值形成均匀模型(图 5b)，并且在随后的每次反演开始时，使用2次B型样条插值得到沿深度方向等速度梯度的新速度网格和线性界面.为了尽可能避免过度拟合，使用尽量粗糙但又不影响拟合误差的网格节点(Scott et al., 2009).采用的网格间距为5 km(水平1)×0.5 km(水平2)×0.5 km(垂直).模型长370 km宽1 km深30 km.使用的速度节点在地壳层为76×3×62，在上地幔为76×3×9，海底面节点数为373(代表节点间距为 1 km)，莫霍面节点数为40(代表节点间距为10 km). 计算所得速度场的不连续可以被解释为界面，使用光滑和连续的多项式深度函数表示.

图 5

Fig. 5

图 5 自动反演方法过程(a)初始模型;(b)参数化模型;(c)最终模型.其他说明同图 4Fig. 5 Initial model(a)，parameterized model(b) and the final model(c)of inversion approach.The other introductions are the same as Fig. 4

自动反演过程不再包括直达水波，其他数据与正演迭代法相同，包括震相类型(Pg1和Pg2统一为地壳内折射震相Pg)和走时不确定性.反演过程中，模型光滑度参数λ_m 从0减小到-9.99(步长为-0.2，负值代表光滑度降低)，每个λ_m值迭代8次，直到得到稳定的模型，这时模型优化率从30%(λ_m 为0时)降低到0.001%(λ_m最小时)，²从884.01降低到1.11(图 6).拾取走时拟合率均大于80%，这样得到 的模型(图 5c)较能好地反映真实地壳结构(Paulatto et al., 2010).

图 6

Fig. 6

图 6 反演参数之间的关系图Fig. 6 Contour plots of the inversion parameters

使用检测板方法进行模型分辨率测试(Paulatto et al., 2010)(图 7).所加扰动为速度值的5%，扰动正弦函数半波长为50 km×2 km，结果表明，模型整体上有很好的纵向分辨率和较好的横向分辨率，边缘部分由于射线较少，分辨率较差.

图 7

Fig. 7

图 7 自动反演获得模型的检测板分辨率测试扰动半波长为50 km×2 km，扰动速度等值线为0.5 km·s-1.(a)理论模型;(b)恢复模型.其他说明同图 4.Fig. 7 Result from the checkerboard resolution test for inversion model Panel(a)shows perturbation added to the final model， and panel(b)shows recovered perturbation. A sinusoidal perturbation with half-wavelength of 50 km×2 km was tested. And the velocity perturbation contour is 0.5 km·s-1. The other introductions are the same as Fig. 4.

NH973-2测线数据，通过振幅补偿、静校正、增益处理、反褶积处理、去除多次波、速度分析、剩余静态校正及频谱分析等叠前处理后，再进行偏移叠加和叠后反褶积及高通滤波处理，获得高质量的偏移叠加剖面(丁巍伟和李家彪，2011).本研究在其工作的基础上，根据上述自动反演获得的速度分布，对多道地震剖面进行时深转换，并与折射地震结果进行对比验证，结果如图 8所示.可以看出，上述3种方法得到的结果互相一致，OBS剖面深部结构可靠性优于多道地震.

图 8

Fig. 8

图 8 NH973-2多道地震时深转换结果(黑线)，OBS973-2测线正演迭代获得的模型(红线)和自动反演获得的模型(绿线)Fig. 8 Comparison of results from time depth conversion profile of NH973-2(in black) and forward(in red)/inversion(in green)models of OBS973-2 profile

获得的正演迭代模型(图 4a)分为7层，最上面一层是海水层，速度为1.5 km·s^-1，沉积层有三 层.沿整个剖面，沉积层厚度较薄，局部地区沉积层 2和沉积层3缺失；整个沉积层的平均厚度为1~2 km，速度从浅部的1.8 km·s^-1向下增大到4.5 km·s^-1； 上地壳厚度变化较大，从洋盆的3~4 km变化到过 渡带的5~7 km，最后增厚至礼乐滩的7~8 km，速度从5.0~5.5 km·s^-1变化到6.4 km·s^-1；下地 壳厚度从洋盆处的3~4 km变化到过渡带的4~5 km，然后增至礼乐滩的9~12 km，速度从6.5 km·s^-1变化到6.9~7.2 km·s^-1；最后一层是上地幔顶部，速度从8.0 km·s^-1增至8.2 km·s^-1(表 3).如无特殊说明，本文讨论的地壳指沉积基底以下的地壳.与之相比，获得的反演模型(图 5c)共分3层，海水层、地壳层和上地幔，可以看出：各层厚度和正演迭代模型及多道地震深度剖面接近(表 3及图 8).但地壳厚度整体偏大约1.5 km，其主要原因为自动反演得到的莫霍面受PmP震相数量和走时的不确定性影响，但总体偏差小于模型厚度(30 km)的5%，表明正反演模型高度相似.在模型120 km 附近地壳中出现高速层，但检测板分析此处恢复较差.与正演方法相比，自动反演方法拟合好(²为1.11，小于正演模型的1.643)，而且效率高(自动计算，省时省力)，所以，在有足够边界条件约束的情况下，应当优先考虑使用自动反演方法.

表 3(Table 3)

表 3 正演迭代与自动反演结果对比(各层平均厚度/km)Table 3 Comparison of results from forward and inversion models(average layer depth/km) 沉积层 上地壳 下地壳 总厚度 正演模型 (阮爱国等，2011) 陆壳区 1.35 8.31 10.20 19.86 过渡区 1.83 5.71 4.80 12.34 洋壳区 1.14 3.10 3.00 7.24 正演模型 (RayInvr方法) 陆壳区 1.35 8.31 11.40 21.06 过渡区 1.83 5.71 4.41 11.95 洋壳区 1.14 3.10 3.26 7.50 反演模型 陆壳区 2.09 7.77 12.69 22.55 过渡区 1.16 6.88 5.63 13.67 洋壳区 1.41 3.45 3.92 8.78

表 3 正演迭代与自动反演结果对比(各层平均厚度/km)Table 3 Comparison of results from forward and inversion models(average layer depth/km)

多道地震时深转换模型(图 8)的海底面与OBS973-2正反演模型的海底面(浅部的黑色和红色实线)稍有差异是由于OBS973-2与NH973-2并不完全重合，两测线相交，距离最大处3 km，但并不影响对构造的解释.COT区域的沉积基底面与自动反演获得的模型的4.8 km·s^-1速度等值线一致，但都比正演迭代获得的模型深度大些，表明沉积层速度最大速度小于4.8 km·s^-1或者时深转换使用速度偏低.使用180~260 km区域内的平均地壳速度计算的195~260 km处的莫霍面与正反演模型都非常一致，进一步表明用自动反演获得的模型的速度结构可以约束得到可靠的多道地震剖面深度剖面.

与阮爱国等(2011)结果对比(表 3)，正演模型结果高度一致，但本文使用RayInvr软件正演得到的速度模型使用 ²进行统计误差分析(表 1)，并计算了模型参数的不确定性(表 2)，对模型的定量化误差分析使得本文结果更具有说服力.

模型中从南东向北西0~200 km为陆壳，平均 厚度为21~23 km，速度从5.5~5.8 km·s^-1变化到6.9~7.2 km·s^-1，上地壳存在横向速度异常和低速区.分别取正演迭代和自动反演获得的模型各自46 km、72 km和170 km处的一维速度结构进行比较(图 9)，可以看出与拉张陆缘速度结构有很好的一致性.在该区域的上地壳中具有速度横向异常，一维模型也存在速度倒转现象和低速区，而多道地震剖面上显示有一系列多米诺状半地堑(丁巍伟和李家彪，2011)，同时有部分断层甚至切穿到基底，甚至到下地壳(图 8)，表明该区域以伸展活动为主，并且其活动一直持至今(丁巍伟和李家彪，2011).

图 9

Fig. 9

图 9 陆壳区正反演结果1D速度曲线对比 (a)位于模型中46 km，红线为RayInvr正演模型，绿线为Jive3D反演模型，黑线区域代表拉张型陆壳速度范围(Christensen and Mooney， 1995)，阴影区域代表平均陆壳速度范围(Christensen and Mooney， 1995)，纵轴0 km处为沉积基底，图中的地壳分层界面使用RayInvr正演模型中的层界面；(b)位于模型中72 km，其他图例同图 9a；(c)位于模型中170 km，其他图例同图 9a.Fig. 9 Comparison of 1D vertical velocity of forward and inversion models in continental crust and other settings (a)Vertical velocity in model 46 km，red line from forward model， and green line from inversion model，black envelope represents extended continental crust(Christensen and Mooney， 1995) and shadow zone represents average continental crust(Christensen and Mooney， 1995). Ordinate 0 km means sediment base and interfaces are from RayInvr forward model;(b)Vertical velocity in model 72 km，the others are same as Fig. 9a.(c)Vertical velocity in model 170 km，the others are same as Fig. 9a.

模型中200~280 km划分为COT，速度结构揭示该区域地壳结构减薄，是介于张裂减薄的陆壳和洋壳之间的一种状态.正演迭代模型揭示该区域地壳厚度为11~12 km(自动反演模型中地壳厚度为 13~14 km).图 10展示了217 km、240 km和280 km 处的1D速度结构，并与ODP1277井蛇纹岩速度曲线(Christensen and Mooney， 1995)和大西洋COT区域的速度范围(Sibuet and Tucholke， 2012)进行了对比.结果显示，上地壳速度基本与大西洋COT区域速度及ODP1277井蛇纹岩速度一致，下地壳速度都比蛇纹岩速度低，下地壳没有发现高速层，说明南海南部陆缘为非火山属性.

图 10

Fig. 10

图 10 洋陆壳过渡区(COT)正反演结果1D速度曲线对比 (a)位于模型中217 km(红线和绿线)，黑线代表ODP1277井蛇纹岩速度曲线(Sibuet and Tucholke， 2012)，阴影区域代表COT区域的速度范围(Sibuet and Tucholke， 2012)，其他同图 9；(b)位于模型中240 km，其他图例同图 10a；(c)位于模型中280 km，其他图例同图 10a.Fig. 10 Comparison of 1D vertical velocity of forward and inversion models in COT and other settings (a)Vertical velocity in model 72 km(red line and green line)，black line represents serpentinisation velocity from ODP1277 site(Sibuet and Tucholke， 2012) and shadow zone represents COT from other settings(Sibuet and Tucholke， 2012). The others are same as Fig. 9;(b)Vertical velocity in model 240 km，the others are same as Fig. 10a;(c)Vertical velocity in model 280 km，the others are same as Fig. 10a.

模型中280~370 km为洋壳，正演迭代模型揭示该区域地壳厚度为7~8 km(自动反演模型中地壳厚度为8~9 km).选取测线292 km隆起区的1D速度曲线与大西洋0~127Ma洋壳速度曲线(White et al., 1992)对比(图 11a)，隆起区下地壳速度偏低，约为6.8 km·s^-1，表明该隆起可能为陆壳(Christensen and Mooney， 1995)，是礼乐滩地块随扩张脊向南移动时的残留，其连线为北西向，进一步验证了礼乐滩地块与其西北部的中沙地块共轭扩张(阮爱国等，2011；李家彪等，2011).另外测线325 km处(图 11b)的速度结构与大西洋0~127 Ma洋壳基本一致，表明本区域洋壳生成的外部条件(扩张速率/地幔温度/熔融物供给量等)可能与大西洋有较大相似性.

图 11

Fig. 11

图 11 洋壳区正反演结果1D速度曲线对比 (a)位于模型中292 km(红线和绿线)阴影区域代表大西洋0~127 Ma洋壳速度范围(White et al., 1992)，其他同图 9；(b)位于模型中325 km，其他图例同图 11a.Fig. 11 Comparison of 1D vertical velocity of forward and inversion models in oceanic crust and other settings (a)Vertical velocity in model 292 km(red line and green line)， and shadow zone represents Atlantic Ocean 0~127 Ma oceanic crust from other settings(White et al., 1992). The others are same as Fig. 9;(b)Vertical velocity in model 325 km，the others are same as Fig. 11a

本文以OBS973-2广角地震测线为素材，在建模过程中分别采用了正演迭代、自动反演计算，并将反演结果用于与之基本重合的多道地震剖面NH973-2时深转换，得到南海南部陆缘礼乐滩及邻区精细地壳结构.主要新认识是：(1)多道地震和水深数据对于广角地震的建模有重要意义.正演-试错迭代法可用于构筑模型基本格架，其特点是直观，易于掌握，便于人为控制主要结构参数，但精细度和光滑度稍差.以正演迭代获得的模型作为初始模型应用于自动反演方法，可以进一步对模型进行细化，得到更光滑更自然的结构，同时提高模型分辨率和精细度.反演获得的模型用于多道地震剖面的时深转化是可行的，结果也是好的.(2)礼乐滩陆缘及附近海盆的地壳结构特点是①沉积层较薄，平均厚度为 1~2 km，速度从1.8 km·s^-1向下递增到4.5~ 4.8 km·s^-1；②陆壳(0~200 km)平均厚度为 21~23 km，速度从5.5~5.8 km·s^-1变化到6.9~ 7.2 km·s^-1，上地壳存在横向速度异常和低速区，对应的多道地震剖面显示存在大量10~20 km的正断层，表明该区域拉张作用强烈；③洋陆过渡带(200~280 km)平均厚度为11~18 km，速度从5.0~5.7 km·s^-1 变化到6.8~7.0 km·s^-1，莫霍面埋深大约12 km，1D速度显示下地壳速度低于蛇纹岩速度，无高速层，表明该区域为非火山型陆缘；④洋壳(280~370 km)平均厚度为7~9 km，速度从4.8~5.8 km·s^-1变化到6.8~7.1 km·s^-1，1D速度显示地形隆起区下地壳速度偏低，莫霍面埋深增大，隆起区可能存在陆壳残留.