雷林林, 劉四新, 傅磊*, 吳俊軍

1 吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院, 長(zhǎng)春 130026 2 中國(guó)石油集團(tuán)東方地球物理勘探有限責(zé)任公司 新興物探開(kāi)發(fā)處， 河北涿州 072751

?

基于全波形反演的探地雷達(dá)數(shù)據(jù)逆時(shí)偏移成像

雷林林1, 劉四新1, 傅磊1*, 吳俊軍2

1 吉林大學(xué)地球探測(cè)科學(xué)與技術(shù)學(xué)院, 長(zhǎng)春 130026 2 中國(guó)石油集團(tuán)東方地球物理勘探有限責(zé)任公司 新興物探開(kāi)發(fā)處， 河北涿州 072751

逆時(shí)偏移成像(RTM)常用來(lái)處理復(fù)雜速度模型，包括陡傾角及橫向速度變化劇烈的模型.與常規(guī)偏移成像方法(如Kirchhoff偏移)相比，逆時(shí)偏移成像能提供更好的偏移成像結(jié)果，近些年逆時(shí)偏移成像越來(lái)越廣泛地應(yīng)用到勘探地震中，它逐漸成為石油地震勘探中的一種行業(yè)標(biāo)準(zhǔn).電磁波和彈性波在動(dòng)力學(xué)和運(yùn)動(dòng)學(xué)上存在相似性，故本文開(kāi)發(fā)了基于麥克斯韋方程組的電磁波逆時(shí)偏移成像算法，并將其應(yīng)用到探地雷達(dá)數(shù)據(jù)處理中.時(shí)間域有限差分(FDTD)用于模擬電磁波正向和逆向傳播過(guò)程，互相關(guān)成像條件用于獲得最終偏移結(jié)果.逆時(shí)偏移成像算法中，偏移成像結(jié)果受初始模型影響較大，而其中決定電磁波傳播速度的介電常數(shù)的影響尤為重要.本文基于時(shí)間域全波形反演(FWI)算法反演獲得了更為精確的地下介電常數(shù)模型，并將其反演結(jié)果作為逆時(shí)偏移成像的初始介電常數(shù)模型.為了驗(yàn)證此算法的有效性，首先構(gòu)建了一個(gè)復(fù)雜地質(zhì)結(jié)構(gòu)模型，合成了共偏移距及共炮點(diǎn)探地雷達(dá)數(shù)據(jù)，分別應(yīng)用常規(guī)Kirchhoff偏移算法及逆時(shí)偏移成像算法進(jìn)行偏移處理，成像結(jié)果顯示由逆時(shí)偏移成像算法得到的偏移結(jié)果與實(shí)際模型具有較高的一致性；此外本文在室內(nèi)沙槽中進(jìn)行了相關(guān)的物理模擬實(shí)驗(yàn)，采集了共偏移距及共炮點(diǎn)探地雷達(dá)數(shù)據(jù)，分別應(yīng)用Kirchhoff和疊前逆時(shí)偏移成像算法進(jìn)行處理，結(jié)果表明疊前逆時(shí)偏移成像在實(shí)際應(yīng)用中能獲得更好的成像效果.

逆時(shí)偏移； 探地雷達(dá)； 全波形反演； 物理模擬實(shí)驗(yàn)

1 引言

偏移是一種地震成像技術(shù)，它在空間域或時(shí)間域?qū)⒂涗浀降耐蜉S歸位到真實(shí)的反射界面.傳統(tǒng)的地震偏移成像方法基于單程波方程，并將多次散射波當(dāng)做噪聲處理，但對(duì)于復(fù)雜地質(zhì)結(jié)構(gòu)及大傾度界面如裂縫，合理的利用多次散射信息至關(guān)重要.傳統(tǒng)波動(dòng)方程偏移方法如Kirchhoff積分偏移(Gray, 1986; Gray and May, 1994; Docherty)難以適應(yīng)橫向速度變化；有限差分波動(dòng)方程偏移(Vidale, 1988; Mufti et al., 1996)近似求解波動(dòng)方程，存在傾角限制；頻率波數(shù)域偏移方法(Loewenthal and Mufti, 1983)是基于常速度模型，當(dāng)?shù)叵陆橘|(zhì)速度存在變化時(shí)會(huì)出現(xiàn)明顯的精度問(wèn)題.有別于傳統(tǒng)的基于單程波方程的偏移方法，逆時(shí)偏移(Baysal et al., 1983; Levin, 1984; Chang and McMechan, 1986，1994; Sun and McMechanz, 2001; Yoon et al., 2003; 劉欣欣等, 2013; 王保利等, 2012)基于全波方程，能夠精確的處理沿任意方向傳播的波場(chǎng)，包括耦合波及橫縱波之間的轉(zhuǎn)換波，從而可更精準(zhǔn)地將波場(chǎng)歸位到其真實(shí)的地下空間位置.雖然逆時(shí)偏移技術(shù)的成功應(yīng)用實(shí)例并不多見(jiàn)，但它在勘探地震中受到越來(lái)越多的關(guān)注，逆時(shí)偏移包括疊前逆時(shí)偏移和疊后逆時(shí)偏移.

探地雷達(dá)電磁波與地震波在運(yùn)動(dòng)學(xué)和動(dòng)力學(xué)存在一定的相似性，此外，采集到的探地雷達(dá)數(shù)據(jù)和地震數(shù)據(jù)都是空間位置的時(shí)域電壓信號(hào)，故而在一定程度上勘探地震中的數(shù)據(jù)處理技術(shù)可以運(yùn)用到探地雷達(dá)數(shù)據(jù)處理中.Fisher(1992)提出用基于張量波動(dòng)方程的逆時(shí)偏移算法對(duì)單通道的共偏移距的探地雷達(dá)數(shù)據(jù)進(jìn)行了偏移處理，處理后的結(jié)果在一定程度顯示了該算法的優(yōu)越性.Feng和Sato(2004)利用疊前逆時(shí)偏移算法，對(duì)地雷探測(cè)中采集到的多偏移距數(shù)據(jù)進(jìn)行處理，一方面強(qiáng)散射波得到壓制，另一方面獲得了地下更為精細(xì)的雷達(dá)剖面圖像.Zhou等(2005)應(yīng)用一種基于遺傳算法的偏移速度分析方法，對(duì)合成及野外實(shí)測(cè)的共炮點(diǎn)雷達(dá)數(shù)據(jù)進(jìn)行分析獲得了地下介電常數(shù)分布，并以其作為疊前逆時(shí)偏移的初始速度模型，取得了較好的結(jié)果.

