連續迭代計算方法在震波走時成像中的應用

2018-06-01 10:51:00李步青李建樓

赤峰學院學報·自然科學版 2018年5期

李步青，李建樓

（1.宿州學院信息工程學院；2.宿州學院資源與土木工程學院，安徽宿州 234000）

地震走時成像技術是一項系統而復雜的工程，由于它在異常體探測中具有重要的作用而不斷得到發展.震波走時成像技術在震波成像中具有獨特的優勢，主要是因為震波走時便于測量和拾取，利用適合的路徑算法能夠快速、高效得到慢度分布特征，從而圈定異常體位置.但是，由于成像算法不同，反演得到的圖像也不同.迭代法是解線性方程組的一個很重要的方法，適合解系數矩陣為稀疏矩陣的大型線性方程組[1,2].本文擬在Matlab軟件中采用直接法和連續迭代法實現震波走時成像數值模型的反演和對比，并利用震波走時成像工程案例反映連續迭代法的成像效果，希望能為震波走時成像初學愛好者提供一定的借鑒.

1 震波走時成像基本理論

震波走時成像是典型的地球物理反演問題，涉及以下幾個方面：①模型參數化，②正演（射線追蹤），③反演，④解的評價[3].首先要建立模型，布置觀測系統，用分塊法或節點法模擬局部小異常，把物理問題轉化為求解一個大型、稀疏的、且往往是病態的線性方程組問題，轉化為矩陣運算問題[4]；為了得到適當的解,網格劃分應該合理，以震源點間距或檢波點間距劃分網格大小較為合適，可提高成像精度，并非網格劃分越小越好[5]；另外，要從數學上剔除病態問題,并把更多的先驗信息引入反演過程[6].模型參數化之后開始正演，即根據射線路徑計算理論走時.計算理論走時可以根據文[7]提出的一種方法,即假定震源到檢波點路徑為直線,從而可以記錄每條射線穿過的單元和統計每個網格單元穿過的射線數目等，從而為反演做準備.正演是為反演做準備的，正演為反演提供了理論震波走時.反演則是利用震波的理論到時和實際到時以及對應的系數矩陣反向求解慢度值.反演方法較多，例如直接反投影法是對每條射線都按照像元中射線穿過長度占整個射線長度的比例，把走時殘差分配到每個像元中，然后像元內走時殘差之和除以像元內射線段總長度即為該像元慢度值，不進行迭代，不存在發散問題，但成像分辨率低[8]；代數重建法將每一條射線走時差平均地分配給通過的單元,而不考慮像元內射線的長短，是按射線依次修改像元慢度值的一類迭代算法，這種方法主要問題是收斂性不好.連續迭代法是在某一輪迭代中, 所有像元的值都用前一輪的迭代結果來修正[9,10].連續迭代法還可以考慮射線彎曲的影響,此時只需在計算理論走時值之前插入射線追蹤過程即可[11].

2 震波透射走時成像數值模型及分析

圖1 震波走時成像數值模型

由圖1可以看出，觀測范圍劃分為4個網格，X1、X2、X3和X4分別表示四個網格的慢度值；使用6個發射點、6個接收點，1發1收方式，合計6條射線，對應的旅行時分別為t1～t6.

每條射線的傳播時間可以用下列式子表示：

寫成矩陣的形式，即為A*X=T

其中，A為系數矩陣，X為慢度矩陣，T為傳播時間矩陣.本例中，

假設我們不知道慢度分布情況，不妨根據經驗設圖1中的慢度均為100，即

令直達波理論到時為T0，則

現在給慢度一個擾動，慢度分布不再均勻，例如

按照最短走時路徑追蹤方法，估算得到的旅行時間為T=[180；100；160；120；198；198].

實際觀測到的旅行時應和估算值差別不大，可令觀測的直達波旅行時為 T=[180；110；170；130；210；200]；根據 A和T對慢度X進行直接法反演，反演結果為

其中，pinv(A)是MATLAB軟件中求矩陣偽逆的函數.

直接法反演后的慢度分布與實際慢度對比及誤差分析如表1：

表1 直接法反演慢度和實際慢度對比

由表1可以看出，利用直接法對慢度分布反演可以定性評價，但定量誤差較大.

現在仍然使用該系數矩陣，采用連續迭代方法進行慢度反演計算.在MATLAB軟件中使用迭代法對慢度逼近的過程代碼如下：

100次迭代反演后的慢度分布和實際慢度對比及誤差分析如表2：

表2 連續迭代反演慢度和實際慢度對比

由表1和表2對比可以看出，直接法計算最大誤差為18.45%，迭代法計算最大誤差為12.45%，迭代法將誤差大大縮小.因此，在觀測系統矩陣系數A不變的情況下，利用迭代運算方法能夠有效降低誤差，逼近真實慢度，提高反演精度.

3 應用案例

3.1 建立觀測系統物理模型

校園內的草地中間有一條素混凝土路面，在馬路的兩側按設計位置分別布置12個震源和12個檢波器，形成一個震波走時觀測系統.實驗現場觀測系統布置圖2：

圖2 路基觀測系統布置

3.2 模型參數化

根據觀測系統建立坐標系，劃分慢度網格單元（即像元），模擬射線路徑，根據炮點和檢波點坐標計算可以得到各單元內每條射線穿過它的長度.模擬的射線路徑如圖3：

圖3 路基探測射線模擬圖

3.3 震波數據采集和處理

震波數據采集使用人工激發的錘擊震源，信號采集儀器為KDZ1114-3型號礦井地質探測儀，檢波器采用TZBS系列100 Hz高阻尼傳感器.設置的主要采樣參數為：采樣頻帶0-4000Hz，采樣間隔120μs，采樣點數1024，固定增益0dB，超前采樣點數0；按照激發1次、12個檢波器同時接收的方式，逐點激發地震波，共可獲得144道地震記錄.每炮記錄（12道）按放炮順序進行拼接，拼接后的地震記錄如圖4：