葛利華,姜 弢*,林 君,徐學純,黃大年,賈海青
1 教育部地球信息探測儀器重點實驗室,長春 130026
2 吉林大學儀器科學與電氣工程學院,長春 130026
3 吉林大學地球科學學院,長春 130061
4 吉林大學地球探測科學與技術學院,長春 130026
深部礦產資源勘探過程中,由于勘探深度大,地面接收到的信號經常會被噪聲淹沒,為提高有效信號的能量,首先需要解決的問題就是如何合理選擇震源.目前用于地震勘探的震源主要分為兩種:一種為炸藥震源,另一種為非炸藥震源.炸藥震源,以其能量強、頻帶寬、分辨率高等優點,在野外油氣勘探中得到廣泛應用,但存在危險、不可控,震源間的一致性差,且當藥量增加到一定程度時激發的地震波能量將不再相應增加[1]等問題.常用的非炸藥震源有電火花震源、石油勘探用車載液壓震源、電磁驅動式可控震源[2].電火花震源適于海上地震勘探.石油勘探用車載液壓震源,不適用于礦區復雜地形條件.輕便電磁驅動式可控震源系統,安全、可控且便攜,有利于復雜地形地震勘探.但是,單個電磁驅動可控震源的輸出力較小,難以獲得地下深部信息.為此,通常使用多個震源組合激發方式,但組合震源主要能量方向垂直向下,在對復雜礦區進行勘探時,并不能有效提高地下目標數據的信噪比.為滿足深部礦產勘探的需求,相控震源技術逐漸受到人們關注.
相控陣列的思想來自于雷達天線系統[3],利用相控陣列的方向特性,使得能量在某一特定方向匯聚,以達到遠距離探測的目的.Arnold[4]首次將相控陣列應用到地震勘探,利用多臺車載液壓震源,進行了相控震源定向波束的野外試驗,驗證了相控陣列在地震勘探中應用的可行性.由于復雜地形礦區的地理因素,考慮用輕便的電磁驅動可控震源代替車載液壓震源,進行山區作業.基于以上考慮,本文重點對電磁驅動可控震源相控陣列在不同礦區模型的地震響應進行數值模擬研究.
地震波場正演模擬是合理解釋地震數據資料的前提,較常用的方法有射線追蹤法、波動方程數值解法和積分方程法.射線追蹤法是建立在高頻近似基礎上的一種地震數值模擬方法[5],用幾何射線的方法解釋波的運動學特征.2004年,吉林大學可控震源課題組[6-8]用幾何射線追蹤法對電磁驅動可控震源相控陣列的方向特性進行了分析和模擬,得到基于電磁驅動可控震源的相控陣列,通過合理設置參數,可以在感興趣的方向上有效提高地震信號的信噪比.基于幾何射線追蹤法的地震波數值模擬,雖然運算速度快,顯示直觀,但只反映了地震波的運動學特征,缺乏對地震波動力學特征的描述,且僅適用于簡單模型.積分方程法基于Huygens原理,建立在波動方程積分表達式的基礎上,計算效率低,主要針對復雜邊界區域的求解.本文采用波動方程數值解法.2006年,王忠仁等[9]用有限差分的方法,對相控陣震源在相關前Chirp信號激勵下的地震響應進行了數值模擬,引入地震波場邊際能量密度的概念,利用地震波場的時間切片技術,對模型空間各個方向上的能量強度進行了定量分析,描述了均勻介質模型下地震波的方向特性.
本文基于電磁驅動可控震源相控技術,采用有限差分方法求解波動方程,仿真計算地下二維空間分別采用單震源、9單元相控源激勵下的波場傳播及地面接收位置的地震響應,給出了基于均勻介質、水平層狀礦區和不規則礦區三種模型下產生的波場快照和波傳播的方向特性圖,通過對比,定量得到相控震源在不同介質模型下激發地震波的方向性特點及對目標數據信噪比的改善能力.
相控震源,借鑒相控天線[10]波束定向原理,通過設定陣列[11]各單元起震的絕對時間或相對時間差,以實現各震動單元產生的波場在某個特定方向上的同相疊加,提高了這一方向上的信號能量,合理設置檢波器的位置,能有效提高地震記錄的信噪比.這一技術已被應用到地震測井等領域[12-14].
相控震源地震勘探原理如圖1所示.