本文開(kāi)發(fā)了電磁波逆時(shí)偏移算法，電磁波場(chǎng)的正向及逆向傳播是基于時(shí)間域有限差分方法，單軸各向異性完全匹配層(UPML)(Fang and Ying, 2009; Feng and Dai, 2011)吸收邊界條件用于吸收邊界的反射波，互相關(guān)成像條件用于獲得空間圖像.應(yīng)用電磁波全波形反演算法，對(duì)合成及實(shí)測(cè)的共偏移距、多偏移距探地雷達(dá)數(shù)據(jù)全波形反演，獲得了與實(shí)際情況更為接近的地下介質(zhì)分布情況.通過(guò)全波形反演獲得的介電常數(shù)和電導(dǎo)率模型作為逆時(shí)偏移的初始模型，逆時(shí)偏移的結(jié)果與基于共偏移距雷達(dá)數(shù)據(jù)Kirchhoff偏移成像結(jié)果對(duì)比發(fā)現(xiàn)，逆時(shí)偏移成像效果明顯好于傳統(tǒng)偏移成像結(jié)果.

2 逆時(shí)偏移基本原理及電磁波逆時(shí)偏移算法

傳統(tǒng)的偏移算法基本是基于單程波方程深度外推，所以在對(duì)多分量地震資料處理過(guò)程中首先要進(jìn)行縱、橫波分離，然后用現(xiàn)有的偏移算法分別處理縱波場(chǎng)和橫波場(chǎng).它的缺陷在于：1)將多分量地震資料當(dāng)作幾個(gè)標(biāo)量波的簡(jiǎn)單疊加，忽略了地震波的矢量特性，從而影響地震多分量資料的處理與解釋精度；2)假如地震波場(chǎng)橫波與縱波沒(méi)有完全分離(目前這種現(xiàn)象廣泛存在)，而在實(shí)際資料的偏移處理過(guò)程中又假設(shè)為標(biāo)量波場(chǎng)，這必然會(huì)使偏移結(jié)果產(chǎn)生許多人為假象，影響處理精度.逆時(shí)偏移與傳統(tǒng)的偏移方法有很大的不同，逆時(shí)偏移可以看成是波場(chǎng)在時(shí)間軸上逆向傳播過(guò)程，當(dāng)波場(chǎng)倒退到零時(shí)刻，則所有反射波與繞射波的能量都回到最初被反射和繞射的空間位置，應(yīng)用成像條件即可得到最終的偏移成像剖面.逆時(shí)偏移的波動(dòng)理論方法基于全波方程，波場(chǎng)可以沿各個(gè)方向傳播，因此該方法不受介質(zhì)傾角限制，適用于速度任意變化的模型.

逆時(shí)偏移算法的實(shí)現(xiàn)概括起來(lái)主要包括三個(gè)步驟：第一是震源波場(chǎng)正向延拓；第二是接收波場(chǎng)的逆時(shí)延拓；第三是應(yīng)用成像條件成像.前兩個(gè)步驟的實(shí)現(xiàn)可以基于相同的數(shù)值方法，本文應(yīng)用時(shí)間域有限差分法(劉四新和曾昭發(fā), 2007; 劉四新等, 2006; Taflove and Hagness, 2005; Liu and Sato, 2005; Liu et al., 2007; Liu et al., 2011)實(shí)現(xiàn)電磁波場(chǎng)在空間的正向和逆向傳播.影響疊前逆時(shí)偏移成像效果的一個(gè)非常重要因素是成像條件(Chattopadhyay and McMechan, 2008)的應(yīng)用，Claerbout(1971)提出的互相關(guān)成像條件是最基礎(chǔ)的成像條件，它采用零延遲互相關(guān)成像條件，該成像條件的不足之處在于對(duì)首波及背景散射波做互相關(guān)會(huì)產(chǎn)生低頻噪聲，其次應(yīng)用該成像條件得到的成像結(jié)果與地下介質(zhì)的反射系數(shù)沒(méi)有對(duì)應(yīng)關(guān)系.鑒于零延遲互相關(guān)成像條件的弊端，后續(xù)逐步發(fā)展了反褶積成像條件(Guitton et al., 2007)、震源補(bǔ)償成像條件(Guitton et al., 2006)、擴(kuò)展成像條件(Sava and Fomel, 2006)以及角度域成像條件(Yoon and Marfurt, 2006)等.標(biāo)準(zhǔn)的逆時(shí)偏移采用零延遲互相關(guān)成像條件，其計(jì)算公式為：

(1)

其中IMAG(x)為地下空間x的成像結(jié)果，S(x,t,si)為源波場(chǎng)，R(x,T-t,si)為接收波場(chǎng)，t和si分別代表時(shí)間和源序號(hào).當(dāng)接收波場(chǎng)傳播到時(shí)間t時(shí)，讀取震源波場(chǎng)t時(shí)刻的場(chǎng)值和對(duì)應(yīng)時(shí)刻接收波場(chǎng)的場(chǎng)值并將兩波場(chǎng)值相乘，最后將計(jì)算區(qū)域?qū)?yīng)網(wǎng)格點(diǎn)相加即可得到最終成像結(jié)果.

本文開(kāi)發(fā)的電磁波逆時(shí)偏移成像算法如圖1所示，對(duì)于合成數(shù)據(jù)，首先給定真實(shí)模型，通過(guò)FDTD模擬電磁波在模型中的傳播過(guò)程，同時(shí)記錄接收波場(chǎng).逆時(shí)偏移成像的第一步是通過(guò)FDTD進(jìn)行激勵(lì)源波場(chǎng)的正向傳播模擬，使用雷克子波作為激勵(lì)源函數(shù)，電磁波最大傳播時(shí)間序列為nmax，為了節(jié)約存儲(chǔ)空間，將每個(gè)時(shí)刻模型空間的外圍波場(chǎng)存儲(chǔ)；第二步是電磁波波場(chǎng)的逆向傳播模擬，該步驟包括兩個(gè)方面的反向傳播過(guò)程：(一)利用包圍模型空間外圍全時(shí)序波場(chǎng)恢復(fù)電磁波場(chǎng)的正向傳播過(guò)程，這樣做的目的其實(shí)就是以時(shí)間換空間；(二)利用接收波場(chǎng)進(jìn)行電磁波場(chǎng)的逆向傳播模擬；第三步是成像條件的應(yīng)用，因?yàn)榛ハ嚓P(guān)條件容易實(shí)現(xiàn)且無(wú)穩(wěn)定性問(wèn)題，本文選取互相關(guān)作為成像條件.波場(chǎng)逆向傳播的過(guò)程中，利用成像條件獲得每個(gè)時(shí)刻的偏移圖像，當(dāng)波場(chǎng)傳播到時(shí)間零點(diǎn)，將各個(gè)時(shí)刻的圖像疊加即可獲得單源逆時(shí)偏移成像結(jié)果.此外，通過(guò)去噪技術(shù)去除低頻影響，再將所有單源偏移結(jié)果疊加得到最終的逆時(shí)偏移成像結(jié)果.

