地球物理学进展  2016, Vol. 31 Issue (5): 2326-2332   PDF    
瞬变电磁虚拟波场成像方法及其对未爆弹探测的试验研究
范涛1, 程建远1, 王保利1, 姚伟华1, 王继矿1, 李貅2, 李博凡3     
1. 中国煤炭科工集团西安研究院有限公司, 西安 710077
2. 长安大学地质工程与测绘学院, 西安 710054
3. 西安科技大学地质与环境学院, 西安 710054
摘要: 未爆弹严重威胁人类的安全,瞬变电磁是一种对其有效的探测方法.探讨了瞬变电磁波场变换理论,实现了对瞬变电磁信号中波场特性的提取,提出了基于虚拟波场波形特征的速度分析反演成像方法,计算了三维数值模型和水槽物理模拟数据的反演结果,取得了良好的成像效果.结合未爆弹探测的工程实例,对瞬变电磁虚拟波场反演成像方法的实用性和有效性进行了验证.研究表明:瞬变电磁虚拟波场反演成像方法能够解决微小低阻异常体精确定位的难题.
关键词瞬变电磁法     波场变换     速度分析     未爆炸弹    
Imaging method of TEM pseudo wave-field and application of unexploded ordnance detection
FAN Tao1 , CHENG Jian-yuan1 , WANG Bao-li1 , YAO Wei-hua1 , WANG Ji-kuang1 , LI Xiu2 , LI Bo-fan3     
1. Xi'an Research Institute, CCTEG, Xi'an 710077, China
2. College of Geology Engineering and Geomatics, Chang'an University, Xi'an 710054, China
3. College of Geology & Environment, Xi'an University of Science and Technology, Xi'an 710054, China
Abstract: Unexploded ordnance today grows into a critical issue of environment, transient electromagnetic is an effective method for its detection. Authors first discussed the theory of transient electromagnetic wave-field transforming, realized the extraction of wave-field characteristics from transient electromagnetic signal, proposed inversion imaging method of velocity analysis based on waveform characteristics of pseudo wave-field, calculated the inversion results of a three-dimensional numerical model data and a physical simulation data, achieved excellent imaging effect. At last, the practicability and effective of the inversion imaging method of TEM pseudo wave-field were proved by the actual exploration of the unexploded ordnance detection. The conclusion show that the inversion imaging method of TEM pseudo wave-field can solve problem of accurate positioning of tiny low resistivity abnormal body.
Key words: transient electromagnetic method     wave-field transformation     velocity analysis     unexploded ordnance    
0 引言

未爆弹主要包括战争及其他军事活动期间,布设和投放的地雷、炸弹等遗留武器,它不仅威胁着人类生命安全,限制了人类正常的活动空间,还降低了土地资源的利用率,加大了土地开发成本,干扰地方正常的社会生产和生活(李小康等,2010).

探测未爆弹的主要地球物理方法包括磁法、电磁法、地质雷达、微重力法和放射性法.其中磁法应用较多,但准确率仅为10%左右(Gasperikova and Beard, 2007),且对非磁性金属物体不适用,放射性法只适用于贫铀弹等带有放射性的未爆弹,地质雷达和微重力法对于应用场地的环境要求相对较高,电磁法相对限制较小,探测精度较高,目前国际上从事未爆弹探测工作的地球物理工作者,都倾向于使用电磁方法(李小康等,2010).

瞬变电磁法属于时间域电磁法,多匝小线圈装置施工灵活方便,可用于埋深较浅(小于10 m)的未爆弹探测.但多匝发射、接收小线圈也会带来较大的电感影响(姜志海等,2007杨海燕等,2007Sijoy and Chaturvedi, 2008李文尧和武中华,2010),造成实测数据畸变严重,使用常规的处理技术解释出的异常区域偏大,不能形成异常圈闭区,难以起到对未爆弹异常精确定位的效果(范涛等,2014),急需研究新的精细反演解释方法.

瞬变电磁拟地震成像技术是一个前沿研究方向(李貅等,2007薛国强等,2008).瞬变电磁满足的微分方程是一个扩散方程,对成像是不利的;地震波场满足的波动方程,则有利于成像.因此,国内外学者通过数学变换,将瞬变电磁场转化为波场(Lee,1987Lee et al., 1989Zhdanov and Traynin, 1995Zhdanov and Portniaguine, 1997),进而应用成熟的地震反演成像方法,解决复杂地质条件下瞬变电磁法的反演问题.