圖1 相控震源系統野外工作原理圖Fig.1 Principle of Phased-array vibrator system
在均勻介質中,速度v為常數,圖2中方向角θ定義為地下目標方向與豎直向下方向夾角,且逆時針方向為正方向,則當相控震源相鄰激震器間距d和延時參數τ一定時,有如下公式[15]


圖2 τ=0.5ms,d=2m,v=2000m/s時主波束方向圖Fig.2 Main-beam direction graph whileτ=0.5ms,d=2m,v=2000m/s
由公式(1),可以計算地震波主波束在均勻介質中的傳播方向,如圖2.
圖2結果說明,相控源激發的地震波在均勻介質中存在方向性,但非均勻介質相控震源波束形態不能確定,需通過求解波動方程獲得.
波動方程數值解法模擬地震響應,主要分為有限差分法[16-17]、有限元法[18-20]和偽譜法[21-22]等.有限元法計算速度慢,占內存大,偽譜法難以消除人工邊界.本文采用計算速度快,易于計算機實現的有限差分法.通過有限差分求解波動方程的方法,模擬相控震源在不規則礦區模型中的地震響應.主要方法如下,首先求解相控震源的每個單元形成的地震波場,然后將各相控單元的波場延時疊加,即得到相控震源的地震波場.
在求解每個單元激勵的地震波場時,需要構建波動方程,包括初始條件、邊界條件以及其它限定條件;波動方程中的震源函數需結合可控震源控制信號特點及可控震源地震勘探原理導出;討論對波動方程進行數值離散后,差分方程的約束條件和參數選取.
首先討論相控震源的單個震動單元波場求解.在非均勻介質中,設定波阻抗為常數,則有全聲波方程:

在模擬區域內,為消除層間多次波,在公式(2)的右側加入補償項[23]

考慮初始時刻,區域內的波場值為0,地震波的傳播速度為0,用方程的形式表示為

在區域的邊界上,將人工截斷邊界處理為吸收邊界[24],此種邊界條件要求震源不能離邊界太近.

另外,對于地面和反射界面,有如下約束方程:

聯立公式(3)—(10)建立波動方程組,求解模擬區域各時刻的波場值.其中,u為網格節點的波場值,D為需求解的區域,g為地面,gr為反射界面,r為界面反射系數,f(t)為震源函數,v為縱波速度,s為常數.
上述求解單個震動單元波場過程,需要給出震源函數.震源函數求解與可控震源激勵信號及可控震源地震勘探原理密切相關.
電磁驅動式可控震源激勵的控制信號X(t)一般采用如公式(11)所示的Chirp信號形式[2]:

其中,A為掃描信號的幅度,f1,f2為掃描起始頻率和終止頻率,t為掃描時間,T為掃描時長.
為加大勘探深度,震源宜采用頻率較低的信號,這里以10~100Hz Chirp為例表示,如圖3a.
可控震源與炸藥震源工作原理不同.炸藥震源工作時,震源信號為沖激信號形式,記為δ(t),當把大地看成一個線性系統,其沖擊響應函數用h(t)表示,則檢波器接收到的地震信號

其中,*為卷積運算符號.
而可控震源工作時,震源控制信號X(t)為調頻信號,檢波器接收到的信號

要得到與炸藥震源相似的具有正確到時信息的脈沖形式地震信號,由可控震源地震原理,可控震源地震數據必須經過相關檢測處理[25],相關檢測后得到的地震數據可表示為

圖3 由可控震源控制信號構造震源地震子波(a)Chirp形式的震源控制信號;(b)控制信號自相關;(c)用于數值模擬的可控震源地震子波.Fig.3 Signal of source(a)Control signal of source in the form of chirp;(b)Control signal autocorrelation;(c)Vibroseis wavelet for numerical simulation.

其中,? 為相關運算符號.X(t)?X(t)為X(t)的自相關函數,如圖3b.對比公式(12)和(14),可控震源的等效震源函數應為X(t)?X(t),即控制信號的自相關函數.容易看出自相關函數主要的能量在中心點位置附近,離中心點位置較遠部分能量很弱,接近為0.計算表明,在自相關函數中心點左右各取32ms時窗內的能量占整個自相關函數96.4%的能量,為計算方便,本文將波動方程中的震源函數定義為