圖1 電磁波逆時(shí)偏移成像算法流程Fig.1 RTM algorithm flow for electromagnetic wave

3 數(shù)值模擬實(shí)驗(yàn)

圖2 相對(duì)介電常數(shù)模型Fig.2 Relative dielectric constant model

文中構(gòu)建了一個(gè)復(fù)雜地質(zhì)結(jié)構(gòu)模型如圖2所示，X方向長(zhǎng)度為1.2 m，Z方向長(zhǎng)度為0.7 m.圖中藍(lán)色區(qū)域代表空氣，暗紅色區(qū)域代表土壤層，其相對(duì)節(jié)電常數(shù)為6，電導(dǎo)率為5 mS·m-1，平均厚度約為0.12 m，該土壤層的底界面起伏變化；黃色區(qū)域代表干沙層，其相對(duì)介電常數(shù)為4，電導(dǎo)率為1 mS·m-1，平均厚度約為0.48 m.由于土壤層與干沙層之間存在大傾角起伏界面，故而該模型中存在變化劇烈的橫向速度.干沙中埋有2個(gè)異常體，其中藍(lán)色方形區(qū)域代表一個(gè)高速異常體，其相對(duì)介電常數(shù)為1；圓心位于(0.8, 0.4)m處圓形暗紅色區(qū)域代表一個(gè)低速異常體，其相對(duì)介電常數(shù)為6，電導(dǎo)率為4 mS·m-1.

3.1 雷達(dá)合成數(shù)據(jù)

典型的雷達(dá)反射測(cè)量中，電場(chǎng)與測(cè)線往往是垂直的，因此我們應(yīng)用二維TM模FDTD求解Maxwell方程組，同時(shí)采用UPML吸收邊界條件來(lái)減少邊界處反射.為確保穩(wěn)定性及減少數(shù)值頻散，設(shè)置網(wǎng)格大小為0.004 m，發(fā)射源為主頻1600 MHz的雷克子波，時(shí)間取樣間隔0.008 ns，采樣點(diǎn)為1024.首先合成了共偏移距剖面探地雷達(dá)數(shù)據(jù)，第一道信號(hào)發(fā)射天線的水平距離位于0.05 m，接收天線與發(fā)射天線間距為0.05 m，發(fā)射天線的移動(dòng)步長(zhǎng)為0.012 m，沿地面總該采集了88道雷達(dá)信號(hào)，水平方向上覆蓋了1.05 m的水平距離.圖3a是合成的共偏移距雷達(dá)剖面，橫坐標(biāo)代表道號(hào)，縱坐標(biāo)代表雷達(dá)波的雙程走時(shí)，從圖可知在2.5～4.5 ns之間的雙曲線同向軸是來(lái)自土壤層和沙層界面的反射波，值得注意的是出現(xiàn)了來(lái)自不規(guī)則界面拐點(diǎn)處的若干散射波，這些散射雜波會(huì)對(duì)地下目標(biāo)的識(shí)別造成困擾.水平位置位于第20道，頂點(diǎn)約在6 ns出現(xiàn)的雙曲線同向軸來(lái)自方形高速異常體的反射波，而水平位置位于第60道，頂點(diǎn)約為4.8 ns處出現(xiàn)的雙曲線同向軸則是來(lái)自圓形低速異常體的反射波.合成的共炮點(diǎn)雷達(dá)數(shù)據(jù)如圖3b所示，水平軸代表道號(hào)，縱軸是雷達(dá)波的雙程走時(shí).共炮點(diǎn)雷達(dá)數(shù)據(jù)集與共偏移距雷達(dá)剖面在某種程度上存在一定的相似性，來(lái)自界面及異常的反射波出現(xiàn)的雙程走時(shí)相當(dāng)，此外反射波的曲線形態(tài)也存在一定的相似性.

3.2 基于全波形反演的初始模型建立

逆時(shí)偏移成像算法中，影響最終成像效果的一個(gè)最重要因素是初始模型.為了獲得一個(gè)較高精確度的初始模型，本文采用時(shí)間域全波形反演技術(shù)，反演地下介質(zhì)的相對(duì)介電常數(shù)分布，從而為逆時(shí)偏移成像提供更為精準(zhǔn)的初始模型.全波形反演技術(shù)試圖利用包括走時(shí)、幅度、相位等全波形信息來(lái)反演地下參數(shù)分布，它比常規(guī)的基于射線的成像方法能提供更加可靠的參數(shù)信息，本文應(yīng)用的時(shí)間域全波形反演技術(shù)細(xì)節(jié)可參考(Meles et al., 2010, 2011)，其主要目的是建立包含觀測(cè)數(shù)據(jù)和模擬數(shù)據(jù)的目標(biāo)函數(shù)(2)式：

(2)

時(shí)間域全波形反演算法中，不管是單參數(shù)分別反演還是介電常數(shù)和電導(dǎo)率聯(lián)合反演，初始介電常數(shù)和電導(dǎo)率模型的選取至關(guān)重要.考慮到如圖2所示的真實(shí)相電常數(shù)模型中存在起伏較大的界面，地下介電常數(shù)模型不易給定，故而這里選取的初始介電常數(shù)模型如圖4a所示，藍(lán)色區(qū)域代表空氣；黃色區(qū)域代表地下介質(zhì)，其相對(duì)介電常數(shù)為4.5.圖4b是時(shí)間域全波形反演25次迭代反演后得到的相對(duì)介電常數(shù)結(jié)果，首先土壤與沙層之間的起伏界面得到了清晰的刻畫(huà)，起伏界面的形狀及位置與真實(shí)模型匹配較好，另外在反演的界面上方相對(duì)介電常數(shù)要高于界面下方的相對(duì)介電常數(shù)；反演圖像中埋于沙層中的兩個(gè)異常體得到了較好的刻畫(huà)，兩個(gè)異常體的水平位置與實(shí)際模型相對(duì)應(yīng).右側(cè)低速異常體的位置與真實(shí)模型一致，而且其內(nèi)部相對(duì)介電常數(shù)值高于周圍背景值；左側(cè)高速異常體的位置與真實(shí)模型一致，其內(nèi)部相對(duì)介電常數(shù)值低于背景值.

