谱反演方法在致密薄层砂体预测中的应用研究

迟唤昭,刘 财,单玄龙,鹿 琪

(1.吉林大学地球科学学院,吉林长春130061;2.吉林大学地球探测科学与技术学院,吉林长春130026)

谱反演方法在致密薄层砂体预测中的应用研究

迟唤昭1,刘 财2,单玄龙1,鹿 琪2

(1.吉林大学地球科学学院,吉林长春130061;2.吉林大学地球探测科学与技术学院,吉林长春130026)

在以致密砂岩为储层的非常规致密油勘探开发中,由于靶区选择和水平井部署的需要,落实致密砂体的横向展布至关重要。松辽盆地北部白垩系扶余油层致密砂体整体薄层发育,特别是T2强反射界面屏蔽问题的存在,一直制约着致密薄层砂体的识别精度。为此,开展了针对大庆长垣扶余油层致密薄砂体预测的谱反演方法应用研究。在简介谱反演方法原理的基础上,为了保证反演算法的计算精度,利用复数域快速匹配追踪算法提高时频转换的精度;通过提取时变子波提高反射系数时频谱的准确性。应用研究结果表明,将谱反演得到的地震反射系数体与宽频零相位子波进行褶积运算,获得的研究区宽频保幅地震数据有效削弱了T2的屏蔽效应,改善了复合波的叠置问题,清晰地反映出了薄层砂体顶、底界面,并很好地保持了振幅与砂体的对应关系。

谱反演;快速匹配追踪;时变子波;致密砂体;T2屏蔽效应

非常规油气资源已成为人们关注的勘探热点之一,其中以致密砂岩为储层的致密油/气是目前最现实的勘探目标[1]。针对致密油储层的精细勘探开发,国内多家油田开始了大规模的水平井部署,这对地震资料的分辨率和薄储层预测方法的应用效果提出了更高的要求。为了提高致密薄砂体地震预测的精度,人们相继提出了一系列技术方法,目前实际应用较广泛的是井控波阻抗反演技术,以及基于地震保幅资料的属性分析技术[2-4]。前者由于井资料的参与,具有分辨率高的优点,但是其缺点同样突出,主要表现在“牛眼”现象严重,井间内插算法无法适应河流相岩相的快速变化;而后者不受井的约束,虽然分辨率不是很高,但是对砂体横向厚度的变化刻画得更加真实,在致密油勘探中已发挥出较好的效果。实践表明,抛弃保幅资料仅依据有井约束反演资料作为水平井轨迹设计与导向,大多数情况下砂体的钻遇率均不尽如人意。以大庆长垣扶余油层为例,存在的地质问题是:3～5m的砂泥岩薄互层发育,古河道多期叠置现象普遍,区内地层岩性的横向非均质性非常强,相变速度很快[5]。如果想进一步提高砂体预测的符合率,并且完全满足水平井钻探的需要(至少在1000m以下成功预测3m厚度的砂岩),常规地震资料的分辨率已变得无能为力。尤其是T2强反射界面屏蔽现象的存在,致使目的层有效信号微弱,不足以清晰刻画薄砂体。因此,提高保幅资料的“真分辨率”乃是大庆长垣扶余油层致密薄砂体横向展布预测需要解决的核心问题。

谱反演是一种基于先验信息和频谱分解技术的高分辨率反射系数成像方法,可以提高小于调谐厚度薄地层的成像精度。Portniaguine等[6]首先提出了有别于常规地震反演理论的谱反演方法,其主要特点是只采用部分频谱资料就可以反演出稀疏反射系数或层厚。Puryear等[7]利用奇偶反射系数分解理论,发展了一种新的谱反演算法,通过对奇偶反射系数权重的重新匹配,利用短时傅里叶变换实现时频转换,理论上可以更稳健地识别小于地震调谐厚度的薄层。目前，国内有关谱反演方法理论研究和应用研究的文献尚为数不多。柴新涛等[8]利用LSQR算法对谱反演目标函数的求解进行了理论探讨;曹鉴华等[9]对谱反演处理技术进行了阐述,并将其应用到实际生产资料中,清晰识别了薄储层反射特征;刘万金等[10]用谱反演方法对三维数据体进行反演及地震属性分析,提高了识别沉积相和小断层的能力；孙雷鸣等[11]提出了线性谱反演约束方法,使反演结果更加精确,提高了地震勘探分辨率的能力。