如圖3c,其中,Δt=32ms.
上面討論了相控震源的單個單元震源函數的求取過程.針對相控震源的不同震動單元,可在相應位置加入震源函數.
有限差分法是將求解區域網格離散化,對離散后的網格節點進行求解,用有限的節點代替連續的求解區域[26].構造差分一般基于Tailor公式展開,用差商代替微商,構建差分方程,然后對差分方程進行求解.本文采用空間域四階,時間域二階的差分格式.
計算中,網格剖分越細,數值解就會越趨近于解析解,但同時也會增大計算量.在模型已知的情況下,如何選擇合適的網格步長,既能滿足要求,又能提高計算效率,是有限差分法解波動方程需要考慮的首要問題.


另外,用有限差分法求解波動方程,時間和空間步長都趨近于零時,差分方程的解趨近于微分方程的解,但也要考慮到計算量問題.一般情況下,當差分方程滿足穩定性條件時,差分方程就收斂于微分方程.鑒于本文采用的空間域四階差分,時間域二階差分,應滿足如下穩定性[23]條件:


根據采樣定理及經驗,一般時間采樣率要求不小于8~10倍的fmax,以及式(16)和(17)的條件約束,模擬中選取空間采樣步長dh=4m,時間采樣步長dt=1/1500≈0.00067s.
以上討論了相控震源的每個單元地震波場的求取過程,根據相控震源原理,將各相控單元的波場延時疊加,即可得到相控震源地震波場.
針對單震源和相控源,分別給出波場快照、方向特性圖,地表檢波器信噪比改善等波場模擬結果.
為展示相控震源具有激發任意方向地震波的能力,相控源選擇了0和1.33ms兩個延時參數.相控震源采用0ms延時,實際相當于組合震源工作方式;延時參數1.33ms為任意選擇,只用于說明延時參數變化對地震波方向的影響.
設計三種模型,即均勻介質模型,水平層狀礦區模型和不規則礦區模型,如圖4所示,其中黑色圓點標記為震源點位置.
計算區域長1020m,深1020m,取水平向右為x軸的正方向,豎直向下為z軸正方向,兩方向的交叉點取為坐標原點,建立直角坐標系.相控震源各單元炮點坐標分別為(168,20)、(176,20)、(184,20)、(192,20)、(200,20)、(208,20)、(216,20)、(224,20)、(232,20),為作對比,將單震源的炮點設在震源陣列的中間位置,坐標為(200,20).
由圖5可見,均勻介質模型中,單震源產生的地震波不存在方向性;相控源激勵的地震波具有明顯方向性,延時參數不同,波束方向不同.
為了描述相控源的方向性,引入方向特性圖.方向特性圖定義為0~0.7s內每隔6.7ms各網格節點處波場能量值的累加.并且規定地震波主波束方向與x軸正向的夾角為方向特性圖中的方向角,記作η.

圖4 (a)均勻介質模型;(b)水平層狀礦區模型;(c)不規則礦區模型Fig.4 (a)Model of homogeneous media;(b)Model of horizontal layers;(c)Model of irregular mining

圖5 單震源、相控源0,1.33ms延時參數下均勻介質模型的波場快照(a)133ms;(b)233ms.Fig.5 Snapshots of seismic wave of the single vibrator,0ms and 1.33ms phased-array vibrator in model of homogeneous media(a)133ms;(b)233ms.

圖6 均勻介質模型方向特性圖(a)單震源;(b)9單元組合源;(c)9單元相控源.Fig.6 Directional character of model of homogeneous media(a)Single vibrator;(b)9-unit vibrator combination;(c)9-unit vibrator phased-array.
由圖6可見,均勻介質模型中,單震源不能產生定向地震波束;相控源激勵下的地震波束有明顯的方向性,方向因延時參數的不同而發生變化,0ms地震震波方向垂直向下,即η為90°,1.33ms波束方向η為70.4°;0ms和1.33ms相控震源分別提高了90°和70.4°方向上信號的能量.
圖7(a-d)分別為地震波在水平層狀礦區模型的不同地層中傳播時的波場快照.由圖7可見,單震源產生的地震波不存在方向性;相控源激勵下的地震波有了明顯的方向性,延時參數不同,方向不同;0ms地震波方向垂直向下,且不同地層中方向不變;1.33ms地震波方向在不同地層中發生偏折;由圖7c見,0ms相控源激勵下第一層的反射波與直達波混雜,降低了分辨率,而1.33ms相控源激勵下第一層的反射波與直達波能明顯分開,提高了分辨率,圖7d中,其它層的反射波也以1.33ms相控源的能量和分辨率最高.
由圖8的方向特性圖可見,水平層狀礦區模型中,單震源不能產生定向的地震波束;0ms相控源產生的地震波主波束方向垂直向下,即方向角η為90°;1.33ms相控源能產生定向的地震波束,但方向隨著介質中速度的變化,主波束方向發生偏折,地震波主波束方向角η由70.4°,依次變為68.5°,66.0°,63.7°,61.1°.