圖5給出了真實(shí)模型、初始模型及反演結(jié)果之間的詳細(xì)數(shù)值比較，黑色實(shí)線代表真實(shí)模型的相對(duì)介電常數(shù)值，藍(lán)色虛線代表初始模型相對(duì)介電常數(shù)值，紅色實(shí)線是時(shí)間域全波形反演迭代25次后的反演結(jié)果.圖5a是經(jīng)過(guò)左側(cè)圓形異常體A1-A2線(見(jiàn)圖2)的相對(duì)介電常數(shù)比較，橫軸為相對(duì)介電常數(shù)值，縱軸對(duì)應(yīng)的空間位置.小于0.1m時(shí)，反應(yīng)在真實(shí)模型中是相對(duì)介電常數(shù)為1空氣，可以發(fā)現(xiàn)三條曲線完全重合.0.1～0.3m之間，給出的初始模型與真實(shí)模型之間存在1.5的差異，而全波形反演結(jié)果基本介于初始模型與真實(shí)模型之間.在0.3m處是界面突變點(diǎn)，反演結(jié)果中在同樣位置出現(xiàn)了拐點(diǎn).0.5m附近是異常體所在的位置，反演結(jié)果的數(shù)值大小及拐點(diǎn)位置與真實(shí)模型基本一致.圖5b是經(jīng)過(guò)右側(cè)圓形異常體B1-B2線(見(jiàn)圖2)的相對(duì)介電常數(shù)比較，從中可以得到與左圖相似的結(jié)論：全波反演的結(jié)果較為準(zhǔn)確的刻畫(huà)了地下相對(duì)介電常數(shù)分布情況，可以為逆時(shí)偏移提供更為精準(zhǔn)的初始模型.

圖4 (a)初始介電常數(shù)模型； (b)全波形反演25次迭代后的介電常數(shù)反演結(jié)果Fig.4 (a) The initial permittivity model for the FWI inversion； (b) the inversion permittivity results after 25 iterations

圖5 真實(shí)模型，初始模型及時(shí)間域全波反演結(jié)果在不同位置的介電常數(shù)比較(a) A1-A2線(見(jiàn)圖2)介電常數(shù)比較; (b) B1-B2線(見(jiàn)圖2)介電常數(shù)比較.Fig.5 The permittivity comparison of true model, initial model and FWI inversion results at different position(a) Permittivity comparison of A1-A2 line (the left black dash line indicated in Fig.2); (b) Permittivity comparison of B1-B2 line (the right black dash line indicated).

3.3 合成數(shù)據(jù)逆時(shí)偏移成像分析

為了評(píng)價(jià)逆時(shí)偏移算法應(yīng)用到探地雷達(dá)數(shù)據(jù)處理中的效果，選取Kirchhoff偏移算法與之比較.對(duì)于逆時(shí)偏移成像算法，我們可以將其處理共偏移距數(shù)據(jù)，亦可以處理共炮點(diǎn)數(shù)據(jù).

對(duì)于如圖3a所示的共偏移距雷達(dá)數(shù)據(jù)，基于真速度模型，用Kirchhoff偏移成像算法獲得的偏移成像結(jié)果圖6a所示.從圖可知，來(lái)自不規(guī)則界面的反射波得到一定程度的抑制，可以大致識(shí)別出土壤與沙層的界面，然而土壤與沙層界面的散射波沒(méi)有完全收斂，仍然存在較嚴(yán)重的雜波干擾.偏移剖面中，來(lái)自兩個(gè)異常體的反射波能量基本上得到了壓制.左側(cè)高速異常體水平位置與真實(shí)模型一致，但是垂直位置與實(shí)際模型存在一定的偏差；右側(cè)低速異常體的位置與實(shí)際模型基本一致.此外需要注意的是成像結(jié)果中出現(xiàn)了一些干擾雜波，這些雜波的存在將使得雷達(dá)數(shù)據(jù)的解釋工作更加復(fù)雜.因此，對(duì)于存在橫向速度變化較大的模型，即便使用真實(shí)速度模型進(jìn)行Kirchhoff偏移，也得不到良好的偏移結(jié)果.

圖6 不同偏移方法應(yīng)用不同初始模型得到的成像結(jié)果(a) 基于共偏移距數(shù)據(jù)的Kirchhoff偏移結(jié)果(真介電常數(shù)模型作為初始模型); (b)基于共偏移距數(shù)據(jù)的RTM結(jié)果(真介電常數(shù)模型作為初始模型); (c) 基于共偏移距數(shù)據(jù)的RTM結(jié)果(地下相對(duì)介電常數(shù)取4作為初始模型); (d) 基于共偏移距數(shù)據(jù)的RTM結(jié)果(地下相對(duì)介電常數(shù)取6作為初始模型); (e) 基于共偏移距數(shù)據(jù)的RTM結(jié)果(FWI反演結(jié)果作為初始模型); (f) 基于共炮點(diǎn)雷達(dá)數(shù)據(jù)集的RTM結(jié)果(真介電常數(shù)模型作為初始模型). 所有圖像都進(jìn)行了振幅歸一化.Fig.6 Imaging results from different migration methods with different initial model(a) The Kirchhoff migration result based on common offset radar data (true dielectric constant model as the initial model); (b) RTM result based on common offset radar data (true dielectric constant model as the initial model); (c) RTM result based on common offset radar data (dielectric constant of 4 in the underground as the initial model); (d) RTM result based on common offset radar data (dielectric constant of 6 in the underground as the initial model); (e) RTM result based on common offset radar data (FWI result as the initial model); (f) RTM result based on common shot radar data (true dielectric constant model as the initial model). All figures are normalized by their maximum amplitude.

