張 浩,張劍鋒
1中國科學院地質與地球物理研究所中國科學院油氣資源研究重點實驗室,北京 100029 2中國科學院研究生院,北京 100049
起伏地表采集數據的三維直接疊前時間偏移方法
張 浩1,2,張劍鋒1
1中國科學院地質與地球物理研究所中國科學院油氣資源研究重點實驗室,北京 100029 2中國科學院研究生院,北京 100049
提出一種可對起伏地表采集的三維地震資料直接進行偏移成像的疊前時間偏移方法和流程.它用兩個等效速度描述近地表和上覆層對地震波傳播的影響,可對炮、檢點不在同一水平面的三維地震資料直接進行疊前時間偏移處理.該方法不對近地表地震波傳播做垂直出、入射假定,因此可適應高速層出露等不存在明顯低、降速帶情況.描述近地表和上覆層的兩個等效速度參數可依據偏移道集的同相軸是否平直來確定,避免了確定近地表速度的困難;而對已知近地表速度的情況,則可進一步修正近地表速度,獲得更好的成像效果.用三維起伏地表的理論數據和中國東部某工區實際數據驗證了所發展方法和處理流程的有效性和實用性.
起伏地表,三維成像,近地表速度,疊前時間偏移
隨著我國陸上油氣勘探向西部和南方地區的深入,復雜地表已成為制約地震資料處理的關鍵因素.由于地表起伏較大、高速的老地層出露,現行的以靜校正處理為基礎的陸上地震資料處理方法和流程遇到困難.浮動基準面技術[1-2]可解決地表起伏較大問題,但高速層出露導致垂直入射和出射的靜校正應用條件不成立,繼續使用靜校正方法將帶來較大的誤差[3].基于波場延拓的基準面重建方法[4]可同時解決地表起伏和高速層出露問題,但這一方法是以獲得準確的近地表速度模型為基礎的;由于起伏地表采集的地震資料信噪比較低,獲得準確的近地表速度模型是一個相當困難的任務[5].此外,基于波場延拓的基準面重建方法導致重建數據的地震道數激增,極大地增加了后續偏移處理的計算量.
就實際中遇到的大量地質目標而言,盡管地表復雜、斷層等構造斷裂發育,但除近地表以外的上覆地層的速度橫向變化并不是很劇烈,這就給疊前時間偏移的應用提供了可能性.疊前時間偏移的主要優點是只用一個等效速度(即均方根速度)描述地震波傳播,可對一類傾角、斷層較為復雜,但速度橫向變化不是很劇烈的構造較好的成像[6].獨立的、單一等效速度使得疊前時間偏移的速度建模變得簡單.本文進一步發展疊前時間偏移方法,通過引入第二個等效參數描述近地表地震波傳播,發展可對炮、檢點不在同一水平面的三維地震資料直接進行疊前時間偏移處理的成像方法,以此避免高速層出露和近地表速度建模這兩個困難.
現行的直接對起伏地表采集數據進行偏移、避免使用靜校正技術的直接偏移方法主要是基于疊前深度偏移算法[7-8].這一方法需要準確的近地表速度和地下構造的層速度模型,而這兩類速度又相互影響,與海上或水平地表情況下確定深度偏移的層速度模型相比,這類方法的速度建模更加困難.文獻[9]盡管在計算地震波走時時引入了平均速度和均方根速度,但這兩個速度是由層速度計算得到的,它們是成像點和炮點(或檢波點)的函數;由于在同一成像點有多個均方根速度存在(隨炮點或檢波點變化),這一速度是不能如疊前時間偏移那樣通過簡單的掃描方式得到的,只能通過先求取層速度模型,進而解析計算得到.因此,文獻[9]仍屬于疊前深度偏移范疇,面臨著與疊前深度偏移相同的困難.董春暉等[10]發展了一個基于雙速度參數的二維直接疊前時間偏移方法,有效地將疊前時間偏移技術推廣應用于起伏地表問題.此外,實際起伏地表問題總是三維的,炮、檢點不在同一方位角平面,導致走時、等效偏移距(浮動基準面上的偏移距)的計算較二維情況的更復雜.本文則在確定兩個等效速度參數、處理三維起伏地表這兩方面提升了文獻[10]的方法,并且將這一方法應用于三維起伏地表采集的實際地震資料.
本文的關鍵有兩點,一是從疊前深度偏移的單程波方程出發,解析給出了三維情況下由描述近地表和上覆層速度的兩個等效速度表達的地震波走時、幅值和等效偏移距,并發展了基于表插值的快速算法;二是發展了基于偏移道集同相軸的雙曲彎曲和局部微彎曲,分別確定兩個等效速度參數的方法和算法流程.文中也應用了隨方位角和構造傾角變化的三維時變偏移孔徑,這使得本文疊前時間偏移方法即有效地壓制了偏移噪音,也減少了偏移計算量.三維起伏地表的理論數據和實際地震資料證明了本文方法和處理流程的正確和有效性.
假設三維非均勻介質可近似為層狀介質,不失一般性,令炮點或檢波點在坐標原點;在波數-頻率域,基于三維深度偏移的相移法,該炮點或檢波點的波場深度延拓可用時間深度表示為

