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

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

1逆時偏移基本原理及在地質(zhì)解釋中的應(yīng)用問題偏移是一種地震圖像技術(shù)，在空間和時間區(qū)域記錄的同向軸回歸到真實的反射界面。傳統(tǒng)的地震偏移成像方法基于波動方程，并將多個散射波用作噪聲處理。然而，對于復(fù)雜地質(zhì)結(jié)構(gòu)和大傾角界面（如裂縫），為了合理使用多個散射信息非常重要。傳統(tǒng)的波動方程偏移方法如吉爾霍夫積分偏移(。探地雷達電磁波與地震波在運動學(xué)和動力學(xué)存在一定的相似性，此外，采集到的探地雷達數(shù)據(jù)和地震數(shù)據(jù)都是空間位置的時域電壓信號，故而在一定程度上勘探地震中的數(shù)據(jù)處理技術(shù)可以運用到探地雷達數(shù)據(jù)處理中．本文開發(fā)了電磁波逆時偏移算法，電磁波場的正向及逆向傳播是基于時間域有限差分方法，單軸各向異性完全匹配層（UPML)(2逆時偏移基本原理及實現(xiàn)傳統(tǒng)的偏移算法基本是基于單程波方程深度外推，所以在對多分量地震資料處理過程中首先要進行縱、橫波分離，然后用現(xiàn)有的偏移算法分別處理縱波場和橫波場．它的缺陷在于：1）將多分量地震資料當作幾個標量波的簡單疊加，忽略了地震波的矢量特性，從而影響地震多分量資料的處理與解釋精度；2）假如地震波場橫波與縱波沒有完全分離（目前這種現(xiàn)象廣泛存在），而在實際資料的偏移處理過程中又假設(shè)為標量波場，這必然會使偏移結(jié)果產(chǎn)生許多人為假象，影響處理精度．逆時偏移與傳統(tǒng)的偏移方法有很大的不同，逆時偏移可以看成是波場在時間軸上逆向傳播過程，當波場倒退到零時刻，則所有反射波與繞射波的能量都回到最初被反射和繞射的空間位置，應(yīng)用成像條件即可得到最終的偏移成像剖面．逆時偏移的波動理論方法基于全波方程，波場可以沿各個方向傳播，因此該方法不受介質(zhì)傾角限制，適用于速度任意變化的模型．逆時偏移算法的實現(xiàn)概括起來主要包括三個步驟：第一是震源波場正向延拓；第二是接收波場的逆時延拓；第三是應(yīng)用成像條件成像．前兩個步驟的實現(xiàn)可以基于相同的數(shù)值方法，本文應(yīng)用時間域有限差分法（其中IMAG(x）為地下空間x的成像結(jié)果，S(x,t,s本文開發(fā)的電磁波逆時偏移成像算法如圖1所示，對于合成數(shù)據(jù)，首先給定真實模型，通過FDTD模擬電磁波在模型中的傳播過程，同時記錄接收波場．逆時偏移成像的第一步是通過FDTD進行激勵源波場的正向傳播模擬，使用雷克子波作為激勵源函數(shù)，電磁波最大傳播時間序列為n3土壤層結(jié)構(gòu)組成文中構(gòu)建了一個復(fù)雜地質(zhì)結(jié)構(gòu)模型如圖2所示，X方向長度為1.2m,Z方向長度為0.7m．圖中藍色區(qū)域代表空氣，暗紅色區(qū)域代表土壤層，其相對節(jié)電常數(shù)為6，電導(dǎo)率為5mS·m3.1合成的共偏移距雷達數(shù)據(jù)典型的雷達反射測量中，電場與測線往往是垂直的，因此我們應(yīng)用二維TM模FDTD求解Maxwell方程組，同時采用UPML吸收邊界條件來減少邊界處反射．為確保穩(wěn)定性及減少數(shù)值頻散，設(shè)置網(wǎng)格大小為0.004m，發(fā)射源為主頻1600MHz的雷克子波，時間取樣間隔0.008ns，采樣點為1024．首先合成了共偏移距剖面探地雷達數(shù)據(jù)，第一道信號發(fā)射天線的水平距離位于0.05m，接收天線與發(fā)射天線間距為0.05m，發(fā)射天線的移動步長為0.012m，沿地面總該采集了88道雷達信號，水平方向上覆蓋了1.05m的水平距離．圖3a是合成的共偏移距雷達剖面，橫坐標代表道號，縱坐標代表雷達波的雙程走時，從圖可知在2.5～4.5ns之間的雙曲線同向軸是來自土壤層和沙層界面的反射波，值得注意的是出現(xiàn)了來自不規(guī)則界面拐點處的若干散射波，這些散射雜波會對地下目標的識別造成困擾．水平位置位于第20道，頂點約在6ns出現(xiàn)的雙曲線同向軸來自方形高速異常體的反射波，而水平位置位于第60道，頂點約為4.8ns處出現(xiàn)的雙曲線同向軸則是來自圓形低速異常體的反射波．合成的共炮點雷達數(shù)據(jù)如圖3b所示，水平軸代表道號，縱軸是雷達波的雙程走時．共炮點雷達數(shù)據(jù)集與共偏移距雷達剖面在某種程度上存在一定的相似性，來自界面及異常的反射波出現(xiàn)的雙程走時相當，此外反射波的曲線形態(tài)也存在一定的相似性．3.2初始相電常數(shù)模型的建立逆時偏移成像算法中，影響最終成像效果的一個最重要因素是初始模型．為了獲得一個較高精確度的初始模型，本文采用時間域全波形反演技術(shù)，反演地下介質(zhì)的相對介電常數(shù)分布，從而為逆時偏移成像提供更為精準的初始模型．全波形反演技術(shù)試圖利用包括走時、幅度、相位等全波形信息來反演地下參數(shù)分布，它比常規(guī)的基于射線的成像方法能提供更加可靠的參數(shù)信息，本文應(yīng)用的時間域全波形反演技術(shù)細節(jié)可參考（其中：x時間域全波形反演算法中，不管是單參數(shù)分別反演還是介電常數(shù)和電導(dǎo)率聯(lián)合反演，初始介電常數(shù)和電導(dǎo)率模型的選取至關(guān)重要．考慮到如圖2所示的真實相電常數(shù)模型中存在起伏較大的界面，地下介電常數(shù)模型不易給定，故而這里選取的初始介電常數(shù)模型如圖4a所示，藍色區(qū)域代表空氣；黃色區(qū)域代表地下介質(zhì)，其相對介電常數(shù)為4.5．圖4b是時間域全波形反演25次迭代反演后得到的相對介電常數(shù)結(jié)果，首先土壤與沙層之間的起伏界面得到了清晰的刻畫，起伏界面的形狀及位置與真實模型匹配較好，另外在反演的界面上方相對介電常數(shù)要高于界面下方的相對介電常數(shù)；反演圖像中埋于沙層中的兩個異常體得到了較好的刻畫，兩個異常體的水平位置與實際模型相對應(yīng)．右側(cè)低速異常體的位置與真實模型一致，而且其內(nèi)部相對介電常數(shù)值高于周圍背景值；左側(cè)高速異常體的位置與真實模型一致，其內(nèi)部相對介電常數(shù)值低于背景值．圖5給出了真實模型、初始模型及反演結(jié)果之間的詳細數(shù)值比較，黑色實線代表真實模型的相對介電常數(shù)值，藍色虛線代表初始模型相對介電常數(shù)值，紅色實線是時間域全波形反演迭代25次后的反演結(jié)果．圖5a是經(jīng)過左側(cè)圓形異常體A1-A2線（見圖2）的相對介電常數(shù)比較，橫軸為相對介電常數(shù)值，縱軸對應(yīng)的空間位置．小于0.1m時，反應(yīng)在真實模型中是相對介電常數(shù)為1空氣，可以發(fā)現(xiàn)三條曲線完全重合．0.1～0.3m之間，給出的初始模型與真實模型之間存在1.5的差異，而全波形反演結(jié)果基本介于初始模型與真實模型之間．在0.3m處是界面突變點，反演結(jié)果中在同樣位置出現(xiàn)了拐點．0.5m附近是異常體所在的位置，反演結(jié)果的數(shù)值大小及拐點位置與真實模型基本一致．圖5b是經(jīng)過右側(cè)圓形異常體B1-B2線（見圖2）的相對介電常數(shù)比較，從中可以得到與左圖相似的結(jié)論：全波反演的結(jié)果較為準確的刻畫了地下相對介電常數(shù)分布情況，可以為逆時偏移提供更為精準的初始模型．3.3理論模型預(yù)處理為了評價逆時偏移算法應(yīng)用到探地雷達數(shù)據(jù)處理中的效果，選取Kirchhoff偏移算法與之比較．對于逆時偏移成像算法，我們可以將其處理共偏移距數(shù)據(jù)，亦可以處理共炮點數(shù)據(jù)．對于如圖3a所示的共偏移距雷達數(shù)據(jù)，基于真速度模型，用Kirchhoff偏移成像算法獲得的偏移成像結(jié)果圖6a所示．從圖可知，來自不規(guī)則界面的反射波得到一定程度的抑制，可以大致識別出土壤與沙層的界面，然而土壤與沙層界面的散射波沒有完全收斂，仍然存在較嚴重的雜波干擾．偏移剖面中，來自兩個異常體的反射波能量基本上得到了壓制．左側(cè)高速異常體水平位置與真實模型一致，但是垂直位置與實際模型存在一定的偏差；右側(cè)低速異常體的位置與實際模型基本一致．此外需要注意的是成像結(jié)果中出現(xiàn)了一些干擾雜波，這些雜波的存在將使得雷達數(shù)據(jù)的解釋工作更加復(fù)雜．因此，對于存在橫向速度變化較大的模型，即便使用真實速度模型進行Kirchhoff偏移，也得不到良好的偏移結(jié)果．同樣，對于如圖3a所示的共偏移距雷達數(shù)據(jù)，將真介電常數(shù)模型作為初始模型參數(shù)，忽略電導(dǎo)率的影響，用逆時偏移成像獲得的成像結(jié)果如圖6b所示，可以發(fā)現(xiàn)來自起伏界面的散射波已被完全壓制住了，土壤與沙層間的界面可以被清晰的識別出來，成像后的界面形狀及位置與真實模型完全一致．來自兩個異常體的反射波聚焦較好，成像結(jié)果顯示兩個異常體的形態(tài)大小，空間位置與真實模型完全一致，此外偏移剖面中存在的雜波干擾較?。傮w來說，逆時偏移成像效果明顯好于Kirchhoff偏移結(jié)果．逆時偏移成像算法中，初始模型的選擇至關(guān)重要，初始模型的建立可以結(jié)合已知的先驗信息、CMP速度分析、或者反射層析成像結(jié)果．為了研究初始模型對成像精度的影響，文中分別選取地下相對介電常數(shù)為4的均勻半空間，地下相對介電常數(shù)為6的均勻半空間以及選取全波形反演的結(jié)果作為初始速度模型．對于共偏移距探地雷達數(shù)據(jù)，當選取地下相對介電常數(shù)為4的均勻半空間作為初始介電常數(shù)模型時，采用逆時偏移技術(shù)獲得的成像結(jié)果如圖6c所示，從圖可知來自起伏界面的繞射波沒有完全歸位，雜波干擾嚴重，來自地下異常體的散射波能量聚焦效果較差．當選取地下相對介電常數(shù)為6的均勻半空間作為初始介電常數(shù)模型時，逆時偏移成像結(jié)果如圖6d所示，來自起伏界面的反射波能量聚焦較好，能清晰的識別起伏界面的位置及形態(tài)，但是來自地下異常體的散射波能量沒有完全聚焦．當選取全波形反演的結(jié)果作為初始介電常數(shù)模型時，逆時偏移成像的結(jié)果如圖6e所示，從圖可知，來自界面的反射波能量聚焦較好，能夠清晰的識別起伏界面的形態(tài)和位置；來自地下的兩個異常體的散射波能量聚焦較好，其空間位置與模型基本一致．由此可知，當初始模型越接近于真實模型時，逆時偏移成像結(jié)果越準確；全波形反演能夠為逆時偏移成像提供更為準確的初始模型，將全波形反演的結(jié)果作為逆時偏移成像的初始模型是行之有效的．對于如圖3b所示的共炮點探地雷達數(shù)據(jù)，采用逆時偏移成像算法對其進行偏移成像，其結(jié)果如圖6d所示，圖中來自起伏界面的反射波能量得到聚焦，起伏界面的空間位置及形態(tài)與實際模型基本一致；來自地下異常體的散射波能量聚焦較好，其空間位置與模型一致．偏移剖面中，在地面處出現(xiàn)了若干個和源位置一致的強能量點；其反射波子波略寬于如圖6b中的反射波子波寬度．4物理實驗和數(shù)據(jù)分析4.1雷達高頻天線本文在實驗室沙槽進行了物理模擬實驗，如圖7a所示，沙槽長8m，寬1.6m，深0.8m，沙槽四壁為混泥土，一個PVC盒以45°傾斜角埋于沙中，其埋深約0.24m．圖7c展示的是PVC盒埋于干沙中的示意圖，PVC盒厚度為0.04m，其中間為空氣，頂端位于（0.75,0.24)m處，干沙相對介電常數(shù)是4,PVC的相對介電常數(shù)是3.1．圖7b展示了本實驗中采用的MALAProEx雷達控制單元及兩套中心頻率為1600MHz的高頻天線．當采集共炮點雷達數(shù)據(jù)時，雷達主機外接兩個高頻模塊，同時將兩套高頻天線分別接于高頻模塊；當采集共偏移距雷達數(shù)據(jù)時，雷達主機外接一個高頻模塊，同時只需要一個高頻天線．首先用1600MHz高頻天線進行了共偏移距測量，測量過程中天線緊貼干沙表面，圖8a所示是實測的共偏移距探地雷達剖面，橫軸為道號，縱軸為雷達波雙程走時．圖中0～1ns之間是天線耦合直達波，頂點位于3.2ns的雙曲線是來自PVC盒子頂端的散射波．同時可以發(fā)現(xiàn)剖面中出現(xiàn)了一些微弱的散射雜波，雜波是由于干沙局部不均勻性造成的．其次，采集了共炮點雷達數(shù)據(jù)集，第一個發(fā)射天線位于起始位置，第一個接收位置距離發(fā)射天線的間距為0.1m，每一個發(fā)射位置對應(yīng)21個接收位置，即單炮雷達數(shù)據(jù)包含了21道信號，接收器間的間隔是0.01m，發(fā)射天線的移動間距為0.05m．圖8b是實測的共炮點探地雷達數(shù)據(jù)集，它包含36個共炮點探地雷達剖面數(shù)據(jù)．這些數(shù)據(jù)將用于時間域全波形反演，以此來反演地下介質(zhì)分布，從而為疊前逆時偏移成像提供初始速度模型．4.2雷達實測數(shù)據(jù)的全波形反演實測數(shù)據(jù)必須經(jīng)過必要的數(shù)據(jù)預(yù)處理，包括對共炮點探地雷達數(shù)據(jù)進行直流偏移和帶通濾波，同時考慮到本文所用的時間域全波形反演算法是二維的，而實測雷達數(shù)據(jù)是三維的，故采用（4.3實測數(shù)據(jù)處理結(jié)果通過時間域全波形反演獲得了地下介電常數(shù)分布情況，基于反演結(jié)果給定的初始介電常數(shù)模型，分別進行了Kirchhoff偏移成像和逆時偏移成像，其結(jié)果見圖11．圖11a是對實測的共偏移距雷達數(shù)據(jù)進行Kirchhoff偏移成像得到的結(jié)果，顯然來自PVC盒頂?shù)碾p曲線同向軸能量收斂到一個點，該聚焦點幾乎與PVC盒頂端位置吻合，然而PVC盒身在偏移圖像中并沒有得到成像，此外，成像結(jié)果中出現(xiàn)了許多散射雜點，這些散射雜點不利于雷達數(shù)據(jù)的解釋．圖11b是對實測共炮點探地雷達數(shù)據(jù)進行逆時偏移成像處理獲得的結(jié)果，從圖可知，成像結(jié)果中傾斜PVC盒子整體（白色橢圓框所示）形態(tài)及傾角與實際模型基本一致，其偏移結(jié)果明顯好于Kirchhoff偏移成像結(jié)果．5逆時偏移與kirchihoff偏移算法的對比本文開發(fā)了電磁波二維逆時偏移成像算法，其中電磁波場的正向傳播和反向傳播均采用時間域有限差分方法，互相關(guān)成像條件用于獲得成像結(jié)果．考慮到全空間波場存儲內(nèi)存開銷較大，文中將成像空間外邊界做為存儲空間，記錄其各個時刻的場值，記錄的場值完全能夠恢復(fù)電磁波場的正向傳播過程．數(shù)值模擬實驗中，構(gòu)建了橫向速度變化較大的模型，合成了共偏移距和共炮點雷達數(shù)據(jù)．對于合成的共偏移距雷達數(shù)據(jù)，同時采用Kirchhoff偏移成像和逆時偏移成像算法，對比成像結(jié)果發(fā)現(xiàn)，逆時偏移成像剖面中，來自起伏界面及地下異常體的空間位置和幾何形態(tài)與實際模型一致，而Kirchhoff偏移剖面中來自起伏界面和異常體的反射波（散射波）能量聚焦不完全，剖面中仍然存在明顯的雜波，可見逆時偏移成像算法明顯好于Kirchhoff偏移算法．然而需要注意的是，逆時偏移成像算法是基于正演計算，每完成一次偏移需要進行兩次以上正演計算，而正演計算非常耗時，故而其計算效率遠低于Kirchhoff偏移算法．采用逆時偏移成