同樣，對(duì)于如圖3a所示的共偏移距雷達(dá)數(shù)據(jù)，將真介電常數(shù)模型作為初始模型參數(shù)，忽略電導(dǎo)率的影響，用逆時(shí)偏移成像獲得的成像結(jié)果如圖6b所示，可以發(fā)現(xiàn)來(lái)自起伏界面的散射波已被完全壓制住了，土壤與沙層間的界面可以被清晰的識(shí)別出來(lái)，成像后的界面形狀及位置與真實(shí)模型完全一致.來(lái)自兩個(gè)異常體的反射波聚焦較好，成像結(jié)果顯示兩個(gè)異常體的形態(tài)大小，空間位置與真實(shí)模型完全一致，此外偏移剖面中存在的雜波干擾較小.總體來(lái)說(shuō)，逆時(shí)偏移成像效果明顯好于Kirchhoff偏移結(jié)果.

逆時(shí)偏移成像算法中，初始模型的選擇至關(guān)重要，初始模型的建立可以結(jié)合已知的先驗(yàn)信息、CMP速度分析、或者反射層析成像結(jié)果.為了研究初始模型對(duì)成像精度的影響，文中分別選取地下相對(duì)介電常數(shù)為4的均勻半空間，地下相對(duì)介電常數(shù)為6的均勻半空間以及選取全波形反演的結(jié)果作為初始速度模型.

對(duì)于共偏移距探地雷達(dá)數(shù)據(jù)，當(dāng)選取地下相對(duì)介電常數(shù)為4的均勻半空間作為初始介電常數(shù)模型時(shí)，采用逆時(shí)偏移技術(shù)獲得的成像結(jié)果如圖6c所示，從圖可知來(lái)自起伏界面的繞射波沒(méi)有完全歸位，雜波干擾嚴(yán)重，來(lái)自地下異常體的散射波能量聚焦效果較差.當(dāng)選取地下相對(duì)介電常數(shù)為6的均勻半空間作為初始介電常數(shù)模型時(shí)，逆時(shí)偏移成像結(jié)果如圖6d所示，來(lái)自起伏界面的反射波能量聚焦較好，能清晰的識(shí)別起伏界面的位置及形態(tài)，但是來(lái)自地下異常體的散射波能量沒(méi)有完全聚焦.當(dāng)選取全波形反演的結(jié)果作為初始介電常數(shù)模型時(shí)，逆時(shí)偏移成像的結(jié)果如圖6e所示，從圖可知，來(lái)自界面的反射波能量聚焦較好，能夠清晰的識(shí)別起伏界面的形態(tài)和位置；來(lái)自地下的兩個(gè)異常體的散射波能量聚焦較好，其空間位置與模型基本一致.由此可知，當(dāng)初始模型越接近于真實(shí)模型時(shí)，逆時(shí)偏移成像結(jié)果越準(zhǔn)確；全波形反演能夠?yàn)槟鏁r(shí)偏移成像提供更為準(zhǔn)確的初始模型，將全波形反演的結(jié)果作為逆時(shí)偏移成像的初始模型是行之有效的.

對(duì)于如圖3b所示的共炮點(diǎn)探地雷達(dá)數(shù)據(jù)，采用逆時(shí)偏移成像算法對(duì)其進(jìn)行偏移成像，其結(jié)果如圖6d所示，圖中來(lái)自起伏界面的反射波能量得到聚焦，起伏界面的空間位置及形態(tài)與實(shí)際模型基本一致；來(lái)自地下異常體的散射波能量聚焦較好，其空間位置與模型一致.偏移剖面中，在地面處出現(xiàn)了若干個(gè)和源位置一致的強(qiáng)能量點(diǎn)；其反射波子波略寬于如圖6b中的反射波子波寬度.

4 物理實(shí)驗(yàn)及數(shù)據(jù)分析

4.1 物理實(shí)驗(yàn)?zāi)Ｐ图皵?shù)據(jù)采集

本文在實(shí)驗(yàn)室沙槽進(jìn)行了物理模擬實(shí)驗(yàn)，如圖7a所示，沙槽長(zhǎng)8 m，寬1.6 m，深0.8 m，沙槽四壁為混泥土，一個(gè)PVC盒以45°傾斜角埋于沙中，其埋深約0.24 m.圖7c展示的是PVC盒埋于干沙中的示意圖，PVC盒厚度為0.04 m，其中間為空氣，頂端位于(0.75, 0.24)m處，干沙相對(duì)介電常數(shù)是4，PVC的相對(duì)介電常數(shù)是3.1.圖7b展示了本實(shí)驗(yàn)中采用的 MALA ProEx雷達(dá)控制單元及兩套中心頻率為1600 MHz的高頻天線.當(dāng)采集共炮點(diǎn)雷達(dá)數(shù)據(jù)時(shí)，雷達(dá)主機(jī)外接兩個(gè)高頻模塊，同時(shí)將兩套高頻天線分別接于高頻模塊；當(dāng)采集共偏移距雷達(dá)數(shù)據(jù)時(shí)，雷達(dá)主機(jī)外接一個(gè)高頻模塊，同時(shí)只需要一個(gè)高頻天線.

首先用1600 MHz高頻天線進(jìn)行了共偏移距測(cè)量，測(cè)量過(guò)程中天線緊貼干沙表面，圖8a所示是實(shí)測(cè)的共偏移距探地雷達(dá)剖面，橫軸為道號(hào)，縱軸為雷達(dá)波雙程走時(shí).圖中0～1 ns之間是天線耦合直達(dá)波，頂點(diǎn)位于3.2 ns的雙曲線是來(lái)自PVC盒子頂端的散射波.同時(shí)可以發(fā)現(xiàn)剖面中出現(xiàn)了一些微弱的散射雜波，雜波是由于干沙局部不均勻性造成的.

其次，采集了共炮點(diǎn)雷達(dá)數(shù)據(jù)集，第一個(gè)發(fā)射天線位于起始位置，第一個(gè)接收位置距離發(fā)射天線的間距為0.1 m，每一個(gè)發(fā)射位置對(duì)應(yīng)21個(gè)接收位置，即單炮雷達(dá)數(shù)據(jù)包含了21道信號(hào)，接收器間的間隔是0.01 m，發(fā)射天線的移動(dòng)間距為0.05 m.圖8b是實(shí)測(cè)的共炮點(diǎn)探地雷達(dá)數(shù)據(jù)集，它包含36個(gè)共炮點(diǎn)探地雷達(dá)剖面數(shù)據(jù).這些數(shù)據(jù)將用于時(shí)間域全波形反演，以此來(lái)反演地下介質(zhì)分布，從而為疊前逆時(shí)偏移成像提供初始速度模型.