若將第m層介質界面作為基準面,可近似式(1)中右端指數項中的相移量為

將式(2)代入(1)式并做空間傅里葉反變換,可得到空間-頻率域波場:

式(3)是一個二重震蕩積分,可利用穩相點原理[11]求得漸進解為

式(4)中,定義

其中px和py是射線參數的零點,可由下式解得:

引入新的變量pr,定義p0x=prcosφ,p0y=prsinφ,其中φ是坐標點(x,y)的方位角,則(7)和(8)式可轉化為單變量的方程:



而代入式(6),則可進一步得到地震波的幅值.
令a2=0,即假設基準面以上是明顯的低速帶,則式(11)簡化為

這就是應用靜校正技術的走時公式.這表明式(11)涵蓋了現行的靜校正方法.若假定介質為上、下層速度分別是V0和Vrms的雙層介質,式(9)和(11)則可由Snell定律直接導出,這從另一個角度證明了式(9)和(11)的正確性;式(9)和(11)既擴展了Snell定律的結果,也給出了非均勻介質情況下等效速度參數的物理含義.
盡管可解析求解[12]式(9),但對三維偏移過程中大量的成像點、炮點、檢波點組合而言,直接計算涉及了太多的計算量.既然3個無量綱系數完全決定了式(9)的四次方程,我們發展了基于表插值的快速算法.3個無量綱系數均有明確的物理意義,就疊前時間偏移方法而言,其對大于70°角度傳播的高角度波場是不感興趣的,因此可令0≤a3≤tan70;即使高速層出露,近地表等效速度也不可能大于深層速度,因此可定義0≤a2≤1.0;疊前偏移成像更多的是關心中、深層,若假定從T2≥T1開始成像計算,則可定義0≤a1≤1.0.這樣,若對這3個無量綱系數等間距離散采樣,預先求解對應的η并建立一個三維表,就可在偏移過程中通過計算3個無量綱系數,直接在表中拾取對應的η.炮點至成像點和檢波點至成像點可用同一個表.這一基于表的算法可極大提高三維偏移的計算效率.
盡管(9)、(11)以及(6)式的走時和振幅值是基于(1)式的層狀介質假設推出的,但通過允許兩個等效速度V0和Vrms以及基準面橫向變化,式(9)、(11)以及(6)式的解可代表三維橫向弱非均勻介質中地震波的走時和振幅值;對較強的速度橫向變化情況,只能采用基于層速度模型的深度域延拓方法;文獻[13]通過引入速度梯度等信息,對較強橫向變速情況下的走時計算進行了討論.
等效速度參數Vrms和基準面深度的橫向變化是由成像點的橫向坐標(x,y)決定的.這樣,Vrms將與基準面上的等效偏移距構成常規疊前時間偏移的雙曲走時關系;如果令a2=0,則本文算法可退化為浮動基準面疊前時間偏移方法.另一個等效參數V0即可隨炮、檢點橫向變化,也可隨成像點橫向變化.
第2節的研究表明,可用兩個等效速度參數描述地震波在近地表和上覆層中的傳播;而如何準確獲得這兩個等效參數,是本文方法實際應用的關鍵問題.本文的兩個等效速度參數在每個成像點是惟一的,這兩個參數完全決定了該成像點處共反射點道集(CRP)的表現,因此可在偏移過程中根據該CRP道集的表現來決定它們;不同于文獻[9]的雙速度參數,它們是成像點和炮點(或檢波點)的函數,若干個不同的參數決定了該成像點處CRP道集的表現,因此不能由單個CRP道集來決定該速度參數.實際上,文獻[9]中的雙速度參數只是用于近似計算地震波的走時,它們是由已知的層速度模型解析計算得到的.
從偏移成像的角度出發,CRP道集中同相軸是否平直是確定兩個等效速度的標準,這也是本文確定等效速度的出發點.但兩個等效速度參數共同影響了同相軸的平直,只有區分這兩者對同相軸彎曲的獨特影響,才能有效地確定這兩個參數.