针对非常规致密油薄储层精细预测的需求,本文在简要介绍谱反演方法原理的基础上,对影响谱反演计算精度的关键技术环节进行了改进:利用复数域快速匹配追踪算法提高了时频转换精度[12-13];通过提取时变子波提高反射系数时频谱的准确性[14-15];通过对谱反演反射系数体进行0.5ms的重新采样,提高地震资料与测井资料的匹配程度。大庆长垣扶余油层致密薄砂体预测的实际应用取得了良好效果。

1 谱反演方法原理

谱反演方法是一种在时频分析及子波提取技术的基础上,利用反射系数奇偶分解理论[16],求解薄层反射系数及其地层厚度的反演方法。在保幅的前提下,可将由于地层的调谐作用而模糊了的薄层信息反演出来,提高地震资料的分辨率。

1.1 多层模型应用

尽管双层模型[17-18]由于容易反演而很有价值,但是实际地层具有多套反射层,因此将该方法推广到多层反射模型更有实际意义。我们可以将反射系数序列r(t)表示为奇部分量和偶部分量脉冲对的和[7],即:

(1)

式中:ro(t)和re(t)分别表示奇部分量与偶部分量;T(t)为地层时间厚度;δe为偶部脉冲对;δo为奇部脉冲对;τ为时间延迟。

对反射系数进行时频转换可得:

(2)

式中:r(t,f)为r(t)的时频谱;tw为时窗的半长度。

i2ro(t)sin[πfT(t)]}dt

(3)

式中:f为频率。

多套地层模型中包括两个以上的反射同相轴,所以在计算的时候有必要使用正好适于层间反射界面的目标函数。如果子波时频谱已知,可以利用最小二乘共轭梯度法[19-21]求解目标函数o(t,re,ro,T)的最小值,进而得到r(t)和T(t),该函数为:

(4)

式中:fL为低频截止值;fH为高频截止值;αe和αo为权重函数,用来调节噪声和分辨率。在αe/αo为高值时,反射率就比较接近Widess模型,分辨率极限值变为λ/8。

1.2 关键技术环节的改进

从谱反演目标函数的构建上可以看出,地震信号与子波的时频转换以及子波的提取是直接影响谱反演算法本身精度的关键技术环节。

1.2.1 时频转换

为了保证地震信号与子波的时频谱的计算精度及时效性,本文利用复数域快速匹配追踪算法进行时频转换，在保证时频分辨率高的同时,确保了较高的计算效率,这已在实际应用中得到了验证。

具体流程如下。

1) 利用多参数控制的Morlet小波函数建立完备原子库,小波函数表达式为:

(5)

式中:μ为位置函数;ωm为频率;φ为相位;σ为尺度。

2) 利用Hilbert变换得到复地震道,进一步计算地震信号的三瞬属性。

3) 初步确定小波匹配原子:瞬时振幅确定匹配起始位置;瞬时相位替代小波相位;瞬时频率平均替代小波主频。之后利用局部动态扫描方法提高匹配计算的效率。

4) 通过迭代匹配,得到一组匹配最佳的时频原子,完成地震信号的稀疏分解,地震信号可表示为:

(6)

式中:an为匹配段的地震信号;m为时频原子;Rn+1为残差。

5)计算时频原子的Wigner分布之和,得到地震信号的时频谱(如图1,单分量信号不存在交叉干扰项,而且时频分辨率很高[22])。

(7)

式中:W(t,f)为Wigner分布;x*为x的共轭;τ为时间延迟。

图1给出了模拟信号(图1a)及其Wigner分布时频谱(图1b)与广义S变换时频谱(图1c)的时频分辨率比较结果。

1.2.2 子波提取

至于地震子波的提取,当然我们最需要的是既具有时变性质又具有空变性质的地震子波。子波的提取方法中基于高阶统计量理论的方法[23]目前研究较多,本文通过提取时变子波提高反射系数时频谱的准确性。具体做法是首先对地震数据开窗分段,在时窗内采用实用性更强的统计子波估计方法,实现时变子波求取(图2)。采用时变子波与前人采用非时变子波进行了对比,结果如图3所示(红色椭圆标记内),可以看出,提取时变子波作谱反演处理后的地震剖面分辨率更高。

图1 模拟信号及其时频谱