4.2 雷達(dá)實(shí)測(cè)數(shù)據(jù)的全波形反演

實(shí)測(cè)數(shù)據(jù)必須經(jīng)過(guò)必要的數(shù)據(jù)預(yù)處理，包括對(duì)共炮點(diǎn)探地雷達(dá)數(shù)據(jù)進(jìn)行直流偏移和帶通濾波，同時(shí)考慮到本文所用的時(shí)間域全波形反演算法是二維的，而實(shí)測(cè)雷達(dá)數(shù)據(jù)是三維的，故采用(Bleistein, 1986)提出的三維轉(zhuǎn)二維技術(shù)，將三維實(shí)測(cè)數(shù)據(jù)轉(zhuǎn)為二維情況下的數(shù)據(jù)，從而消除它們?cè)诜群拖辔簧系牟町?由于在沙中僅埋有一個(gè)PVC盒子，故而在時(shí)間域全波形反演中，采用相對(duì)介電常數(shù)為4均勻模型作為的初始相對(duì)介電常數(shù)模型(圖9a).以實(shí)測(cè)雷達(dá)數(shù)據(jù)和計(jì)算雷達(dá)數(shù)據(jù)的均方根作為反演迭代次數(shù)的目標(biāo)函數(shù)，如圖10所示，經(jīng)過(guò)15次迭代后，目標(biāo)函數(shù)的誤差改變小于1%.圖9b是迭代15次后的相對(duì)介電常數(shù)反演結(jié)果，圖中可清晰辨別出傾斜掩埋的PVC盒子，其形狀和位置與實(shí)際模型基本一致，傾角也非常接近45°，此外反演結(jié)果中小散射點(diǎn)清晰可見(jiàn).

圖7 物理模擬實(shí)驗(yàn)(a) PVC盒子傾斜掩埋在沙中; (b) 實(shí)驗(yàn)用到的儀器：Mala ProEx控制單元及兩套1600 MHz的高頻天線;(c) PVC盒子埋于干沙中的示意圖.Fig.7 Laboratory experiment setup(a) The PVC box tilted buried in sand; (b) Instruments used to carry out radar measurements consisting of Mala ProEx control unit and two sets of 1.6 GHz antenna; (c) A sketch of the PVC box embedded in sand.

圖8 物理實(shí)驗(yàn)實(shí)測(cè)數(shù)據(jù)(a) 實(shí)測(cè)的共偏移距剖面; (b) 實(shí)測(cè)的包含36條共炮點(diǎn)雷達(dá)剖面的雷達(dá)數(shù)據(jù)集，所有圖像均進(jìn)行了歸一化.Fig.8 Physical experiments radar data(a) The measured common-offset profile; (b) The measured radar gather consisting of 36 common-shot radar profiles, all figures are normalized by their maximum amplitude.

圖9 FWI的初始介電常數(shù)模型(a) 相對(duì)介電常數(shù)定為4; (b) 15次反演迭代后的介電常數(shù)結(jié)果.Fig.9 The initial permittivity model for the FWI(a) The dielectric constant is set to 4; (b) The inversion permittivity results after 15 iterations.

圖10 以實(shí)測(cè)雷達(dá)數(shù)據(jù)和計(jì)算雷達(dá)數(shù)據(jù)間均方根值作為迭代次數(shù)的函數(shù)(均方根值經(jīng)過(guò)以做歸一化，而且15次迭代后均方根改變值小于1%)Fig.10 RMS values as a function of the iteration number for the observed radar data and the calculated data (The RMS is normalized to its maximum and after 15 iterations; the RMS misfit changes less than 1%)

圖11 (a)實(shí)測(cè)共偏移距剖面Kirchhoff偏移結(jié)果; (b)實(shí)測(cè)共炮點(diǎn)雷達(dá)數(shù)據(jù)得到的逆時(shí)偏移成像結(jié)果，兩圖均已歸一化處理Fig.11 (a) Kirchhoff migration result from measured common-offset profile, (b) RTM result from measured common-shot radar gather; all figures are normalized by their maximum amplitude

4.3 探地雷達(dá)實(shí)測(cè)數(shù)據(jù)的逆時(shí)偏移結(jié)果

通過(guò)時(shí)間域全波形反演獲得了地下介電常數(shù)分布情況，基于反演結(jié)果給定的初始介電常數(shù)模型，分別進(jìn)行了Kirchhoff偏移成像和逆時(shí)偏移成像，其結(jié)果見(jiàn)圖11.圖11a是對(duì)實(shí)測(cè)的共偏移距雷達(dá)數(shù)據(jù)進(jìn)行Kirchhoff偏移成像得到的結(jié)果，顯然來(lái)自PVC盒頂?shù)碾p曲線同向軸能量收斂到一個(gè)點(diǎn)，該聚焦點(diǎn)幾乎與PVC盒頂端位置吻合，然而PVC盒身在偏移圖像中并沒(méi)有得到成像，此外，成像結(jié)果中出現(xiàn)了許多散射雜點(diǎn)，這些散射雜點(diǎn)不利于雷達(dá)數(shù)據(jù)的解釋.圖11b是對(duì)實(shí)測(cè)共炮點(diǎn)探地雷達(dá)數(shù)據(jù)進(jìn)行逆時(shí)偏移成像處理獲得的結(jié)果，從圖可知，成像結(jié)果中傾斜PVC盒子整體(白色橢圓框所示)形態(tài)及傾角與實(shí)際模型基本一致，其偏移結(jié)果明顯好于Kirchhoff偏移成像結(jié)果.

5 結(jié)論與討論

本文開(kāi)發(fā)了電磁波二維逆時(shí)偏移成像算法，其中電磁波場(chǎng)的正向傳播和反向傳播均采用時(shí)間域有限差分方法，互相關(guān)成像條件用于獲得成像結(jié)果.考慮到全空間波場(chǎng)存儲(chǔ)內(nèi)存開(kāi)銷較大，文中將成像空間外邊界做為存儲(chǔ)空間，記錄其各個(gè)時(shí)刻的場(chǎng)值，記錄的場(chǎng)值完全能夠恢復(fù)電磁波場(chǎng)的正向傳播過(guò)程.