圖1 等效速度參數對走時影響分析示意圖Fig.1 Illustration of equivalent velocity parameter analysis using 2Dmodel including topography

圖2 等效速度參數對走時殘差影響Fig.2 Influence of equivalent velocity parameters on residual traveltime
本文用一個簡單的二維模型來考察兩個速度參數變化對走時的影響.如圖1所示,起伏地表按正弦函數變化,其波長為1000m,波鋒與波谷間的高程為200m;成像點在模型中心,對應該成像點的基準面如圖1中虛線所示,較地表高程曲線的波谷低0.1s(時間深度),而成像點在基準面下1.5s處;基準面上、下的介質是均勻的,速度分別為1800m/s和2500m/s.圖2a給出了由(11)式求得的地面上炮點至成像點再反射到地表檢波點的地震波走時隨偏移距(水平距離)的變化曲線.若令V0=1800m/s,Vrms=2000,2500,3000m/s,利用式(11)做動校正,可得剩余走時曲線(如圖2b),而圖2c則給出了隨基準面上的等效偏移距變化的剩余走時曲線;圖2b中的點線是最佳擬合剩余走時的雙曲線,而圖2c的點線與剩余走時曲線的實線重合.觀察圖2b和圖2c可發現,速度參數Vrms影響了同相軸的整體彎曲程度,而等效偏移距使得這種彎曲更好地展現了如常規時間偏移的雙曲特征.若令V0=1600,1800,2000m/s,Vrms=2500m/s,利用式(11)做動校正,可得剩余走時曲線(如圖2d).觀察圖2d可發現,不準確的速度參數V0在產生同相軸的小幅度彎曲的同時,導致了同相軸的局部繞曲.因此,CRP道集中同相軸的局部微彎曲是確定V0的指標,而整體的雙曲彎曲可指示速度參數Vrms的正確與否.實際應用中,同相軸的局部微彎曲可通過疊加道波形的寬窄來判斷,而雙曲彎曲使得可用常規的NMO方法拾取剩余速度.
上述討論表明,基準面上的等效偏移距對確定速度參數是重要的;而三維起伏地表情況,特別是炮、檢點不在同一方位角平面上時,等效偏移距的計算較二維情況變的更復雜.圖3給出了三維情況下等效偏移距的示意圖,圖中S和G點分別代表起伏地表上的炮點和檢波點,AOB平面代表基準面,I代表成像點.采用常規靜校正方法時,地震道的偏移距不變,如圖中的PQ;而本文方法考慮了地震波的實際傳播路徑,基準面上的真實(等效)偏移距應是MN.
已由式(9)或(10)求得由S傳至M再傳至I點的地震波的射線參數pr.對M傳至I點的地震波應用式(10)(可認為T1為零),地震波的射線參數pr不變,則簡單得到


式中(xs,ys)是炮點水平坐標,(xg,yg)是檢波點水平坐標,(x0,y0)是成像點水平坐標,ηs和ηg分別是根據炮點和檢波點信息由表中拾取的η.式(14)表明利用基于表插值的快速算法可在得到走時和幅值的同時,得到等效偏移距;圖3中曲線SMI和GNI代表非均勻介質中地震波的傳播路徑,因此式(14)可適用于非均勻介質.