圖7 單震源、相控源0、1.33ms延時參數下(從左至右)水平層狀礦區模型波場快照(A)133ms;(B)233ms;(C)333ms;(D)367msFig.7 Snapshots of seismic wave of the single vibrator,0ms and 1.33ms phased-array vibrator in model of horizontal layers(A)133ms;(B)233ms;(C)333ms;(D)367ms.
圖9為不規則礦區模型的波場快照.圖9a為地震波穿過第一層介質,進入圍巖,圖9b為地震波開始進入不規則礦體,圖9c為地震波下行波波前即將傳出不規則礦體.由圖9(b,c)可見,單震源激勵下的地震波不存在方向性;相控源激勵下的地震波存在方向性,延時參數不同,方向不同;0ms地震波方向垂直向下,且不同地層中方向不變;1.33ms地震波存在方向性,方向在不同地層中發生偏折;由圖9(b,c)可見,不規則礦體上界面處,單震源激勵下的地震波能量弱,不能看到明顯的反射波,0ms相控源能看到反射波,但分辨率低,1.33ms相控源反射波能量加強,分辨率高.

圖8 水平層狀礦區模型方向特性圖(a)單震源;(b)0ms相控源;(c)1.33ms相控源.Fig.8 Directional character of model of horizontal layers(a)Single vibrator;(b)0ms phased-array vibrator;(c)1.33ms phased-array vibrator.

圖9 單震源、相控源0、1.33ms延時參數下(從左至右)不規則礦區模型波場快照(A)133ms;(B)333ms;(C)367ms.Fig.9 Snapshots of seismic wave of the single vibrator,0ms and 1.33ms phased-array vibrator in model of irregular mining(A)133ms;(B)333ms;(C)367ms.

圖10 不規則礦區模型方向特性圖(a)單震源;(b)0ms相控源;(c)1.33ms相控源.Fig.10 Directional character of model of irregular mining(a)Single vibrator;(b)0ms phased-array vibrator;(c)1.33ms phased-array vibrator.
由圖10的方向特性圖可見,不規則礦區模型中,單震源不能產生定向的地震波束;相控源激勵下的地震波束有明顯的方向性,方向因延時參數的不同而發生變化;0ms相控源產生的地震波束垂直向下,即方向角η為90°,且在不同地層中方向不變;1.33ms相控源產生的定向地震波束方向角η在首層介質內為70.4°,圍巖內為65.4°,礦體內為63.1°;相比于單震源,0ms和1.33ms相控源在不規則礦體的分界面處有更高的分辨率,三者中1.33ms相控源的最高.
在噪聲強度相同的情況下,設兩個信號S1,S2,則后者相對于前者信噪比的改善[27]為

由公式(18)可以定量分析相控源相對單個可控震源在接收點處信噪比的改善能力.需要注意的是,地表檢波器排列位置不同,接收數據的信噪比也不相同.
圖11a給出了水平層狀礦區模型中,在地面888m處檢波器接收到第一層界面的反射波波形.由公式(18)計算得出,相對于單震源,1.33ms和0ms相控源信噪比分別改善了14.1dB和10.2dB.
圖11b給出了不規則礦區模型中,地面244m處檢波器接收到的不規則礦體上分界面的反射波波形.相對于單震源,1.33ms和0ms相控源信噪比分別改善了18.9dB和14.9dB.
上述結果表明,在本文給出的兩種延時參數下,相控源均具有改善信噪比的能力.實際上,由于延時參數直接影響了地震波場在不同方向上的能量分布,可以預期,很多個但不是任一延時參數均能夠改善地下目標信噪比.理論上來說,針對某個檢波點或檢波器排列,是存在一個最優的延時參數,可以令定向地震波照射到地下目標體,使得檢波器接收的地下目標數據具有最優信噪比.由于延時參數最優選取是一個較復雜的問題,其涉及內容較多,本文暫不討論.