數(shù)值模擬實(shí)驗(yàn)中，構(gòu)建了橫向速度變化較大的模型，合成了共偏移距和共炮點(diǎn)雷達(dá)數(shù)據(jù).對(duì)于合成的共偏移距雷達(dá)數(shù)據(jù)，同時(shí)采用Kirchhoff偏移成像和逆時(shí)偏移成像算法，對(duì)比成像結(jié)果發(fā)現(xiàn)，逆時(shí)偏移成像剖面中，來(lái)自起伏界面及地下異常體的空間位置和幾何形態(tài)與實(shí)際模型一致，而Kirchhoff偏移剖面中來(lái)自起伏界面和異常體的反射波(散射波)能量聚焦不完全，剖面中仍然存在明顯的雜波，可見(jiàn)逆時(shí)偏移成像算法明顯好于Kirchhoff偏移算法.然而需要注意的是，逆時(shí)偏移成像算法是基于正演計(jì)算，每完成一次偏移需要進(jìn)行兩次以上正演計(jì)算，而正演計(jì)算非常耗時(shí)，故而其計(jì)算效率遠(yuǎn)低于Kirchhoff偏移算法.采用逆時(shí)偏移成像處理共炮點(diǎn)雷達(dá)合成數(shù)據(jù)，得到的成像剖面與該算法處理共偏移距雷達(dá)數(shù)據(jù)獲得的成像剖面基本一致，不同之處在于前者成像剖面中的子波寬度略寬于后者.

影響逆時(shí)偏移成像效果的關(guān)鍵因素之一是包括介電常數(shù)和電導(dǎo)率在內(nèi)的初始模型的選擇，數(shù)值實(shí)驗(yàn)表明更精準(zhǔn)的初始模型能獲得更高精度的偏移成像結(jié)果.為了獲得更精確的初始模型，本文通過(guò)時(shí)間域全波形反演算法對(duì)探地雷達(dá)數(shù)據(jù)進(jìn)行全波形反演，反演后的地下相對(duì)介電常數(shù)模型作為逆時(shí)偏移成像的初始速度模型.數(shù)值模擬實(shí)驗(yàn)及實(shí)際物理模擬實(shí)驗(yàn)表明，通過(guò)全波形反演獲取的信息作為逆時(shí)偏移成像的初始模型是行之有效的，同時(shí)逆時(shí)偏移成像的效果在一定程度上檢驗(yàn)了全波反演結(jié)果的準(zhǔn)確性.

Baysal E, Kosloff D D, Sherwood J W C. 1983. Reverse time migration.Geophysics, 48(11): 1514-1524.

Bleistein N. 1986. Two-and-one-half dimensional in-plane wave propagation.GeophysicalProspecting, 34(5): 686-703.

Chang W F, McMechan G A. 1986. Reverse-time migration of offset vertical seismic profiling data using the excitation-time imaging condition.Geophysics, 51(1): 67-84.

Chang W F, McMechan G A. 1994. 3-D elastic prestack, reverse-time depth migration.Geophysics, 59(4): 597-609.

Chattopadhyay S, McMechan G A. 2008. Imaging conditions for prestack reverse-time migration.Geophysics, 73(3): S81-S89.

Claerbout J F. 1971. Toward a unified theory of reflector mapping.Geophysics, 36(3): 467-481.

Docherty P C. 1991. A brief comparison of some Kirchhoff integral formulas for migration and inversion.Geophysics, 56(8): 1164-1169.

Fang N S, Ying L A. 2009. Stability analysis of FDTD to UPML for time dependent Maxwell equations.ScienceinChinaSeriesA, 52(4): 794-816.

Feng D S, Dai Q W. 2011. GPR numerical simulation of full wave field based on UPML boundary condition of ADI-FDTD.NDT&EInternational, 44(6): 495-504.Feng X, Sato M. 2004. Pre-stack migration applied to GPR for landmine detection.InverseProblems, 20(6): S99-S115.

Fisher E, McMechan G A, Annan A P, et al. 1992. Examples of reverse-time migration of signal-channel ground-penetrating radar profiles.Geophysics, 57(4): 577-586.

Gray S H. 1986. Efficient traveltime calculations for Kirchhoff migration.Geophysics, 51(8): 1685-1688.

Gray S H, May W P. 1994. Kirchhoff migration using eikonal equation traveltimes.Geophysics, 59(5): 810-817.

Guitton A, Kaelin B, Biondi B. 2006. Least-squares attenuation of reverse-time-migration artifacts.Geophysics, 72(1): S19-S23.

Guitton A, Valenciano A, Bevc D, et al. 2007. Smoothing imaging condition for shot-profile migration.Geophysics, 72(3): S149-S154. Levin S A. 1984. Principle of reverse-time migration.Geophysics, 49(5): 581-583.

Liu S X, Sato M. 2005. Transient radiation from an unloaded, finite dipole antenna in a borehole: experimental and numerical results.Geophysics, 70(6): K43-K51.

Liu S X, Zeng Z F, Xu B. 2006. FDTD simulation of ground penetrating radar signal in 3-Dimensional dispersive medium.JournalofJilinUniversity(EarthScienceEdition) (in Chinese), 36(1): 123-127.

Liu S X, Zeng Z F, Deng L. 2007. FDTD simulations for ground

penetrating radar in urban applications.JournalofGeophysicsandEngineering, 4(3): 262-267.

Liu S X, Zeng Z F.2007.Numerical simulation for Ground Penetrating Radar wave propagation in the dispersive medium.ChineseJ.Geophys. (in Chinese), 50(1): 320-326.

Liu S X, Wu J J, Zhou J F, et al. 2011. Numerical simulations of borehole radar detection for metal ore.IEEEGeosci.RemoteSens.Lett., 8(2): 308-312.

Liu X X, Yin X Y, Liu B H. 2013. Prestack reverse-time migration of elastic waves in porous media.ProgressinGeophys. (in Chinese), 28(6): 3165-3173.

Loewenthal D, Mufti I R. 1983. Reversed time migration in spatial frequency domain.Geophysics, 48(5): 627-635.