圖3 等效偏移距示意圖Fig.3 Illustration of equivalent offset on the datum
實際應用中等效參數V0的估計有兩種情況:一是地表調查、微測井、折射或層析靜校正等方法已給出了較準確的近地表速度場;二是尚沒有做上述工作.對于第一種情況,令等效參數V0橫向通過較小的百分比掃描,根據各水平位置處剩余動校后疊加道波形的寬窄和道集的情況綜合確定各水平位置處的最佳百分比參數,對百分比參數橫向平滑,用平滑后的百分比參數修正近地表速度,即得到等效參數V0.對第二種情況,可以估計幾個近地表速度對最小偏移距數據進行偏移,觀察最小偏移距的成像情況確定一個初始的等效參數V0,再將等效參數V0做速度百分比掃描,根據各水平位置處剩余動校后疊加道波形的寬窄和道集情況確定各水平位置處新的V0,橫向平滑后即是最終的近地表速度場.評估疊加道波形的寬窄時應綜合考慮幾個不同深度的同相軸.
等效參數Vrms可根據求得的等效偏移距,在偏移形成的CRP道集上,利用常規的NMO技術拾取.詳見第5節的數據處理流程.
基于式(10)的查表法僅給出炮點至成像點和檢波點至成像點的走時和振幅值;而偏移的目的是進一步得到地下界面的反射系數圖像.若將單個地震道看作是僅有一個接收道的單炮記錄,則根據(4)式和基于(10)式的表插值算法,可分別得到炮域偏移的下行波場和反傳波場:

式中假設震源是一時間脈沖,f(ω)是接收信號的傅里葉變換.將式(15)代入波動方程疊前深度偏移的反褶積成像條件[14],有成像結果

式中F′(t)是f(ω)對應的時域函數的一階導數.式(16)的成像結果表明,對任一單道數據,對數據求一階導數;對成像區域的每一成像點,計算走時ts+tg和成像權系數Ag/As;在數據的一階導數上拾取ts+tg時刻的振幅值并乘上權系數,即完成該道數據的成像.
對全部地震道做上述操作,將成像結果累加,即完成了全部數據的成像.不同于忽略權系數的現行方法,式(16)的權系數Ag/As實現了正確補償地震波的幾何擴散效應,使得深層成像更清晰.
本文偏移流程的關鍵是在偏移過程中確定兩個等效速度V0和Vrms,偏移流程如下:
1.綜合地表起伏情況、可能的地表信息和最小偏移距剖面,確定一個光滑的定浮動基準面;
2.對地表調查、微測井、折射或層析靜校正等方法已給出了較準確的近地表速度場情況,使用該速度場作為初始V0;對無此類信息情況,根據經驗確定一均勻速度作為初始;
3.根據V0,在幾個典型CDP位置,將CMP道集靜校正到浮動基準面,利用NMO確定初始;
4.對未經靜校正處理的數據,依據初始V0和Vrms做式(16)的偏移計算,計算等效偏移距,形成等效偏移距下的共反射點道集;
5.利用初始Vrms,對共反射點道集做反動校,再做NMO動校拾取速度,對拾取的速度做空間平滑,作為新的Vrms;
6.利用初始V0和新的Vrms,對V0做百分比掃描,對較疏的CDP位置和主要的同相軸區域做局部偏移計算,利用第3節討論的方法確定新的V0;
7.在新的V0和Vrms下重新進行偏移計算,形成等效偏移距下的共反射點道集.若剩余動校正量較小,進行剩余動校正、切除和疊加,得到最終偏移剖面;若剩余動校正量較大,利用Vrms做反動校,再做動校,形成再次更新的Vrms,對V0和Vrms重復上述過程,得到最終偏移剖面;
8.由于浮動基準面一般是非水平的,基于浮動基準面的偏移成像將扭曲地下構造圖像.計算浮動基準面與選定的水平基準面的垂直走時,用這一時間修正偏移圖像,可得到基于水平基準面的偏移剖面.
為了驗證本文方法的正確性和工業應用的可能性,分別應用本文方法處理了三維起伏地表采集的理論模型數據和我國東部某工區的實際三維起伏地表采集資料,前者證明了本文理論和方法的正確性,后者表明了本文方法具有實際工業應用價值.
三維起伏地表模型數據由東方地球物理公司提供.該數據近地表高程變化明顯,地表不存在明顯的降低速帶,圖4是模擬該數據集用到的速度模型在inline和crossline方向的切片,實線表示地表高程,最大高差達到590m,虛線是本文采用的浮動基準面.該理論數據集共3940炮,每炮10條線接收,每條線201道,每炮共2010道,中間放炮,兩邊接收.inline方向道間距40m,crossline方向線距100m,最小偏移距50m,最大偏移距5100m,記錄總長度為5s,采樣率為4ms.圖5是該三維數據的典型單炮記錄,從中可以看到:炮記錄中的波場已被嚴重扭曲,同相軸不再是水平地表情況下的雙曲線.圖6給出了使用本文方法得到的最終偏移結果.5個界面的反射同相軸被很好的成像,成像結果與速度模型在形狀上不完全對應是因為成像結果是時間剖面,這也是時間偏移的特點.由于利用了成像權系數,不用AGC深層構造也得到了較好的成像;第5個同相軸較弱是該層的速度對比很小(見圖4速度模型),陡界面成像不好是因為地震記錄的偏移距較小,導致不能接收到陡界面的反射信號所致.