图2 时变子波及其频谱

图3 采用提取的非时变子波(a)及时变子波(b)作谱反演处理后的地震剖面

2 实际应用效果分析

大庆长垣YP1井区非常规致密油储层——扶余油层以曲流河多期叠置沉积为主。虽然该区岩性比较简单,以砂泥岩为主,但储层物性较差,属低孔、低渗储层(孔隙度低于8%的储层占45.8%,渗透率低于10-3μm2的储层占74.7%)。另外,该区河道改道迁移现象普遍存在,不但砂体错叠连片,而且与泥岩交互形成的薄互层(一般小于5m)发育。在这种特定的沉积环境及物性条件下,致使砂层整体较薄,地震资料上层间干涉严重,薄层叠置与单砂体地震反射特征相似,难以区分。亟需能够清晰刻画河道砂体边界的地震预测技术,以满足致密油薄储层精细勘探需求。为此,本文尝试将谱反演技术应用于研究区致密薄层砂体的精细预测,以期为大庆长垣非常规致密油的勘探开发提供技术支撑。

与其它叠后反演算法一样,谱反演的最终结果也是反射系数体,反射系数的极性与大小是砂体顶、底界面波阻抗差的反映。因此,利用谱反演技术获得的高分辨率反射系数体,可以对薄砂体的顶、底界面进行追踪。抽取F42-F43-F44-F45井连井测线的反射系数剖面(图4)与常规保幅地震剖面(图5)对比,可以看出,反射系数资料更好地指明了砂泥界面。在水平井随钻导向时,将谱反演的反射系数资料与上、下伽马曲线结合应用,将可有效防止水平井轨迹“出层”现象的发生。另外,我们对反射系数进行了0.5ms重新采样,其目的主要是提高地震资料与测井资料的匹配程度,虽然不能增加地震资料的新信息,却能提高反演结果的分辨率[24]。

图4 F42-F43-F44-F45井连井测线反射系数剖面

图5 F42-F43-F44-F45井连井测线常规保幅地震剖面

但是,谱反演作为一种无井反演方法,在提高分辨率方面并没有井的参与控制,因此其分辨率仍受限于保幅地震资料的分辨率。同时,保幅地震资料的信噪比也是影响谱反演效果的一个关键因素。

反射系数资料虽然能清晰反映砂岩顶、底界面,但这是建立在工区有足够井资料、对砂体地震响应特征充分认识的基础上的。反射系数资料归根结底只是反射界面波阻抗差的反映,而人们更关心的是工区内砂体横向上厚薄的展布情况。波阻抗是直接与岩性相关的信息,但是在谱反演反射系数资料的基础上,采用道积分方法得到的相对波阻抗资料分辨率很差,并不适用于3～5m薄砂体的识别。因为积分算法是一种“平滑”作用,平滑掉了谱反演出来的细节信息。

本文根据Chopra等[25]提出的方法,通过对谱反演获取的地震反射系数体与宽频零相位子波进行褶积运算,获得宽频保幅地震数据。将该方法应用到大庆长垣致密砂体预测中,取得了非常好的应用效果,为研究区非常规致密油薄储层的精细刻画提供了一种新的技术手段。

图6 谱反演宽频资料(a)和常规保幅资料(b)与合成记录的匹配对比

图6是谱反演宽频资料和常规保幅资料与合成记录的匹配对比情况。通过3个标记处(图中标注1,2,3)的对比可以看出,相对于常规保幅资料,谱反演宽频资料与合成记录的符合度更高,说明其在保真性、分辨率方面均好于常规保幅资料,对薄砂体有更强的识别能力。尤其在标记1处提高的分辨率更有意义:标记1处上方为T2地震反射层位(它是青山口组的底,扶余油层之上)。其上覆的青山口组一直被公认为下伏扶余油层的油源层,两者构成了“上生下储”的源储形式。从非常规油气勘探的思路可知,扶余油层上部层组更接近油源,其含油性应该最好。但是,一直以来由于T2强反射界面对下部地层的屏蔽效应,致使扶余上部油层组反射能量很弱,无法通过地震响应特征准确认识其分布规律。所以,扶余油层上部层组的勘探程度最低,这直接影响了对大庆长垣资源潜力的认识。谱反演宽频资料在一定程度上削弱了T2的屏蔽效应,使该部分储层的地震响应得到了加强。