Lee等建立了满足时域扩散方程的电磁场与虚拟波场的关系式(Lee et al., 1989),这一变换不仅可以将对应地电模型的波动方程模拟结果变换成时域电磁响应,更重要的是它的逆变换意味着可以将高导电媒质中电磁场所满足的扩散方程转换为无耗媒质中电磁场所满足的无阻尼波动方程.

近年来,利用这一关系式开展的瞬变电磁波场变换研究取得了较大进展.国外Tae J应用电磁波的走时实现了变换后的波场层析成像(Tae J et al., 2002).国内陈本池、李金铭等利用奇异值分解实现了瞬变电磁波场变换,并利用有限差分开展了二维拟地震解释(陈本池等,1999);李貅等提出了分段波场变换的最优化算法、虚拟波场的三维曲面延拓成像方法,并利用合成孔径算法提高了成像精度(李貅等, 2005, 2010, 2012);薛国强等运用波场变换算法转换出虚拟波场子波,研究了子波宽度压缩技术(Xue et al., 2007薛国强等,2011);朱宏伟等研究了从虚拟波场数据中提取有效成像信息的方法(朱宏伟等,2010);范涛研究了瞬变电磁虚拟波场的速度提取问题,使实测数据拟地震成像成为可能(范涛,2011);张军等研究了波场变换的高分辨算法(张军等,2011);戚志鹏等改进了李貅的波场变换算法,实现了预条件正则化的全时域反变换(戚志鹏等,2013);程久龙等在煤矿井下瞬变电磁探测中初步应用了波场变换方法,使用波形来反映地电界面(程久龙等,2013).目前,国内外学者对瞬变电磁波场变换方面的主要研究工作集中在数值模拟方面,实际应用研究还很少.

笔者基于满足扩散方程的瞬变电磁场与满足波动方程的波场之间的转换关系,提取了瞬变电磁信号中的虚拟波场信息,利用其波形特性进行波场速度分析反演,增强了瞬变电磁法对未爆弹的识别能力.

1 方法原理 1.1 波场变换

扩散方程到波动方程转换的表达式为

(1)

将(1)式离散为

(2)

(2)式写为矩阵形式为

(3)

其中 A=[aijwj]m×nU=[uj]n×1为虚拟子波,F=[hi]m×1为接收的瞬变场时域信号.

将(3)式转化为

(4)

只要 A是列满秩矩阵,A T A就是对称正定矩阵,因此可以利用共轭梯度法.当然,A T A的条件数也会大于 A,方程病态性会加剧,因此令

(5)

针对矩阵 A (v)可构造预条件子M (v),形成新方程为

(6)

矩阵 M (v)-1 A (v)接近单位阵.用共轭梯度法进行内部迭代,则外部正则化迭代会很快收敛,实现全时域的波场变换.

1.2 速度分析

经过前一步工作,瞬变电磁数据已转化为类似于地震自激自收的波形数据.

此时可将地质模型用m(x)表示,m是多种参数组合的参数集,考虑到需要求取的地质模型参数只有速度,为简化计算,这里可假设介质参数仅包含纵波速度vp(x),方程可用如下声波波动方程来进行描述,公式为

(7)

式中,t为时间,d obs(xg, t; xs)为接收到的地震数据,s (xs, t)为震源的函数,xg为接收点空间位置,xs为震源空间位置.

欲求解该反问题,可采用非线性共轭梯度法.首先梯度表达式可写为

式中, u (x, t)代表波场;λ是将虚拟波形记录与模拟合成记录的残差记录在接收点,从最大时间向零时间反传得到的波场记录(Tarantola et al., 1984).

之后可利用本次计算的负梯度方向-▽Em与前一次的共轭梯度方向进行线性组合作为下一次迭代的搜索方向,即

-▽Em0为初始速度模型计算得到的负梯度,-▽Emk为第k次迭代时模型mk计算的负梯度,γ(k-1)γ(k)分别为k-1次和第k次的共轭梯度,βk为线性加权系数,其计算公式为

最终,迭代格式可写为

(8)

式中,ε为一较小的数.

2 正演模拟 2.1 数值模拟