圖4 三維理論數據速度模型,地表高程(實線),浮動基準面(虛線)Fig.4 Velocity model,surface elevation(solid)and floating datum(dot)of 3Dsynthetic data

圖5 三維理論數據典型單炮記錄Fig.5 Typical shot record of 3Dsynthetic data

圖6 三維理論數據直接疊前時間偏移剖面Fig.6 Direct prestack time migration result of the 3Dsynthetic data
該理論模型盡管簡單,但由于地表不存在明顯的降低速帶,應用常規的靜校正加浮動基準面方法是不能正確成像的.
本節進一步給出了應用本文方法處理某工區三維實際數據的例子.圖7是該工區地表高程的等值線圖,可見地表起伏變化是較劇烈的.三維資料共1710炮,每炮16線接收,每線142道,道距40m,線距120m,中間放炮,記錄長度為6s,采樣率為4ms.該三維資料經過壓制面波等去噪處理,但沒有進行靜校正處理.圖8顯示了典型的單炮記錄,地表起伏導致的同相軸錯動明顯可見.圖9給出了在兩個等效速度參數進行更新過程中,典型共反射點道集中同相軸的變化情況,隨著速度參數的更新,同相軸變得更加平直和清晰.

圖7 三維實際數據地表高程Fig.7 Surface elevation of the 3Dfield data
我們沒有得到該工區地表速度的任何資料,因此憑經驗選取1800m/s作為初始.圖10是更新前后在一條inline線上的速度對比情況,圖11是最終得到的在一條inline線上的分布圖.圖12是在一條inline線上的偏移疊加剖面,同相軸連續性較好,斷層斷點清晰,原始資料中的同相軸錯動被消除,這表明本文方法在偏移過程中正確地處理了起伏地表影響,實現了直接偏移.
本文偏移處理是在沒有任何近地表速度資料的情況下進行的,這表明本文方法具有了自主建立近地表和上覆層速度的能力.

圖8 典型三維單炮記錄中單條測線的接收信號Fig.8 Seismic signature of one receiver line with respect to 3Dfield data

圖9 CRP道集隨著速度更新的變化情況Fig.9 Changes in CRP gather while updating the velocities

圖10 Inline方向速度參數V0更新情況Fig.10 Comparison between initial and updated V0parameter in inline

圖11 Inline方向速度參數Vrms更新結果Fig.11 Updated parameter Vrmsin inline
本文提出了一種可對起伏地表采集的三維地震資料直接進行偏移成像的疊前時間偏移方法和流程.該方法可適應高速層出露等不存在明顯低、降速帶的情況,能在偏移過程中自主建立(或修正)近地表速度模型和偏移使用的疊加速度模型.三維理論模型數據證明了方法的正確性和其處理近地表不存在明顯低、降速帶情況的能力.未經靜校正處理的三維實際陸上資料的直接偏移結果表明,本文方法具有很好的工業化應用能力.即使不考慮高速層出露,僅從避免了初至拾取、近地表速度建模等繁雜的靜校正流程而言,本文方法就有很好的實用價值.本文方法為我國陸上油氣勘探面臨的復雜地表問題指示了一個具有工業應用價值的有效途徑.