基于测井曲线的合成记录只反映出“点”的信息量,评价谱反演宽频资料对横向砂体展布的刻画能力更为重要。图7显示的是常规保幅资料,图8显示的是谱反演宽频资料。常规保幅资料是大庆长垣扶余油层地震技术攻关的结果,其良好的保幅性已经在大庆长垣多个地震工区的水平井部署中得到验证。但是常规保幅资料也存在严重的T2屏蔽效应(图7上部绿色椭圆标记内)、地震波复合叠加(图7中红色椭圆和蓝色方框标记内)。而在图8所示的谱反演宽频资料的绿色椭圆标记内,薄砂体的反射能量增强,砂体的连通性加强,在一定程度上削弱了T2屏蔽效应;在图8中红色椭圆和蓝色方框标记内,复合波已被打开,而且分开的同相轴与砂体一一对应,这将更有利于五级层序的精细解释。图9a和图9b分别为常规保幅资料与谱反演宽频资料的频谱,对比可见宽频数据体具有更宽的频带范围,并不是简单地拓展高频,而是切实提高了地震资料的“真分辨率”,从而保证了薄砂体和厚砂体均具有相应的地震响应。

图7 P64-P41-P611-P611X1井连井剖面常规保幅资料

图8 P64-P41-P611-P611X1井连井剖面谱反演宽频资料

图9 图7常规保幅资料(a)和图8谱反演宽频资料(b)的频谱

图10 YP1井区扶余油层中部层位砂体预测结果

图10a和图10b分别为基于常规保幅资料和谱反演宽频资料的YP1井区扶余油层中部层位砂体预测结果。对比图10a和图10b可以看出,谱反演宽频资料在横向上与常规保幅资料保持了基本的一致性,这说明前者同样具有良好的保幅性;而在细节上宽频资料的河道砂体边界更清晰,尤其是对地质家更关心的地质体(“点砂坝”)的刻画能力更强,这说明其具有更高的横向分辨率。

综上所述,谱反演技术可以将由于地层的调谐作用而模糊了的薄层信息(反射系数)反演出来,其宽频资料在横向保幅的基础上提高了对砂体的识别能力及细分研究单元层位解释的精度。选取YP1井目的层段(河流相沉积)进行谱反演方法验证,钻前预测砂岩厚度1523m,根据随钻自然伽马测井曲线计算砂岩厚度1447m,砂体预测正确率达到95%。

3 结论

本文针对大庆长垣扶余油层致密薄砂体预测的谱反演方法应用研究取得如下认识。

1) 在介绍谱反演算法的基础上,本文利用复数域快速匹配追踪算法提高时频转换的精度;通过提取时变子波提高了反射系数时频谱的准确性。

2) 谱反演高分辨率反射系数体,可以对砂体顶、底界面进行准确追踪,这将大大提高非常规致密油勘探中水平井的钻探成功率。

3) 为了更好地发挥地震技术横向分辨率高的优点,利用谱反演高分辨率反射系数体与宽频零相位子波褶积,得到保幅的宽频地震资料。这对削弱大庆长垣扶余油层的T2屏蔽效应,克服地震波复合叠加,提高细分研究单元的解释精度有很大帮助。

[1] 邹才能,朱如凯,吴松涛,等.常规与非常规油气聚集类型、特征、机理及展望——以中国致密油和致密气为例[J].石油学报,2012,33(2):173-187 Zou C N,Zhu R K,Wu S T,et al.Types,characteristics,genesis and prospects of conventional and unconventional hydrocarbon accumulations:taking tight oil and tight gas in China as an instance[J].Acta Petrolei Sinica,2012,33(2):173-187

[2] 杨立强,邬长武,董宁.基于模拟退火算法的随机反演技术在砂体预测中的应用[J].地球物理学进展,2013,28(1):287-292 Yang L Q,Wu C W,Dong N.Application of stochastic inversion based on simulated annealing algorithm for sand body prediction[J].Progress in Geophysics,2013,28(1):287-292

[3] 沈财余,江洁,赵华,等.测井约束地震反演解决地质问题能力的探讨[J].石油地球物理勘探,2002,37(4):372-376 Shen C Y,Jiang J,Zhao H,et al.Approaching ability of using well-log constrained inversion to solve geologic problems[J].Oil Geophysical Prospecting,2002,37(4):372-376