为验证算法对复杂三维模型的反演效果,设计了如图 1所示的模型,该模型为3层介质中夹有1个未爆弹的三维模型,模型参数如图 1中所示.其中未爆弹设置为不规则梭型体,以尽量真实的模拟弹体,长度1.5 m,宽度0.25 m,左侧弹体厚度小,右侧弹体厚度大,底界面埋深3 m.模型正演采用回线源激发瞬变场的方式,使用孙怀凤开发的三维时域有限差分正演程序(孙怀凤,2013孙怀凤等,2013)取得正演数据.回线发射框为2 m×2 m,发射电流为1A,发射框和测线布置如图 1b,每条测线上有19个测点,测点间距0.25 m.

图 1 模型示意图 (a)模型立体图;(b)模型俯视图. Figure 1 Schematic diagram of the model

图 1中的测点数据进行反演,可得如图 2所示的反演结果.对比(a) (b)两图可以看出,常规处理结果尽管能反映出未爆弹的低阻特征,但异常形态表现为一向深部凹陷的巨大等值线扭曲,无法反映模型设置的异常厚度,且形态为左右对称,与模型有较大差异,同一深度上不易提取出弹体的细节信息;虚拟波场反演成像结果中,未爆弹的厚度、位置均与模型设置吻合较好,且左右的弹体长度、厚度变化差异也有较好体现,同一深度上对弹体的细节信息体现更优.同时可以看到,由于未爆弹低阻异常体的存在,中间的灰岩层位在处理结果上并未表现出高阻特征,但埋深和厚度与模型设置基本一致.

图 2 模型反演结果 (a)常规处理结果图; (b)虚拟波场反演成像图. Figure 2 Inversion results of the model
2.2 物理模拟

为进一步验证算法对实际数据的反演效果,进行了物理模拟试验.以铜棒模拟未爆弹,物理模型设置如图 3所示,模型参数如表 1所示.采用多匝小回线源激发、中心磁探头接收的施工方式取得数据,回线发射框为0.4 m×0.4 m×10匝,发射电流为1.5 A,测线布置如图 3,测线上有13个测点.

图 3 模型示意图 Figure 3 Schematic diagram of the model

表 1 模型参数表 Table 1 Model parameter table

图 3中水槽上方的测点数据进行反演,可得如图 4所示的反演结果.对比(a) (b)两图可以看出,常规处理结果尽管能反映出铜棒的低阻特征,但异常范围偏大,无法反映出铜棒厚度;虚拟波场反演成像结果中,铜棒的低阻异常与实际吻合更好,厚度与模型基本一致,对空气与水面的界面也有很好的反映.

图 4 模型反演结果 (a)常规处理结果图; (b)虚拟波场反演成像图. Figure 4 Inversion results of the model
3 探测实例

为某兵工厂探测两个不同大小的未爆弹(兵工厂事先埋设),弹体目标如图 5所示,30弹直径30 mm,长度约15 cm,57弹直径57 mm,长度约25 cm.

图 5 目标弹体(30弹和57弹) Figure 5 Target ordnance

探测工作结束后,经实际挖掘确认,30弹的弹体埋藏位置与测线布置如图 6所示,瞬变电磁探测采用多匝小回线源激发、中心磁探头接收的施工方式取得数据,回线发射框为2 m×2 m×10匝,发射电流为3 A,布置了3条测线,每条测线15个测点.

图 6 30弹探测工作示意图 Figure 6 Diagram of detection work

常规处理结果与虚拟波场反演成像结果如图 7所示.图 7a中有多个异常,难以判断炮弹位置,且低阻异常最强的位置为X=0.4 m处,与实际情况不相符;图 7b中则可以看到,仅X=2.6 m处有一明显的低阻异常,与实际情况吻合.

图 7 30弹L2线对比结果 (a)常规处理结果图; (b)虚拟波场反演成像图. Figure 7 Results of L2

虚拟波场反演成像三维空间结果如图 8所示,对比图 6可以看出,主异常三维空间分布与未爆弹实际空间位置吻合度很高.

图 8 30弹反演结果三维成像图 Figure 8 3D imaging results of inversion

57弹的弹体埋藏位置与测线布置如图 9所示,瞬变电磁探测采用多匝小回线源激发、中心磁探头接收的施工方式取得数据,回线发射框为2 m×2 m×10匝,发射电流为3A,布置了3条测线,每条测线16个测点.