圖11 單震源、0、1.33ms相控源激勵下檢波器分別接收到的反射波(a)888m水平層狀礦區模型;(b)244m不規則礦區模型Fig.11 Reflected wave recorded in the receiving point for the single vibrator,0ms and 1.33ms phased-array vibrator(a)Model of horizontal layers,888m;(b)Model of irregular mining,244m.
根據模擬結果,本文得到兩個結論.一是相控震源在復雜介質條件下仍然具有激發定向地震波的能力,除0延時參數外,相控震源形成的地震波主波束方向會隨介質變化發生偏折;二是不同延時參數下,目標方向的地震波能量分布發生變化,可能存在的多個延時參數均能改善目標數據信噪比,但信噪比改善能力不同.以上結論說明,相控震源技術應用于復雜礦產資源勘探是可行的而且有益的.
同時,我們也清楚地看到,如果能設計合適的延時參數,令地下目標方向上地震波能量最強,可以最大程度改善采集數據質量.有關最優延時參數設計是一個較復雜的問題,理論上存在最優延時參數,但針對復雜模型或實際勘探環境,延時參數選擇與多種因素有關,包括相控震源工作參數、地下介質模型及關注的主要目標區域、檢波器排列位置等因素,這部分內容研究將作為后續研究工作的重點.另外,相控震源波場模擬,對相控震源地震采集過程的觀測系統設計也具有重要指導作用.
(References)
[1]張智,劉財,邵志剛.地震勘探中的炸藥震源藥量理論與實驗分析.地球物理學進展,2003,18(4):724-728.Zhang Z,Liu C,Shao Z G.Theory and experimentation of the charge sizes in seismic-sources for seismic exploration seismic-sources for seismic exploration.Progress in Geophysics(in Chinese),2003,18(4):724-728.
[2]林君.電磁驅動可控震源地震勘探原理及應用.北京:科學出版社,2004.Lin J.Seismic Exploration and Application of Vibrator Driven by Electromagnetic (in Chinese).Beijing:Science Press,2004.
[3]Garrod A.Digital modules for phase array radar.//IEEE International Radar Conference.Alexandria,VA:IEEE,1995:726-731.
[4]Arnold M E.Beam forming with vibrator arrays.Geophysics,1977,42(7):1321-1338.
[5]裴正林,牟永光.地震波傳播數值模擬.地球物理學進展,2004,19(4):933-941.Pei Z L,Mu Y G.Numerical simulation of seismic wave propagation.Progress in Geophysics (in Chinese),19(4):933-941.
[6]姜弢,林君,楊冬等.相控震源定向地震波信號分析.地球物理學報,2008,51(5):1551-1556.Jiang T,Lin J,Yang D,et al.Analysis of directional seismic signal based on phased-array vibrator system.Chinese J.Geophys.(in Chinese),2008,51(5):1551-1556.
[7]姜弢,林君,陳祖斌等.相控震源對水平層狀地下介質的高信噪比檢測.儀器儀表學報,2006,27(11):1396-1372.Jiang T,Lin J,Chen Z B,et al.High signal-to-noise ratio detection of horizontal layer underground medium with phased-array vibroseis. Chinese Journal of Scientific Instrument(in Chinese),2006,27(11):1369-1372.
[8]姜弢,林君,陳祖斌等.相控震源地震波定向技術.吉林大學學報(信息科學版),2004,22(3):181-184.Jiang T,Lin J,Chen Z B,et al.Seismic beam-forming in phased array of vibroseis.Journal of Jilin University(Information Science Edition)(in Chinese),2004,22(3):181-184.
[9]王忠仁,林君,馮聲涯.地震勘探中相控陣震源的方向特性研究.地球物理學報,2006,49(7):1191-1197.Wang Z R,Lin J,Feng S Y.Directivity of phase array vibrators in seismic exploration.Chinese J.Geophys.(in Chinese),2006,49(4):1191-1197.
[10]Konishi Y.Phased array antennas.IEICE Transactions on Communications,2003,E86B(3):954-967.
[11]朱崇燦.天線.武漢:武漢大學出版社,1996.Zhu C C.Antenna(in Chinese).Wuhan:Wuhan University Press,1996.
[12]喬文孝,車小花,鞠曉東等.聲波測井相控圓弧陣及其輻射指向性.地球物理學報,2008,51(3):939-946.Qiao W X,Che X H,Ju X D,et al.Acoustic logging phased arc array and its radiation directivity.Chinese J.Geophys.(in Chinese),2008,51(3):939-946.
[13]Che X H,Qiao W X.Numerical simulation of an acoustic field generated by aphased arc array in a fluid-filled borehole.Petroleum Science,2009,6(3):225-229.
[14]喬文孝,杜光升,陳雪蓮.相控線陣聲波輻射器在聲波測井中應用的可行性分析.地球物理學報,2002,45(5):714-722.Qiao W X,Du G S,Chen X L.Feasibility of application of phased array in acoustic well-logging.Chinese J.Geophys.(in Chinese),2002,45(5):714-722.
[15]姜弢.基于相控震源的地震波定向方法研究[博士論文].吉林:吉林大學儀器科學與電氣工程學院,2006.Jiang T.Study on seismic beam-forming method based on phased-array vibrator system[Ph.D thesis](in Chinese).Jilin:College of Instrumentation Electrical Engineering,2006.
[16]Mohanty R K,Venu G.High accuracy cubic spline finite difference approximation for the solution of one-space dimensional non-linear wave equations.Applied Mathematics and Computation,2011,218(8):4234-4244.
[17]Liu Y,Sen M K.Numerical modeling of wave equation by a truncated high-order finite-difference method.Earthq.Sci.,2009,22(2):205-213.
[18]Phongthanapanich S,Dechaumphai P.An explicit finite volume element method for solving characteristic level set equation on triangular grids.Acta Mechanica Sinica,2011,27(6):911-921.
[19]RenéM.An efficient finite element time-domain formulation for the elastic second-order wave equation: Anon-split complex frequency shifted convolutional PML.Int.J.Numer.Meth.Engng.,2011,88(10):951-973.
[20]Kahkomen S,Glowinski R,Rossi T,et al.Solution of Time-Periodic wave equation using mixed finite elements and controllability techniques.Journal of Computational Acoustics,2011,19(4):335-352.
[21]Carlos S,Escolano J,Garriga A.Semi-empirical boundary conditions for the linearized acoustic Euler equations using Pseudo-Spectral Time-Domain methods.Applied Acoustics,2011,72(4):226-230.
[22]李淅龍.煤礦井下反射地震勘探技術初步研究與應用[碩士論文].西安:西安大學,2010.Li X L.The Preliminary Study Andapplication of Reflection Seismic Exploration Technique Applied in Coal Mine Underground[Master thesis](in Chinese).Xi′an:Xi′an University of Science and Technology,2010.
[23]單延明,吳清嶺,于承業等.復雜地質構造波動方程數值模擬.物探化探計算技術,2002,24(1):22-26.Shan Y M,Wu Q L,Yu C Y,et al.Wave equation modeling for complex media.Computing Techniques for Geophysical and Geochemical Exploration (in Chinese),2002,24(1):22-26.
[24]Reynolds A C.Boundary conditions f or the numerical solution of wave propagation problems.Geophysics,1978,43(6):1099-1110.
[25]姜弢,林君,陳祖斌等.可控震源相控地震的相關檢測技術.儀器儀表學報,2005,26(4):336-339.Jiang T,Lin J,Chen Z B,et al.The correlation detection method of signals in phased array of vibroseis system.Chinese Journal of Scientific Instrument (in Chinese),2005,26(4):336-339.
[26]楊瑩.二維地震波場有限差分法數值模擬研究[碩士論文].北京:中國地質大學,2009.Yang Y.The research on Two-dimensional Finite-difference seismic wave field numerical simulation[Master thesis](in Chinese).Beijing:China University of Geosciences,2009.
[27]姜弢,林君,李桐林等.相控震源對地震信號信噪比的改善研究.地球物理學報,2006,49(6):1819-1825.Jiang T,Lin J,Li T L,et al.Boosting signal to noise ratio of seismic signals using the phased-array vibrator system.Chinese J.Geophys.(in Chinese),2006,49(6):1819-1825.