2. 中国石油大学(北京)油气资源与探测国家重点实验室 北京 102249;
3. 中国海洋石油研究中心 北京 100027
裂陷盆地的伸展量是构造物理的一个重要参数,也是盆地分析中的常用量。目前,最常用的伸展量计算方法是“地壳厚度变化计算法”(Allen and Allen,1992)和“断层平衡复原法”(Gibbs,1983; Davison,1986)。这两种方法计算结果的差异,国内外有较多研究(Pinet et al.,1987; Artyushkov et al.,1991; Bott,1992; 任建业等,1997; 刘绍文等,2001; Ziegler and Cloetingh,2004; Zhou et al.,2012)。“地壳厚度变化计算法”的结果较接近实际,“断层平衡复原法”的结果明显偏小。但“断层平衡复原法” 因方法简单而经常使用(刘绍文等,2001)。对于造成偏差的原因,以往人们多归因于犁式断层上盘的变形、 裂陷作用期间下地壳物质的损失、 上-下地壳变形机制的差异或计算误差等因素(任建业等,1997),也有作者将偏差归结为盆内下盘上升导致裂谷肩的侵蚀作用产生了测量误差(Kusznir and Ziegler,1992)。国内外理论界大多认为地震波因为不能分辨小断层而造成的累计伸展量是更重要的原因(Walsh et al.,1991)。 本文试图从另一个角度研究造成这两类伸展量差异的原因。
1 盆地伸展量计算的技术脉络吸收盆地变形的一种重要形式是地层的断裂,另一种重要形式是地层的韧性变形,但关于盆地韧性变形的量化计算还鲜有报道。地层韧性变形是连续的,即使是在断裂迅速发展的断陷期,韧性变形依然存在,而且还可能远大于裂后拗陷期。
“断层平衡复原法”是盆地伸展量计算最常用的方法,剖面数据直观,计算方法简单,数据容易获得,能够体现盆地内各个单元的伸展量差异,适用于有剖面资料的盆地研究。而“地壳厚度变化计算法”难以准确量化,因为它涉及岩石圈、 岩石密度、 压实校正、 热流等多参数,最大的问题是它不能体现盆地内各个单元的伸展量差异,而只适用于盆地尺度总的拉伸因子。“断层平衡复原法”主要针对某个断层的恢复,由于其它因素影响,盆地总的伸展量难以精确获得,“地壳厚度变化计算法”则恰恰相反。在过去,这两种方法的分别应用往往产生差异,使得盆地伸展量的量化计算结果争议颇多。本文认为,这两种方法的巨大差异是由于盆地韧性伸展所致,并在后文验证。由正断层滑动导致反射界面的错开进而导致地层伸展,称为“(断层)滑移伸展”,对应的伸展量为“(断层)滑移伸展量”,由滑移伸展量计算的伸展速率为“(断层)滑移伸展速率”; 由地层本身韧性变形导致的地层伸展称为“韧性伸展”,对应的伸展量为“韧性伸展量”,由韧性伸展量计算的伸展速率为“韧性伸展速率”。相对于韧性伸展来讲,滑移伸展属于岩石脆性伸展。
假设盆地伸展后的剖面(两个钉线间的水平距离)长度为L,盆地伸展前地层界限的长度为L0,地层界限序列k对应的盆地伸展前地层界限的长度为L0k,k对应的盆地伸展后地层界限的长度为Wk,如图 1所示。
Wk和L0k的差值表示为ak,Wk和L的差值表示为bk:
${a_k} = {W_k} - {L_{0k}}$ | (1) |
${b_k} = L - {W_k}$ | (2) |
其中,ak为反射界面k对应的韧性伸展量; bk为反射界面k对应的断层滑移伸展量。
2 盆地伸展量模型关于“断层平衡复原法”,文献(Gibbs,1983; Davison,1986)中已有表述与讨论。首先,假设条件如下: 1)水深约等于0,适用于陆相湖泊相、 滨浅海相的盆地; 2)无外来地体(如火山岩)侵入; 3)剖面边部某处设为“钉线”,“钉线”处不存在明显的水平层间差异滑动,以此来保障地体平衡; 4)不存在构造挤压反转。在应用“断层平衡复原法”时,首先要获得深度剖面,并据此计算目标剖面的相关指标。在深度剖面中,各个反射界面的参数表示如图 2所示。
如图 2所示,假设存在m个地层,基底为第m层,反射层k被断裂所分割,形成反射层单元长度分别为lki,角标k表示反射界面的垂向编号,介于0~m之间,角标i表示反射界面水平编号,lki在水平方向上投影长度为lvpki,产生断裂滑移,水平的反射层滑开的水平距离为δki,k+1层对应的水平反射层滑开的水平距离为δ(k+1)i,对于目前长度为L的剖面,假设它在拉张之初的原始长度为L0,则反射层k对应的水平投影总长度和反射层滑开总距离之和为目前研究区长度L,即:
$\sum\limits_i^{1 - n} {{l_{vpki}}} + \sum\limits_i^{1 - n} {{\delta _{ki}}} = L$ | (3) |
$\sum\limits_i^{1 - n} {{l_{ki}}} = {W_k}$ | (4) |
$\sum\limits_i^{1 - n} {{l_{ki}}} + \sum\limits_i^{1 - n} {{\delta _{ki}}} > L$ | (5) |
在韧性变形的条件下,如果不存在构造的挤压反转,那么在盆地伸展过程中,地层只发生水平伸展变长,根据图 2,则理论上有:
${L_{0k}} < \sum\limits_i^{1 - n} {{l_{ki}}} = {W_k} < L$ | (6) |
盆地伸展前目标地层的原始长度L0k如何确定?首先要知道海岸处的理论地壳厚度(理论Moho面深度)Hm,一般可以按照Hm大约为28~30 km计算,图 2中的实际研究剖面每个反射界面k和k+1界面之间地层的单层面积为
${S_k} = L \cdot {d_k}$ | (7) |
dk为反射界面k和k+1界面之间地层的单层当量厚度。由公式(7),可以依据Sk反求dk。
反射界面k到Moho面之间的总的当量地层厚度为hk:
${h_k} = \sum\limits_{k = m}^k {{d_k} = {H_0} - \sum\limits_{k = 0}^k {{d_k}} } $ | (8) |
其中,基底层对应的第m个单层厚度dk,也就是dm(d(k=m)),要由Moho面的深度(H0)来计算:
${d_m} = {H_0} - \sum\limits_{k = 0}^{m - 1} {{d_k}} $ | (9) |
需注意,公式(9)参与求和的是基底反射界面以上的地层,因此计算仅到达m-1编号的地层。Moho面目前深度为H0,H0有两个来源,一个是应用岩石圈研究获得的目前实际深度,另一个是应用地层、 密度数据和Airy重力均衡理论计算的理论Moho面深度(Allen and Allen,1992),这里不再赘述。鉴于压实和密度计算较为复杂,读者可参考孔隙度和深度的相关关系,如公式(10)的关系(张荣虎等,2011)进行密度压实计算:
$\phi = 0.40359{e^{ - 0.0002H}}$ | (10) |
式中,H为地层深度/m。密度压实的计算简易模型如下:
$\rho = 1000 \cdot \left[ {\phi + \left({1 - \phi } \right)\cdot 2.70} \right]$ | (11) |
式中,2.70为砂粒的相对密度; $\phi $为地层孔隙度,无因次; ρ为岩石密度/kg/m3。公式(10)和公式(11)提供了根据不同时期目标地层k对应的深度来计算目标地层密度的方法。公式(10)中系数0.403 59表示沉积岩一般的初始沉积孔隙度为40.359%,该模型假设每个时期新沉积的地层孔隙度(40%)符合石油地质领域的统计规律。
每个反射界面k形成时期tk对应时间的地层顶界应该在海面(湖面)附近,根据地质平衡剖面技术,下部地体在地震剖面上符合面积守恒原理,tk时期的实际厚度应该等于该时间理论Moho面深度Hm-tk,根据下式可计算出L0k:
${H_{m - tk}} \cdot {L_{0k}} = {h_k} \cdot L$ | (12) |
计算出L0k后,可以依据图 1和公式(1)及公式(2)计算ak和bk。那么,各个时期的运动量如何分离?一系列ak和bk均随着k的增大而增大,下层的变形含有其上每一层的累积变形,因此得出:
${a_{\left({k + 1} \right)}} - {a_k} = \Delta {a_k}$ | (13) |
式中,Δak为反射界面k对应的下部地层在此地层沉积期t(k+1)-tk之间总的韧性变形量。表示为韧性变形速率Tbpk:
${T_{bpk}} = \Delta {a_k}/\left[ {{t_{\left({k + 1} \right)}} - {t_k}} \right]$ | (14) |
反射界面k下部地层在沉积期t(k+1)-tk之间总的当量滑移变形量Δbk和反射层k对应的平均滑移变形速率Tbsk分别为:
${b_{\left({k + 1} \right)}} - {b_k} = \Delta {b_k}$ | (15) |
${T_{bsk}} = \Delta {b_k}/\left[ {{t_{\left({k + 1} \right)}} - {t_k}} \right]$ | (16) |
需要注意的是,因为地层数量比反射界面数量多一个,所以Sk、 hk、 ak、 bk均应该从S0、 h0、 a0、 b0开始计算。到此处,可以发现最下层ak、 bk对应值am、 bm分别为盆地总的韧性伸展量和总的断层滑移伸展量(断滑伸展量),两者是由不同的物理机理导致的。
盆地的平均伸展系数βb:
${\beta _b} = L/{L_{0m}}$ | (17) |
新方法计算的盆地总伸展率ebt为:
${e_{bt}} = {\beta _b} - 1$ | (18) |
盆地的平均伸展系数βb可以分成两部分,分别是断层滑移伸展和韧性变形伸展对应的伸展系数,断层滑移产生的拉张因子βbsm表示如下:
${\beta _{bsm}} = L/\left({{L_{0m}} + {a_m}} \right)$ | (19) |
传统方法计算的盆地伸展率ebst(断层滑移导致的盆地伸展率)为:
${e_{bst}} = {\beta _{bsm}} - 1$ | (20) |
韧性变形产生的拉张因子βbpm表示如下:
${\beta _{bpm}} = \left({{L_{0m}} + {a_m}} \right)/{L_{0m}}$ | (21) |
盆地韧性伸展率ebpt为:
${e_{bpt}} = {\beta _{bpm}} - 1$ | (22) |
对于某一个具体剖面,定义一个韧性伸展量和断层脆性滑移导致的伸展量之间的比值,称为总韧滑比λt,则:
${\lambda _t} = {a_m}/{b_m}$ | (23) |
上文为计算(断层)滑移伸展量和韧性伸展量的主体理论,在实际反射界面k对应的目前线段长度Wk统计工作中,盖层发展阶段经常会出现Wk大于剖面长度L的情况(颜丹平等,2003),因为在拗陷盖层发展阶段很少有大尺度的断裂滑动,由“连接两点的线段中,直线段长度最短”的几何基本原理可推知: 在没有断裂的情况下,作为非直线的Wk必然大于直线Lk,Lk为Wk在地表接受沉积初形成时对应的空间跨度,可见一般断陷盆地因持续伸展使Lk小于现今剖面长度L。这是一个很有意义的现象,使我们不得不思考以下问题。
1)首先涉及Wk统计工作中采用的精度问题,或者选择合适的Wk曲线粗糙度的问题。 如果应用非常细致的地震剖面,会发现实际的反射界面是粗糙不平的,现代海岸地貌和湖岸也是粗糙不平的,浅水区域的沉积大多如此,这是盆地沉积时的自然禀赋,在反射界面形成之初就具有的粗糙度影响了反射界面形成之初的Wk对应的长度Wko,此时的原始线长大于地理空间的跨度Lk。假设对地表沉积面某个识别精度i对应的沉积时间原始线段长度Wko满足下式:
${W_{ko}} = {\mu _i} \cdot {L_k}$ | (24) |
式中,μi是随着识别精度变化而变化的系数。对地表沉积面的识别精度越高,所得到的Wko值越是比Lk大得多,μi系数就越大。这说明,地震解释精确到一定小的尺度后,增加的只是这个曲线Wko的“粗糙度”而已,这就产生了第2个问题。
2)地震波不能分辨的小断层是影响伸展量计算差异的原因吗?地震波不能分辨的小断层如果被高效地识别出来,其实践意义也不过相当于增加了曲线Wko的“粗糙度”而已,识别出越多的微小断层,Wko的“粗糙度”越高,Wko也就比对应Lk大更多。在盆地断陷期地层Wk统计工作中,如果能够识别出越多的微小断层,那么将导致Wko增大,计算出的伸展量反而减小,因此,地震波不能分辨的小断层的累计伸展量并不是“断层平衡复原法”计算的伸展量小于“地壳厚度变化计算法”计算的累计伸展量的原因,是韧性变形承担的伸展量导致“断层平衡复原法”计算的伸展量小于“地壳厚度变化计算法”计算的累计伸展量。亚分辨正断层、 微破裂的确造成了小量的伸展,但是过于精细的地震解释增加的Wko曲线粗糙度导致的线长增加却可能远大于亚分辨正断层、 微破裂导致的伸展量。
3)既然盖层内的Wk比Lk大,那么Wk和现今剖面长度L的倒置关系还有意义吗?笔者认为不能否定“断层平衡复原法”计算伸展量的理论意义和实用价值。Wk小于现今剖面长度L的部分,的确是由于断裂的滑移伸展导致的,韧性伸展的存在在理论上也很容易理解,在反射界面形成之初就具有的粗糙度的确影响了Wko。只要还原出地表沉积面对应的沉积时间原始线段长度Wko,就可以计算出真实的滑移伸展量。
在实际的地层计算中,可以计算出盖层中诸多Wk中的最大值Wkmax,如果Wkmax大于现今剖面长度L,那么可以近似地认为Wkmax就是地表沉积面对应时间原始线段长度Wko,虽然会有小量误差,但整体还是适用的。在计算中获得了Wkmax,用Wkmax替换公式(24)中的Wko,并去掉μi中的对应识别精度符号i,则获得如下公式:
${W_{k\max }} = \mu \cdot L$ | (25) |
式中,μ是调整系数,大于等于1。实际地震剖面中首先获得Wk序列,然后依据下式计算出修正的数据序列Wkm:
${W_k}/\mu = {W_{km}}$ | (26) |
${W_{km}} - L = {E_k}$ | (27) |
式中,Ek为反射界面k对应的总的断层滑移伸展量,也就是传统的剖面中反射截面线长累加所能得到的伸展量。在统计Wk序列时,尽量保证各个反射界面的解释精度保持一致,曲线的粗糙度相差不多,计算的结果才准确。以上就完成了本文的理论主体的计算。
4 模型应用 4.1 伸展总量指标的应用将本文定量模型用于渤海湾盆地的计算,计算结果如表 1所示。表 1中,除歧口凹陷数据来源于Zhou et al.(2012)外,其它凹陷剖面数据均来源于陆克政等(1997)。在具体剖面的计算中,首先要选取“钉线”作为计算剖面的始末点,始末点间的距离为剖面长度L。表 1中的剖面长数据即为图 2中的剖面长度L,层长即为始末点间反射界面(不含断层面部分)的总长度${W_k} = \sum\limits_i^{1 - n} {{l_{ki}}} $。表 1说明,渤海湾盆地总伸展量指标显示韧性伸展量占渤海湾总伸展量的40%左右,与“断层平衡复原法”和“地壳厚度变化计算法”计算的巨大差异一致,韧性变形伸展量需要引起足够重视。渤海湾盆地是以断层滑移(脆性伸展)为主要伸展模式的盆地,地壳脆性相对较强,韧性/脆性伸展比(表 1 中的总韧滑比)平均为0.655,整体来说差异不大。
将上述理论用于渤海湾歧口凹陷某剖面(剖面来源: Zhou et al.,2012),经速度校正后进行伸展量的量化计算。剖面如图 3所示,反射层和地层的相关数据如表 2所示。
表 2中,断陷发展初期的理论L0m=41.63 km,这就是目标研究尺度L对应的伸展前原始长度,反射界面线长Wm=47.885 km,这是“断层平衡复原法”的内在缺陷导致的,而恰恰是这种差异保证了本文耦合模型存在的意义。02反射界面下部的Ed地层已经逐渐转化为断拗发展期,其上部的Ng地层的反射界面线长W1=55.652 km,小于W2=56.359 km,这是由于此时期(拗陷阶段初期)的“削峰填谷”结束造成的,本文作为特殊数据予以处理,不在计算分析之列。
表 2中,渤海湾歧口凹陷计算的韧性变形伸展量为6.255 km,而断层滑移伸展量为9.275 km,韧性变形伸展量大约为总伸展量的40%,这与一般发现的“断层平衡复原法”和“地壳厚度变化计算法”的巨大差异是相一致的。这个40%的比例可能也有普遍意义,诚然,不同盆地或一个盆地内的不同凹陷的韧性/断滑伸展量的比例差别应该是较大的。
渤海湾歧口凹陷不同盆地发展期的韧性伸展速率和断层滑移伸展速率分别如图 4所示。从图 4可以得出结论,断陷期的韧性伸展速率明显大于拗陷期的韧性伸展速率,这对于我们理解韧性伸展与应力的关系很有意义。岩石在差应力作用下发生塑性变形是岩石的固有性质,本文认为断陷期更低的水平应力使得驱动韧性变形的垂直—水平压力差更大,进而导致断陷期的韧性变形也相应的更大。
渤海湾歧口凹陷的韧性伸展速率和断层滑移伸展速率之间显示出某种指数式的正比关系(图 5),裂陷期不仅断层滑移剧烈,韧性变形速度也很大。由图 4对比可知,相对于断层滑移速率来说,韧性变形则要相对稳定很多。
渤海湾歧口凹陷韧性伸展速率的地质指导意义在于实现了韧性变形伸展和断层滑移伸展的定量分离,也印证了韧性变形伸展量有希望占到总伸展量的40%甚至更高,韧性变形伸展量需要引起足够重视。
5 讨 论由图 4对比还可以得到一个有意义的推论,就是断陷期1和断陷期2在剖面上显示为一个连续的断层发展期,但是断陷期2断层发育更小些,断层滑移伸展量上却非常明显小于断陷期1,笔者认为这说明了以下几点: 1)断层滑移伸展量和断层活动指数是两个不同的参数,断层活动指数是一个绝对的变化,断层滑移伸展量是由断层滑移保证的盆地伸展形式,更具有成因方面的暗示,其理论意义还有待深入的探讨; 2)微观上的断层滑移(断层发育指数)大并不都意味着伸展一定快,反过来快速的滑移伸展也未必一定表现为很大的断层发育指数;(3)断层的水平活动指数和垂直活动指数似乎代表了不同的伸展含义,水平活动指数更多指示的是伸展,而垂直活动指数更多指示的是区域的差异沉降。
由图 5的变形速率可知,在歧口凹陷中,Y=0.007 8e24.266x,由于公式并不是基于大量盆地的统计学拟合所得,因此并没有广泛的意义。但是,公式的指数形式可能有些意义,这似乎反映了韧性和断层滑移两者之间有某种关系,但需要慎重地认识和应用。图 5中的纵横坐标究竟谁决定谁,这是一个有意义的问题。
图 5的关系可以认为两个变形速度均是侧压系数的函数,在弹性力学体系中,侧压系数是泊松比的函数,可以认为是常数,而考虑到岩石的塑性变形时,岩石的侧压系数可以是变化的。从机理上讲,如果岩石能够实现通过韧性变形吸收应变,就不必发生很大的断层滑动,因此,笔者倾向于塑性变形能力决定断层的滑移量,反过来断层的滑移速度越大,形成的平衡侧压系数也越小,水平方向的塑性变形能力也就越强,塑性变形也就越快。
对于基底埋深在3 000 m以上的地层,水深一般不会超过200 m。因此中国除南海以外的东部地区的裂陷盆地本模型都可以应用,而西部的前陆盆地还需要专门的挤压韧性变形计算模型,对于存在严重挤压反转的盆地,还需要在本模型基础上做出提升和处理。因为统计本身的Wk也有一定的误差,但是L0k计算的误差主要是影响韧性和滑移的比例分配关系,而对于总的变形量,则作为一个系统误差而变化不大。在更多的盆地进行此类分离计算是有意义的,断裂滑移发展为主还是韧性变形伸展为主的临界值也很重要,可能是划分断陷期和拗陷期的依据之一。本模型应用于浅水盆地效果较好,而在形如南海这样的深水盆地,在做伸展恢复计算时就需要知道各个发展时期的古水深,古水深的研究难度一般较大,古水深数据的可靠性也较差,因此,深水盆地计算时需要进一步提升本文的模型。
6 结 论(1)本文的韧性伸展计算模型实现了韧性伸展量和断层滑移伸展量的定量分离,“断层平衡复原法”和“地壳厚度变化计算法”盆地伸展量计算的差距是岩石韧性变形所致,而非亚分辨断层所致。
(2)渤海湾盆地韧性伸展量约为总伸展量的40%,韧性变形伸展量需要引起足够重视。
(3)断陷期的韧性伸展速率明显大于拗陷期的韧性伸展速率。断层滑移速率和韧性变形速率呈某种指数式正比关系,裂陷期不仅断层滑移剧烈,而且韧性变形速率也很大。相对于断层滑移速率来说,韧性变形要相对稳定很多。
[1] | 刘绍文,王良书,刘波. 2001. 拉张盆地伸展量的分形分析--以渤海盆地为例.地质论评, 47(3):229-233. |
[2] | Liu Shaowen, Wang Liangshu and Liu Bo. 2001. Fractal analysis of extensional quantity of the extensional basin:Taking the Bohai Basin as an example. Geological Review, 47(3):229-233. |
[3] | 陆克政,漆家福,戴俊生等. 1997. 渤海湾新生代含油气盆地构造模式.北京:地质出版社. 1-251. |
[4] | Lu Kezheng, Qi Jiafu, Dai Junsheng et al. 1997. Tectonic Model of Cenozoic Petroliferous. Beijing:Geological Publishing House. 1-251. |
[5] | 任建业,张俊霞. 1997. 确定裂陷盆地伸展量的分数维法.地质科技情报, 16(3):102-106. |
[6] | Ren Jianye and Zhang Junxia. 1997. Fractal method for estimating extension amount of fault basin. Geological Science and Technology Information, 16(3):102-106. |
[7] | 颜丹平,田崇鲁,孟令波等. 2003. 伸展构造盆地的平衡剖面及其构造意义--以松辽盆地南部为例.地球科学, 28(3):275-280. |
[8] | Yan Danping, Tian Chonglu, Meng Lingbo et al. 2003. Balanced geological section for extensional tectonic basin and its implication:An example from southern Songliao Basin. Earth Science, 28(3):275-280. |
[9] | 张荣虎,姚根顺,寿建峰等. 2011. 沉积、成岩、构造一体化孔隙度预测模型.石油勘探与开发, 38(2):145-151. |
[10] | Zhang Ronghu, Yao Genshun, Shou Jianfeng et al. 2011. An integration porosity forecast model of deposition, diagenesis and structure. Petroleum Exploration and Development, 38(2):145-151. |
[11] | Allen P A and Allen J R. 1992. Basin Analysis:Principles and Applications. Oxford:Blackwell Scientifics Publications. 1-451. |
[12] | Artyushkov E V, Baer M A and Letnilov F A. 1991. On the mechanism of graben formation. Tectonophysics, 197(2):99-115. |
[13] | Bott M H P. 1992. Modelling the loading stresses associated with active continental rift systems. Tectonophysics, 215(1-2):99-115. |
[14] | Davison Ⅰ. 1986. Listric normal fault profiles:Calculation using bed-length balance and fault displacement. Journal of Structural Geology, 8(2):209-210. |
[15] | Gibbs A D. 1983. Balanced cross-section construction from seismic sections in areas of extensional tectonics. Journal of Structural Geology, 5(2):153-160. |
[16] | Kusznir N J and Ziegler P A. 1992. The mechanics of continental extension and sedimentary basin formation:A simple-shear/pure-shear flexural cantilever model. Tectonophysics, 215(2):117-131. |
[17] | Pinet B, Montadert L, Mascle A et al. 1987. New insights on the structure and the formation of sedimentary basins from deep seismic profiling in western Europe. In:Brooks J and Glennie K W(Eds.). Petroleum Geology of North West Europe. London:Graham & Trotman. 11-31. |
[18] | Walsh J, Watterson J and Yielding G et al. 1991. The importance of small-scale faulting in regional extension. Nature, 305:391-393. |
[19] | Zhou L H, Fu L X, Lou D et al. 2012. Structural anatomy and dynamics of evolution of the Qikou sag, Bohai Bay Basin:Implications for the destruction of North China craton. Journal of Asian Earth Sciences, 47:94-106. |
[20] | Ziegler P A and Cloetingh S. 2004. Dynamic processes controlling evolution of rifted basins. Earth-Science Reviews, 64(1):1-50. |
2. State Key Laboratory of Petroleum Resources and Prospecting, China University of Petroleum, Beijing 102249;
3. Research Center, CNOOC, Beijing 100027