图 9 57弹探测工作示意图 Figure 9 Diagram of detection work

常规处理结果与虚拟波场反演成像结果如图 10所示.图 10a中浅部有低阻盖层,深部以高阻为主,在H=1.7 m到H=2.5 m深度出现低阻异常分布,其中X=1.2 m处异常与实际情况较相符,但除此外同一深度层内还有多个异常分布,为准确确定未爆弹位置带来干扰;图 10b中则可以看到,仅X=1.2 m处有一明显的低阻异常,与实际情况吻合.

图 10 57弹L2线弹对比结果 (a)常规处理结果图; (b)虚拟波场反演成像图. Figure 10 Results of L2

三维空间结果如图 11所示,对比图 9可以看出,主异常三维空间分布与未爆弹实际空间位置吻合度很高.

图 11 57弹反演结果三维成像图 Figure 11 3D imaging results of inversion
4 结论 4.1

瞬变电磁场满足的扩散方程与地震波场满足的波动方程之间存在数学转换关系式,因此可以提取瞬变电磁法资料中的波场特性进行拟地震反演解释.转换过程中压制了电磁波传播过程中与频散和衰减有关的特性,增强了与传播有关的特性,能提高瞬变电磁方法对地下电性界面的分辨能力.

4.2

数值模拟和物理模拟数据的处理效果说明这种基于虚拟波场的反演成像方法,相比常规瞬变电磁处理技术,克服了低阻屏蔽作用,减小了电法体积效应对反演结果的影响,对异常体的底界面反映清晰,能够精细定位微小低阻异常体的准确位置.探测实例通过反演结果与常规处理结果、实际挖掘揭露结果的对比,验证了该方法对微小异常体精确定位的效果,可以应用于未爆弹探测实际工作中.