圖12 實際數據直接疊前時間偏移疊加剖面的inline切片Fig.12 Stack section of direct prestack time migration in inline of 3Dfield data
(References)
[1] Sun C W,Jing P G,Pu Y.Seismic imaging in the mountainous area of Southern China:Statics correction compared with migration from topography.SEG expanded abstract,2009.
[2] 劉國峰,劉洪,李博等.山地地震資料疊前時間偏移方法及其GPU實現.地球物理學報,2009,52(12):3101-3108.Liu G F,Liu H,Li B,et al.Method of prestack time migration of seismic data of mountainous regions and its GPU implementation.Chinese J Geophys.(in Chinese),2009,52(12):3101-3108.
[3] Cox M.Static Corrections for Seismic Reflection Surveys.Tulsa:Society of Exploration Geophysics,1999.
[4] Bevc D.Flooding the topography:Wave-equation datuming of land data with rugged acquisition topography.Geophysics,1997,62(5):1558-1569.
[5] Kelamis P G,Erickson K E,Verschuur D J,Berkhout A J.Velocity-independent redatuming:A new approach to the near-surface problem in land seismic data processing.The Leading Edge,2002,21(8):730-735.
[6] Black J L,Brzostowski M A.Systematics of time-migration errors.Geophysics,1994,59(9):1419-1434.
[7] Reshef M.Depth migration from irregular surfaces with depth extrapolation methods.Geophysics,1991,56(1):119-122.
[8] Rajasekaran S,McMechan G A.Prestack processing of land data with complex topography.Geophysics,1995,60(6):1875-1886.
[9] Wiggins J W.Kirchhoff integral extrapolation and migration of nonplanar data.Geophysics,1984,49(8):1239-1248.
[10] 董春暉,張劍鋒.起伏地表下的直接疊前時間偏移.地球物理學報,2009,52(1):239-244.Dong C H,Zhang J F.Prestack time migration including surface topography.Chinese J.Geophys.(in Chinese),2009,52(1):239-244.
[11] 李世雄.波動方程的高頻近似與辛幾何.北京:石油工業出版社,2001.Li S X.Wave Equation with High Frequency Approximation and Symplectic Geometry(in Chinese).Beijing:Petroleum Industry Press,2001.
[12] Zhang J F,Verschuur D J,Wapenaar C P A.Depth migration of shot records in heterogeneous,transversely isotropic media using optimum explicit operators.Geophysical Prospecting,2001,49(3):287-299.
[13] 張廉萍,劉洪,李幼銘.單程波李代數深度積分的精度分析和算法改進.地球物理學報,2010,53(11):2739-2746.Zhang L P,Liu H,Li Y M.The precision analysis and algorithm improvement in Lie algebra depth integral of oneway seismic wave.Chinese J.Geophys.(in Chinese),2010,53(11):2739-2746.
[14] 張宇.振幅保真的單程波方程偏移理論.地球物理學報,2006,49(5):1410-1430.Zhang Y.The theory of true amplitude oneway wave equation migration.Chinese J.Geophys.(in Chinese),2006,49(5):1410-1430.
3Dprestack time migration including surface topography
ZHANG Hao1,2,ZHANG Jian-Feng1
1 Key Laboratory of Petroleun Resources Research,Institute of Geology and Geophysics,Chinese Academy of Sciences,Beijing100029,China 2 Graduate University,Chinese Academy of Sciences,Beijing100049,China
We present a 3Dprestack time migration scheme and workflow that can image land seismic data without datum static corrections.The scheme mimics the influences of the nearsurface and overlay on wave propagation using two equivalent velocities.As a result,it is not necessary for wave path to be vertical incidence in the near-surface,as in conventional static corrections.This leads to the proposed scheme which can handle the more complex surface cases,such as the exposure of high velocity in the near-surface.The proposed scheme can update the two equivalent velocities with respect to the flats and consistency of the events.A synthetic 3D dataset and a field 3Ddataset without datum static corrections are used to demonstrate the scheme and the workflow.
Topography,3Dimaging,Near-surface velocity,Pre-stack time migration
P631收修定稿2010-12-16,2012-01-17收修定稿
國家自然科學基金重點項目(40930422)和國家科技重大專項(2011ZX05008-006)資助.
張浩,男,1983年生,中國科學院地質與地球物理研究所博士研究生,主要從事疊前偏移成像研究.E-mail:zhanghao@mail.iggcas.ac.cn
張浩,張劍鋒.起伏地表采集數據的三維直接疊前時間偏移方法.地球物理學報,2012,55(4):1335-1344,
10.6038/j.issn.0001-5733.2012.04.029.
Zhang H,Zhang J F.3Dprestack time migration including surface topography.Chinese J.Geophys.(in Chinese),2012,55(4):1335-1344,doi:10.6038/j.issn.0001-5733.2012.04.029.
10.6038/j.issn.0001-5733.2012.04.029
(本文編輯 汪海英)