Meles G, Greenhalgh S, van der Kruk J, et al. 2011. Taming the non-linearity problem in GPR full-waveform inversion for high contrast media.JournalofAppliedGeophysics, 73(2): 174-186. Meles G A, van der Kruk J, Greenhalgh S A, et al. 2010. A new vector waveform inversion algorithm for simultaneous updating of conductivity and permittivity parameters from combination crosshole/borehole-to-surface GPR Data.IEEETransactionsonGeoscienceandRemoteSensing, 48(9): 3391-3407.

Mufti I R, Pita J A, Huntley R W. 1996. Finite-difference depth migration of exploration-scale 3-D seismic data.Geophysics, 61(3): 776-794.

Sava P, Fomel S. 2006. Time-shift imaging condition in seismic migration.Geophysics, 71(6): S209-S217.

Sun R, McMechanz G A. 2001. Scalar reverse-time depth migration of prestack elastic seismic data.Geophysics, 66(5): 1519-1527. Taflove A, Hagness S. 2005. Computational Electrodynamics - The Finite-Difference Time-Domain Method, 3rd Ed. London: Artech House.Vidale J E. 1988. Finite-difference calculation of travel times.Bull.Seismol.Soc.Am., 78(6): 2062-2076.

Wang B L, Gao J H, Chen W C, et al. 2012. Efficient boundary storage strategies for seismic reverse time migration.ChineseJ.Geophys. (in Chinese), 55(7): 2412-2421, doi: 10.6038/j.issn.0001-5733.2012.07.025.

Yoon K, Marfurt K J. 2006. Reverse-time migration using the Poynting vector.ExplorationGeophysics, 37(1): 102-107.

Yoon K, Shin C, Suh S, et al. 2003. 3D reverse-time migration using the acoustic wave equation: An experience with the SEG/EAGE data set.TheLeadingEdge, 22(1): 38-41.

Zhou H, Sato M, Liu H J. 2005. Migration velocity analysis and prestack migration of common-transmitter GPR data.IEEETransactionsonGeoscienceandRemoteSensing, 43(1): 86-91.

附中文參考文獻(xiàn)

劉四新, 曾昭發(fā), 徐波. 2006. 三維頻散介質(zhì)中地質(zhì)雷達(dá)信號(hào)的FDTD數(shù)值模擬. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版), 36(1): 123-127. 劉四新, 曾昭發(fā). 2007. 頻散介質(zhì)中地質(zhì)雷達(dá)波傳播的數(shù)值模擬. 地球物理學(xué)報(bào), 50(1): 320-326.

劉欣欣, 印興耀, 劉百紅. 2013. 孔隙介質(zhì)彈性波疊前逆時(shí)偏移方法. 地球物理學(xué)進(jìn)展, 28(6): 3165-3173.

王保利, 高靜懷, 陳文超等. 2012. 地震疊前逆時(shí)偏移的有效邊界存儲(chǔ)策略. 地球物理學(xué)報(bào), 55(7): 2412-2421, doi: 10.6038/j.issn.0001-5733.2012.07.025.

(本文編輯 汪海英)

Reverse time migration applied to GPR data based on full wave inversion

LEI Lin-Lin1, LIU Si-Xin1, FU Lei1*, WU Jun-Jun2

1CollegeofGeo-explorationScienceandTechnology,JilinUniversity,Changchun130026,China2BGPInc.,ChinaNationalPetroleumCorporation,ZhuozhouHebei072751,China

Reverse-time migration (RTM) is used for subsurface imaging to handle complex velocity models including steeply dipping interfaces and dramatic lateral variations, and promises better imaging results compared to traditional migration methods such as Kirchhoff migration algorithm. RTM has been increasingly applied to seismic surveys for hydrocarbon resource explorations. Based on the similarity of kinematics and dynamics between electromagnetic waves and elastic waves, we develop a pre-stack RTM method and apply it to processing ground penetrating radar (GPR) data.The finite-difference time domain (FDTD) numerical method is used to simulate the electromagnetic wave propagation including forward and backward extrapolation. The cross-correlation imaging condition is used to obtain the final image. In order to provide a velocity model with relatively higher accuracy as the initial velocity model for RTM, we apply a full waveform inversion (FWI) in the time domain to estimate the subsurface velocity structure based on reflection radar data. For testing the effectiveness of the algorithm, we have constructed a complex geological model; and synthesized common-offset radar data and common-shot profile (CSP) radar reflection data. All data are migrated with the traditional Kirchhoff migration method and pre-stack RTM method separately. The migration results from pre-stack RTM show better coincidence with the true model. Furthermore, we have performed a physical experiment in a sandbox where a polyvinyl chloride (PVC) box is buried in the sand. The obtained common-offset radar data and common-shot radar data are migrated by using the Kirchhoff migration method and per-stack RTM algorithm separately. The per-stack RTM result shows that the RTM algorithm can obtain better imaging results.In the numerical experiments, velocity variation exists in the geological model due to the irregular interface between soil and sand. The imaging result from pre-stack RTM is better than that of the traditional Kirchhoff migration algorithm because pre-stack RTM can handle the model with horizontal velocity variation. In the RTM imaging, the shape and position of the interface and their anomalies match well with the true model, while the imaging of Kirchhoff migration cannot do so. As for the physical experiment, the Kirchhoff migration algorithm and pre-stack RTM have been applied to process the measured radar data. The imaging results from RTM show its advantage compared with the Kirchhoff migration method.

Reverse time migration; GPR; Full waveform inversion; Physical model experiment

雷林林, 劉四新, 傅磊等. 2015. 基于全波形反演的探地雷達(dá)數(shù)據(jù)逆時(shí)偏移成像.地球物理學(xué)報(bào)，58(9):3346-3355,

10.6038/cjg20150927.

Lei L L, Liu S X, Fu L, et al. 2015. Reverse time migration applied to GPR data based on full wave inversion.ChineseJ.Geophys. (in Chinese),58(9):3346-3355,doi:10.6038/cjg20150927.

10.6038/cjg20150927

P631

2014-03-14，2015-08-02收修定稿

國(guó)家自然科學(xué)基金項(xiàng)目(40874073,41074076),國(guó)家高技術(shù)研究發(fā)展計(jì)劃項(xiàng)目(2013AA064603)及中石油東方地球物理公司中青年科技創(chuàng)新基金項(xiàng)目(11-06-2013)共同資助.

雷林林，女，1987年生，博士生，主要從事地球物理數(shù)據(jù)處理與正反演研究.E-mail:leilljlu@126.com

*通訊作者 傅磊，男，1985年生，講師.E-mail:fuleijlu@163.com