致谢 感谢审稿专家提出的修改意见和编辑部的大力支持!
参考文献
[] Chen B C, Li J M, Zhou F T .1999b. Quasi wave equation migration of transient electromagnetic field[J]. Oil Geophysical Prospecting , 34 (5) : 546–554.
[] Chen B C, Zhou F T, Li J M .1999a. The wavefield transformation study of transient electromagnetic field[J]. Geophysical & Geochemical Exploration , 23 (3) : 195–201.
[] Cheng J L, Qiu H, Ye Y T, et al .2013. Research on wave-field transformation and data processing of the mine transient electromagnetic method[J]. Journal of China Coal Society , 38 (9) : 1646–1650.
[] Fan T .2011. A 3D continuous velocity analysis of TEM fictitious wave-field and its application to tunnel advanced prediction[J]. Geophysical & Geochemical Exploration , 35 (2) : 243–247.
[] Fan T, Zhao Z, Wu H, et al .2014. Research on inductance effect removing and curve offset for mine TEM with multi small loops[J]. Journal of China Coal Society , 39 (5) : 932–940.
[] Gasperikova E, Beard L P .2007. Special issue on geophysics applied to detection and discrimination of unexploded ordnance[J]. Journal of Applied Geophysics, 61 (3-4) : 165–167. DOI:10.1016/j.jappgeo.2006.06.007
[] Jiang Z H, Yue J H, Liu S C .2007. Mine transient electromagnetic observation system of small multi-turn coincident configuration[J]. Journal of China Coal Society , 32 (11) : 1152–1156.
[] Lee K H, Liu G, Morrison H F .1989. A new approach to modeling the electromagnetic response of conductive media[J]. Geophysics, 54 (9) : 1180–1192. DOI:10.1190/1.1442753
[] Lee S, McMechan G A, Aiken C L V .1987. Phase-field imaging:The electromagnetic equivalent of seismic migration[J]. Geophysics, 52 (5) : 678–693. DOI:10.1190/1.1442335
[] Lee T J, Suh J H, Kim H J, et al .2002. Electromagnetic traveltime tomography using an approximate wavefield transform[J]. Geophysics, 67 (1) : 68–76. DOI:10.1190/1.1451344
[] Li W Y, Wu Z H .2010. Accurate formula of self-inductance coefficient of rectangular coil in transient electromagnetic methods[J]. Geology and Exploration , 46 (1) : 160–164.
[] Li X, Qi Z P, Xue G Q, et al .2010. Three dimensional curved surface continuation image based on TEM pseudo wave-field[J]. Chinese Journal of Geophysics , 53 (12) : 3005–3011. DOI:10.3969/j.issn.0001-5733.2010.12.025
[] Li X, Xue G Q, Guo W B .2007. Research progress in TEM pseudo-seismic imaging[J]. Progress in Geophysics , 22 (3) : 811–816. DOI:10.3969/j.issn.1004-2903.2007.03.023
[] Li X, Xue G Q, Liu Y A, et al .2012. A research on TEM imaging method based on synthetic-aperture technology[J]. Chinese Journal of Geophysics , 55 (1) : 333–340. DOI:10.6038/j.issn.0001-5733.2012.01.034
[] Li X, Xue G Q, Song J P, et al .2005. An optimize method for transient electromagnetic field-wave field conversion[J]. Chinese Journal of Geophysics , 48 (5) : 1185–1190. DOI:10.3321/j.issn:0001-5733.2005.05.029
[] Li X K, Yang L, Liu L .2010. Overview on unexploded ordnance problem and a solution:The geophysical scheme[J]. China Mining Magazine , 19 (S1) : 187–191.
[] Qi Z P, Li X, Wu Q, et al .2013. A new algorithm for full-time-domain wave-field transformation based on transient electromagnetic method[J]. Chinese Journal of Geophysics , 56 (10) : 3581–3595. DOI:10.6038/cjg20131033
[] Sijoy C D, Chaturvedi S .2008. Calculation of accurate resistance and inductance for complex magnetic coils using the finite-difference time-domain technique for electromagnetics[J]. IEEE Transactions on Plasma Science, 36 (1) : 70–79. DOI:10.1109/TPS.2007.914693
[] Sun H F. 2013. Three-dimensional transient electromagnetic responses of water bearing structures in tunnels and prediction of water inrush sources (in Chinese)[Ph. D. thesis]. Jinan:Shandong University.
[] Sun H F, Li X, Li S C, et al .2013. Three-dimensional FDTD modeling of TEM excited by a loop source considering ramp time[J]. Chinese Journal of Geophysics , 56 (3) : 1049–1064. DOI:10.6038/cjg20130333
[] Tarantola A .1984. Inversion of seismic reflection data in the acoustic approximation[J]. Geophysics, 49 (8) : 1259–1266. DOI:10.1190/1.1441754
[] Xue G Q, Li X, Di Q Y .2008. Research progress in TEM forward modeling and inversion calculation[J]. Progress in Geophysics , 23 (4) : 1165–1172.
[] Xue G Q, Li X, Qi Z P, et al .2011. Study of sharpening the TEM pseudo-seismic wave-form[J]. Chinese Journal of Geophysics , 54 (5) : 1384–1390. DOI:10.3969/j.issn.0001-5733.2011.05.028
[] Xue G Q, Yan Y J, Li X .2007. Pseudo-Seismic wavelet transformation of transient electromagnetic response in engineering geophy exploration[J]. Journal of Geophysical Research Letters, 34 (16) . DOI:10.1029/2007GL031116
[] Yang H Y, Yue J H, Hu W W, et al .2007. The characteristics of the early signal in tem affected by self-induction of Multiturn coil[J]. Computing Techniques for Geophysical and Geochemical Exploration , 29 (2) : 96–98.
[] Zhang J, Li X, Zhao Y, et al .2011. A technology research of high-resolation imaging for the transient electromagnetic pseudo wave field[J]. Progress in Geophysics , 26 (3) : 1077–1084. DOI:10.3969/j.issn.1004-2903.2011.03.038
[] Zhdanov M S, Portniaguine O .1997. Time-domain electromagnetic migration in the solution of inverse problems[J]. Geophysics Journal International, 131 (2) : 293–309. DOI:10.1111/gji.1997.131.issue-2
[] Zhdanov M S, Traynin P N, Portniaguine O .1995. Resistivity imaging by time domain electromagnetic migration (TDEMM)[J]. Exploration Geophysics, 26 (3) : 186–194. DOI:10.1071/EG995186
[] Zhu H W, Li X, Zhang J, et al .2010. Information collecting technology in 3-D pseudo-seismic imaging of transient electromagnetics[J]. Progress in Geophysics , 25 (5) : 1648–1656. DOI:10.3969/j.issn.1004-2903.2010.05.017
[] 陈本池, 李金铭, 周凤桐.1999b. 瞬变电磁场拟波动方程偏移成像[J]. 石油地球物理勘探, 34 (5) : 546–554.
[] 陈本池, 周凤桐, 李金铭.1999a. 瞬变电磁场的波场变换研究[J]. 物探与化探, 23 (3) : 195–201.
[] 程久龙, 邱浩, 叶云涛, 等.2013. 矿井瞬变电磁法波场变换与数据处理方法研究[J]. 煤炭学报, 38 (9) : 1646–1650.
[] 范涛.2011. TEM虚拟波场三维连续速度分析及其在隧道超前预报中的应用[J]. 物探与化探, 35 (2) : 243–247.
[] 范涛, 赵兆, 吴海, 等.2014. 矿井瞬变电磁多匝回线电感影响消除及曲线偏移研究[J]. 煤炭学报, 39 (5) : 932–940.
[] 姜志海, 岳建华, 刘树才.2007. 多匝重叠小回线装置的矿井瞬变电磁观测系统[J]. 煤炭学报, 32 (11) : 1152–1156.
[] 李文尧, 武中华.2010. 瞬变电磁法矩形线圈自感的精确表达式[J]. 地质与勘探, 46 (1) : 160–164.
[] 李小康, 杨磊, 刘磊.2010. 未爆弹药问题及其地球物理解决方案综述[J]. 中国矿业, 19 (S1) : 187–191.
[] 李貅, 戚志鹏, 薛国强, 等.2010. 瞬变电磁虚拟波场的三维曲面延拓成像[J]. 地球物理学报, 53 (12) : 3005–3011. DOI:10.3969/j.issn.0001-5733.2010.12.025
[] 李貅, 薛国强, 郭文波.2007. 瞬变电磁法拟地震成像研究进展[J]. 地球物理学进展, 22 (3) : 811–816. DOI:10.3969/j.issn.1004-2903.2007.03.023
[] 李貅, 薛国强, 刘银爱, 等.2012. 瞬变电磁合成孔径成像方法研究[J]. 地球物理学报, 55 (1) : 333–340. DOI:10.6038/j.issn.0001-5733.2012.01.034
[] 李貅, 薛国强, 宋建平, 等.2005. 从瞬变电磁场到波场的优化算法[J]. 地球物理学报, 48 (5) : 1185–1190. DOI:10.3321/j.issn:0001-5733.2005.05.029
[] 戚志鹏, 李貅, 吴琼, 等.2013. 从瞬变电磁扩散场到拟地震波场的全时域反变换算法[J]. 地球物理学报, 56 (10) : 3581–3595. DOI:10.6038/cjg20131033
[] 孙怀凤. 2013.隧道含水构造三维瞬变电磁场响应特征及突水灾害源预报研究[博士论文].济南:山东大学. http://www.cnki.com.cn/Article/CJFDTotal-YSLX201503026.htm
[] 孙怀凤, 李貅, 李术才, 等.2013. 考虑关断时间的回线源激发TEM三维时域有限差分正演[J]. 地球物理学报, 56 (3) : 1049–1064. DOI:10.6038/cjg20130333
[] 薛国强, 李貅, 底青云.2008. 瞬变电磁法正反演问题研究进展[J]. 地球物理学进展, 23 (4) : 1165–1172.
[] 薛国强, 李貅, 戚志鹏, 等.2011. 瞬变电磁拟地震子波宽度压缩研究[J]. 地球物理学报, 54 (5) : 1384–1390. DOI:10.3969/j.issn.0001-5733.2011.05.028
[] 杨海燕, 岳建华, 胡文武, 等.2007. 多匝回线的自感对瞬变电磁早期信号的影响特征[J]. 物探化探计算技术, 29 (2) : 96–98.
[] 张军, 李貅, 赵莹, 等.2011. 瞬变电磁虚拟波场高分辨成像技术研究[J]. 地球物理学进展, 26 (3) : 1077–1084. DOI:10.3969/j.issn.1004-2903.2011.03.038
[] 朱宏伟, 李貅, 张军, 等.2010. 瞬变电磁法三维拟地震成像信息提取技术[J]. 地球物理学进展, 25 (5) : 1648–1656. DOI:10.3969/j.issn.1004-2903.2010.05.017