[4] 杨占龙,彭立才,陈启林,等.地震属性分析与岩性油气藏勘探[J].石油物探,2007,46(2):131-136 Yang Z L,Peng L C,Chen Q L,et al.Seismic attributes analysis and lithological reservoir exploration[J].Geophysical Prospecting for Petroleum,2007,46(2):131-136

[5] 李雪英,文慧俭,陈树民,等.等厚薄互层时频特征的正演模拟[J].地球物理学报,2013,56(3):1033-1042 Li X Y,Wen H J,Chen S M,et al.Forward modeling studies on the time-frequency characteristics of isopachous thin interbedding[J].Chinese Journal of Geophysics,2013,56(3):1033-1042

[6] Portniaguine O,Castagna J P.Spectral inversion:lessons from modeling and Boonesville casse study[J].SEG Technical Program Expanded Abstracts,2005,24:1638-1641

[7] Puryear C I,Castagna J P.Layer-thickness determination and stratigraphic interpretation using spectral inversion:theory and application[J].Geophysics,2008,73(2):37-48

[8] 柴新涛,李振春,韩文功,等.基于LSQR算法的谱反演方法研究[J].石油物探,2012,51(1):11-18 Chai X T,Li Z C,Han W G,et al.Spectral inversion method analysis based on the LSQR algorithm[J].Geophysical Prospecting for Petroleum,2012,51(1):11-18

[9] 曹鉴华,邱智海,郭得海,等.叠后地震数据的谱反演处理技术及其应用浅析[J].地球物理学进展,2013,28(1):387-393 Cao J H,Qiu Z H,Guo D H,et al.Post-stack seismic spectral-inversion technique and its application[J].Progress in Geophysics,2013,28(1):387-393

[10] 刘万金,周辉,袁三一,等.谱反演在地震属性解释中的应用[J].石油地球物理勘探,2013,48(3):423-428 Liu W J,Zhou H,Yuan S Y,et al.Application of spectral inversion in the seismic attribute interpretation[J].Oil Geophysical Prospecting,2013,48(3):423-428

[11] 孙雷鸣,曾维辉,方中于.地震薄层反射系数谱反演算法研究及应用[J].物探化探计算技术,2014,36(4):462-470 Sun L M,Zeng W H,Fang Z Y.Thin-bed reflectivity inversion and seismic application[J].Computing Techniques for Geophysical and Geochemical Exploration,2014,36(4):462-470

[12] 武国宁,曹思远,孙娜.基于复数道地震记录的匹配追踪算法及其在储层预测中的应用[J].地球物理学报,2012,55(6):2027-2034 Wu G N,Cao S Y,Sun N.Matching pursuit method based on complex seismic traces and its application of hydrocarbon exploration[J].Chinese Journal of Geophysics,2012,55(6):2027-2034

[13] 张繁昌,李传辉.基于正交时频原子的地震信号快速匹配追踪[J].地球物理学报,2012,55(1):277-283 Zhang F C,Li C H.Orthogonal time-frequency atom based fast matching pursuit for seismic signal[J].Chinese Journal of Geophysics,2012,55(1):277-283

[14] 梁光河.地震子波的时变与短程微曲多次反射[J].石油物探,1997,36(3):18-27 Liang G H.The time-variant seismic wavelet and the peg-leg multiple[J].Geophysical Prospecting for Petroleum,1997,36(3):18-27

[15] 王晓华,秦义龙,黄真萍,等.基于频谱分析的时变子波反褶积[J].石油物探,1998,37(增刊):109-112 Wang X H,Qing Y L,Huang Z P,et al.Time-variant wavelet deconvolution on spectrum analysis[J].Geophysical Prospecting for Petroleum,1998,37(S1):109-112

[16] Lopez J A,Partyka G,Haskell N L,et al.Identification of deltaic facies with 3-D seismic coherency and spectral decomposition cube[J].Abstract of Istanbul ’97 International Geophysical Conference and Exposition,1997,7-10

[17] Partyka G,Gridley J,Lopez J.Interpretational applications of spectral decomposition in reservoir characterization[J].The Leading Edge,1999(3):353-360[18] Marfurt K J,Kirlin R L.Narrow-band spectral analysis and thin-bed tuning[J].Geophysics,2001,66(4):1274-1283

[19] 吴小平,徐果明.利用共轭梯度法的电阻率三维反演研究[J].地球物理学报,2000,43(3):420-426 Wu X P,Xu G M.Study on 3-D resistivity inversion using conjugate gradient method[J].Chinese Journal of Geophysics,2000,43(3):420-426

[20] 李冰,刘洪,李幼铭.三维地震数据离散光滑插值的共轭梯度法[J].地球物理学报,2002,45(5):691-699 Li B,Liu H,Li Y M.3-D seismic data discrete smooth interpolation using conjugate gradient[J].Chinese Journal of Geophysics,2002,45(5):691-699

[21] 霍志周,熊登,张剑峰.预条件共轭梯度法在地震数据重建方法中的应用[J].地球物理学报,2013,56(4):1321-1330 Huo Z Z,Xiong D,Zhang J F.Application of the preconditioned conjugate gradient method to reconstruction of seismic data[J].Chinese Journal of Geophysics,2013,56(4):1321-1330

[22] 吴小羊,刘天佑.基于时频重排的地震信号Wigner-Ville分布时频分析[J].石油地球物理勘探,2009,44(2):201-205 Wu X Y,Liu T Y.Time-frequency analysis on Wigner-Ville distribution of seismic signal based on time-frequency rearrangement[J].Oil Geophysical Prospecting,2009,44(2):201-205

[23] 戴永寿,王俊岭,王伟伟,等.基于高阶累积量ARMA模型线性非线性结合的地震子波提取方法研究[J].地球物理学报,2008,51(6):1851-1859 Dai Y S,Wang J L,Wang W W,et al.Seismic wavelet extraction via cumulant based ARMA model approach with linear and nonlinear combination[J].Chinese Journal of Geophysics,2008,51(6):1851-1859

[24] 钱勇先,王慎中.重采样法高分辨率处理[J].石油地球物理勘探,1994,29(4):498-503 Qian Y X,Wang S Z.High-resolution data Processing by resembling method[J]. Oil Geophysical Prospecting,1994,29(4):498-503

[25] Chopra S,Castagna J,Xu Y.Thin-bed reflectivity inversion and some applications[J].The First Break,2009,27(5):17-24

(编辑：顾石庆)

Application of spectral inversion for tight thin-bed sand body prediction

Chi Huanzhao1,Liu Cai2,Shan Xuanlong1,Lu Qi2

(1.CollegeofEarthSciences,JilinUniversity,Changchun130061,China;2CollegeofGeo-ExplorationScienceandTechnology,JilinUniversity,Changchun130026,China)

Due to the need of selecting target areas and deploying horizontal wells,it is very important to understand the lateral distribution of tight sand bodies in unconventional tight sandstone reservoir exploration and development.The tight sandstone is characterized by the development of Cretaceous thin-beds for Fuyu reservoir in the north of Songliao Basin.And the identification accuracy of the thin-bed sand bodies is restricted by the shielding effect from the strong reflection interface T2.Therefore,we study spectral inversion technique for the prediction of the tight thin-bed sand bodies in Fuyu formation,Daqing Placanticline.In order to ensure the accuracy of spectral inversion,fast matching pursuit algorithm in complex domain is used to improve the precision of the time-frequency conversion and the accuracy of time-frequency spectrum of reflection coefficient is increased by the time-varying wavelet extraction.The results of spectral inversion show that:the seismic reflection coefficient volumes from spectral inversion are convoluted with the broadband zero-phase wavelet,thus obtained broadband preserved amplitude seismic data effectively suppress the shielding effect from T2in target area.The results also indicate the improvements of composite reflections overlay,the top and the bottom interface of the thin-bed sandstone detected clearly,and the favorable correspondence between amplitude attributes and sand bodies.

spectral inversion,fast matching pursuit,time-varying wavelet,tight sandstone,shielding effect from T2

2014-07-22;改回日期：2014-10-30。

迟唤昭(1985—),男,博士在读,主要从事油气地质与地震方法方面研究。

鹿琪(1972—),女,博士,教授,主要从事探地雷达和地震数据处理研究。

国家重点基础研究发展计划(973计划)(2013CB429805,2009CB219301)、国家自然科学基金(41174080,41340039)和公益性行业科研专项(201011078)共同资助。

P631

A

1000-1441(2015)03-0337-08

10.3969/j.issn.1000-1441.2015.03.013