丁健華，宛新林，賀曉華（安徽建筑大學(xué)土木工程學(xué)院，安徽合肥，230601）

?

一種用于波速測(cè)試的正演和反演算法

丁健華，宛新林，賀曉華
（安徽建筑大學(xué)土木工程學(xué)院，安徽合肥，230601）

摘要：波速測(cè)試的正演和反演算法是通過(guò)地表或井間觀測(cè)到的地震波走時(shí)，推算地下或井間的介質(zhì)速度結(jié)構(gòu)和地震波的傳播路徑。其步驟如下：（1）數(shù)據(jù)采集。（2）將地下介質(zhì)離散成一定大小的網(wǎng)格并初始化慢度（速度的倒數(shù)）場(chǎng)，采用射線追蹤技術(shù)推導(dǎo)各網(wǎng)格節(jié)點(diǎn)的射線路徑和最小旅行時(shí)。（3）建立線性方程組。（4）反演求解線性方程組，并修正初始慢度模型。（5）重復(fù)（2）-（4）步驟，不斷修正慢度場(chǎng)，直到滿足精度要求為止。本文采用最短路徑射線追蹤正演算法和同時(shí)迭代重建（SIRT）反演算法在二維平面上進(jìn)行求解計(jì)算，并且把得到的結(jié)果與傳統(tǒng)單孔法比較，從而了解本文所采用的波速測(cè)試方法的可靠性、實(shí)用性及其精度情況。

關(guān)鍵詞：最短路徑射線追蹤；SIRT算法；單孔法

0　引　言

土層波速，特別是剪切波速的研究具有重要意義，因?yàn)橥翆蛹羟胁ㄋ訇P(guān)系到場(chǎng)地類別的劃分、關(guān)系到砂土液化現(xiàn)象和土層的動(dòng)剪切模量等。目前土層波速測(cè)試主要是分為實(shí)驗(yàn)室測(cè)試和現(xiàn)場(chǎng)原位測(cè)試。前者主要應(yīng)用于場(chǎng)地復(fù)雜、無(wú)法進(jìn)行原位測(cè)量的重要工程?，F(xiàn)場(chǎng)測(cè)試波速的方法主要有單孔法、跨孔法和表面波法等。單孔法具有設(shè)備簡(jiǎn)單，操作方便，是應(yīng)用的較多的原位測(cè)試技術(shù)，但是單孔法假設(shè)地下介質(zhì)為水平層狀結(jié)構(gòu)和地震波沿直線路徑傳播，這都忽略了地下介質(zhì)非均勻分布的影響，從而導(dǎo)致計(jì)算結(jié)果的誤差。因此能不能模擬出地震波的真實(shí)傳播路徑，能不能得到地下非均勻介質(zhì)各點(diǎn)處的波速（即波速場(chǎng)），這就是本文所研究的內(nèi)容。目前關(guān)于這方面的研究很多，也取得了很大進(jìn)展。其中大部分方法都采用了將地下介質(zhì)離散成網(wǎng)格并在網(wǎng)格上布置節(jié)點(diǎn)，然后，分兩步求解：一是正演計(jì)算，即計(jì)算每個(gè)網(wǎng)格節(jié)點(diǎn)的最小旅行時(shí)和追蹤節(jié)點(diǎn)的地震波傳播路徑；二是反演求解，修正模型。在求解過(guò)程中需不斷地重復(fù)這兩步，直到前后兩次的慢度差滿足精度要求為止。網(wǎng)格中的正演求解是通過(guò)射線追蹤的方法來(lái)實(shí)現(xiàn)的。

目前，用于射線追蹤的算法有很多，如傳統(tǒng)的試射法［1］和彎射法［2］、最短路徑法、有限差分求解程函方程法［3］、旅行時(shí)線性插值法（LTI）［4］等。反演求解線性方程組的方法有最小二乘法（QR）、代數(shù)重建法（ART）、同時(shí)迭代重建法（SIRT）、分解法（LSQR）［5］等多種算法。各種算法都有各自的優(yōu)缺點(diǎn)，本文采用的正演算法是最短路徑射線追蹤算法，它是一種穩(wěn)定性較好的算法；反演算法采用了SIRT算法，它能夠較為準(zhǔn)確的反應(yīng)地下介質(zhì)的波速分布情況。文中采用了MATLAB編程技術(shù)，在二維平面上進(jìn)行地震波波速的正演和反演計(jì)算，得到所有網(wǎng)格節(jié)點(diǎn)的波速和地震波傳播路徑，取接收點(diǎn)所在水平線上的所有節(jié)點(diǎn)波速平均值作為該層的波速值，然后，將得到的波速結(jié)果與單孔法比較，以確定該方法的可靠性、實(shí)用性及其精度。

1　射線追蹤

1.1射線追蹤方法原理

最短路徑射線追蹤方法的理論基礎(chǔ)是光學(xué)中的Fermat原理和圖論中的最短路徑理論［6］。將地下非均勻介質(zhì)離散成若干個(gè)大小相等的小單元體，并認(rèn)為這些小單元是均勻的，然后在各單元邊界上設(shè)置一些離散節(jié)點(diǎn)，求解出這些節(jié)點(diǎn)處的速度，地下介質(zhì)的速度模型就可近似的由這些節(jié)點(diǎn)速度代替。在網(wǎng)格中，地震波從激發(fā)點(diǎn)到接收點(diǎn)一般有很多條路徑。每條從激發(fā)點(diǎn)到接收點(diǎn)的路徑都是由相互連接的節(jié)點(diǎn)序列組成的，相鄰節(jié)點(diǎn)之間的旅行時(shí)為節(jié)點(diǎn)之間直線距離與兩相鄰節(jié)點(diǎn)的平均慢度之積，整條路徑上總的旅行時(shí)即為該路徑所有相鄰節(jié)點(diǎn)射線路徑旅行時(shí)之和。根據(jù)Fermat原理，把從激發(fā)點(diǎn)到接收點(diǎn)傳播時(shí)間最短的路徑近似為地震波傳播路徑，該路徑即為最短路徑射線追蹤得到射線路徑。

1.2網(wǎng)格模型的建立

網(wǎng)格根據(jù)需要可以有多種形狀，在二維平面中可以是三角形、平行四邊形、矩形等，三維空間中可以采用四面體、立方體單元等。本文研究的是二維平面，在二維網(wǎng)格模型［7］中一般采用矩形網(wǎng)格，將地下介質(zhì)離散化為m×n個(gè)大小相同的矩形網(wǎng)格，地下波速場(chǎng)就可以近似的以這些網(wǎng)格單元的波速表示。在網(wǎng)格邊界上，節(jié)點(diǎn)的設(shè)置存在著多種方法。如可以將節(jié)點(diǎn)僅布置在角點(diǎn)上，也可將接節(jié)點(diǎn)只規(guī)則的布置于角點(diǎn)之間的邊界上，當(dāng)然，還可以在角點(diǎn)與邊界上同時(shí)布置節(jié)點(diǎn)，本文采用將節(jié)點(diǎn)僅布置于角節(jié)點(diǎn)上的方法。網(wǎng)格上每個(gè)節(jié)點(diǎn)與周圍節(jié)點(diǎn)連接方式取決于有效半徑r：r=1時(shí)，網(wǎng)格節(jié)點(diǎn)可以與周圍8個(gè)節(jié)點(diǎn)相連接，r=2時(shí)，網(wǎng)格節(jié)點(diǎn)可以與周圍16個(gè)節(jié)點(diǎn)相連接，r=3時(shí)，網(wǎng)格節(jié)點(diǎn)可以與周圍16個(gè)節(jié)點(diǎn)相連接，以此類推，還可以有更多的連接方式，這樣的連接方式可以避免在同一方向上射線經(jīng)過(guò)多個(gè)網(wǎng)格節(jié)點(diǎn)的情況。實(shí)際上并非r越大越好，從圖1中可以看出，r越大，射線就可能會(huì)經(jīng)過(guò)更厚的地下介質(zhì)，由于地下介質(zhì)的非均勻性，射線在厚的介質(zhì)中沿直線傳播的可能性很小。所以應(yīng)根據(jù)實(shí)際情況選擇合適的有效半徑，本文采用r=2的連接方式。

圖1　網(wǎng)格模型

1.3最短路徑射線追蹤的算法

由惠更斯原理［8］可以知道，波陣面上的點(diǎn)可以視為次級(jí)波源，向外發(fā)射次級(jí)波。因此可以將網(wǎng)格上的節(jié)點(diǎn)作為次級(jí)波源。再由Fermat原理和圖論中的Dijkstra算法［9］求出地震波從震源傳播到所有網(wǎng)格節(jié)點(diǎn)上的最小旅行時(shí)［10］。實(shí)際上，整個(gè)射線追蹤過(guò)程就是不斷向前追蹤次級(jí)波源的過(guò)程。具體步驟如下：

（1）從震源p開(kāi)始，計(jì)算與震源相連的節(jié)點(diǎn)A1-A16的旅行時(shí)，從中找出旅行時(shí)最小的節(jié)點(diǎn)A1作為次級(jí)波源（如圖2所示），并記錄A1節(jié)點(diǎn)的最小旅行時(shí)。

（2）計(jì)算與A1相連的所有節(jié)點(diǎn)的旅行時(shí)（震源P除外），當(dāng)相連節(jié)點(diǎn)為已算過(guò)旅行時(shí)節(jié)點(diǎn)時(shí)，需重復(fù)計(jì)算該節(jié)點(diǎn)旅行時(shí)，并將算出的最小的旅行時(shí)替換為該節(jié)點(diǎn)的旅行時(shí)，然后在所有已算過(guò)旅行時(shí)的節(jié)點(diǎn)（除震源P和已做過(guò)次級(jí)波源的節(jié)點(diǎn)A1）中找出旅行時(shí)最小的節(jié)點(diǎn)A13，把A13作為次級(jí)波源（如圖3所示），并記錄A13節(jié)點(diǎn)的最小旅行時(shí)。

圖2　從震源開(kāi)始追蹤次級(jí)波源

圖3　從次級(jí)波源繼續(xù)搜索下一次級(jí)波源

（3）重復(fù)步驟（2），不斷地搜索下一個(gè)次級(jí)波源，直到網(wǎng)格所有節(jié)點(diǎn)都做了次級(jí)波源，都計(jì)算出最小旅行時(shí)為止。在搜索次級(jí)波源過(guò)程中應(yīng)同時(shí)記錄每一次級(jí)波源的前一級(jí)次級(jí)波源。

（4）從接收點(diǎn)開(kāi)始逆向追蹤得到從震源到各接收點(diǎn)的射線路徑。

2　生成線性方程組

反演線性方程組是以矩陣形式給出，AS=T，其中A是m行n列的矩陣，其中m為射線路徑的條數(shù)，即接收點(diǎn)的個(gè)數(shù)，n等于網(wǎng)格節(jié)點(diǎn)總數(shù)。A中的i行j列元素有四種可能：若第i條射線路徑未通過(guò)j節(jié)點(diǎn)，則該元素為零；若j節(jié)點(diǎn)為第i條射線的起始節(jié)點(diǎn)，則該元素為射線路徑起始節(jié)點(diǎn)到路徑第二節(jié)點(diǎn)連線長(zhǎng)度的一半；若j節(jié)點(diǎn)為射線路徑起始節(jié)點(diǎn)與末端節(jié)點(diǎn)之間的點(diǎn)，則此元素為該節(jié)點(diǎn)與路徑后一節(jié)點(diǎn)射線長(zhǎng)度和該節(jié)點(diǎn)與路徑前一節(jié)點(diǎn)射線長(zhǎng)度和的一半；若j節(jié)點(diǎn)為路徑末端，則該元素為路徑末端節(jié)點(diǎn)和前一節(jié)點(diǎn)射線長(zhǎng)度的一半。S為每個(gè)節(jié)點(diǎn)慢度構(gòu)成的n維向量，S（j）為對(duì)應(yīng)于A中第j列節(jié)點(diǎn)j的慢度。網(wǎng)格中的節(jié)點(diǎn)可以用矩陣索引的方式唯一表示，例如i行第j個(gè)元素可以用來(lái)表示第i行第j個(gè)節(jié)點(diǎn)，亦可用一個(gè)列向量唯一表示節(jié)點(diǎn)，兩者之間存在著對(duì)應(yīng)關(guān)系，即第i行第j列的節(jié)點(diǎn)與向量中的第k個(gè)元素對(duì)應(yīng)，k=（i-1）*N+j ，N為節(jié)點(diǎn)矩陣的列數(shù)。

3　SIRT算法反演方程組

SIRT算法是Simultaneous Iterative Reconstruction Techniques的簡(jiǎn)稱，譯為同時(shí)迭代重建技術(shù)。算法步驟如下：

（2）利用初始慢度值求出每個(gè)方程的旅行時(shí)估計(jì)值：

（i=1，2，…，m。m為射線路徑條數(shù)。）

（4）計(jì)算每個(gè)方程對(duì)網(wǎng)格中節(jié)點(diǎn)的慢度修正值：如第i個(gè)方程對(duì)第j 個(gè)節(jié)點(diǎn)的慢度修正值為

（5）求慢度修正值。將每個(gè)方程對(duì)節(jié)點(diǎn)的慢度修正值求平均，就得到節(jié)點(diǎn)的慢度修正值，由于有很多節(jié)點(diǎn)沒(méi)有射線通過(guò)，其慢度值便沒(méi)有得到修正，對(duì)于這樣的節(jié)點(diǎn)可通過(guò)插值的方法獲得其修正。具體方法是一該節(jié)點(diǎn)為中心，通過(guò)某一半徑r范圍內(nèi)的已修正的節(jié)點(diǎn)的慢度修正值對(duì)其加以修正：

（其中，di為節(jié)點(diǎn)i到所求節(jié)點(diǎn)的距離，k為所求節(jié)點(diǎn)半徑r范圍內(nèi)的已修正節(jié)點(diǎn)個(gè)數(shù)）。

（8）重復(fù)上述步驟，直到前后兩次的慢度差滿足精度要求為止。

4　理論模型與算例

為了驗(yàn)證本文方法具有比傳統(tǒng)單孔法更高的計(jì)算精度，進(jìn)行了水平層狀理論模型下的波速仿真計(jì)算。在此基礎(chǔ)上，將本文方法用于工程實(shí)例計(jì)算，進(jìn)一步驗(yàn)證本文方法在實(shí)際中的可靠性與穩(wěn)定性。

4.1水平層狀模型

模型參數(shù)如下：水平方向1.5m，豎向深度10m，模型分為三層：第一層深2m，波速為180m/s；第二層深3m，波速為220m/s；第三層深5m，波速為260m/s。采用單孔法進(jìn)行原位測(cè)試，激發(fā)點(diǎn)位于（0，0）處，接收點(diǎn)坐標(biāo)分別為（1.5，1），（1.5，2），…，（1.5，10）處（模型如圖4所示）。分別用單孔法和最短路徑射線追蹤與SIRT反演算法求得地下介質(zhì)波速，波速結(jié)果如表1所示。

圖4　水平層狀模型

表1　波速結(jié)果

表2　波速和接收點(diǎn)走時(shí)

4.2 工程實(shí)例

本文以蕪湖市某擬建工程為例。該工程總占地面積約2萬(wàn)m2，總建筑面積約為60000m2，其中地下建筑面積約為10000m2。根據(jù)現(xiàn)場(chǎng)鉆孔和場(chǎng)地情況，選取了現(xiàn)場(chǎng)9號(hào)鉆孔進(jìn)行了剪切波試驗(yàn)。該試驗(yàn)將激發(fā)點(diǎn)置于離孔邊1.5米處，接收點(diǎn)分別位于孔內(nèi)1m、2m、3m、4m、5m、6m、7m、8m、9m、10m、11m、12m、13m、14m、15m、16m、17m、18m、19m、20m處。分別用傳統(tǒng)的單孔法和最短路徑射線追蹤與SIRT反演算法求得地下介質(zhì)波速和接收點(diǎn)走時(shí)，兩種方法得到的結(jié)果，如表2所示。

5　結(jié)論

（1）由表1數(shù)據(jù)可知，與傳統(tǒng)單孔法相比，本文所采用的方法具有更高的精度。另外，從表2可看出由該方法得到的接收點(diǎn)的走時(shí)近似等于現(xiàn)場(chǎng)拾取的走時(shí)，由Fermat原理知道該算法得到的射線路徑逼近真實(shí)的地震波傳播路徑。

（2）從表中數(shù)據(jù)可看出最短路徑法與SIRT算法結(jié)合使用得到的波速與單孔法測(cè)試結(jié)果總體相差不大，它們都能夠得到滿足工程精度要求的波速結(jié)果。表中數(shù)據(jù)還表明文中所用算法程序具有穩(wěn)定性和可靠性。

（3）理論上，網(wǎng)格劃分的越細(xì)，越能逼近地震波傳播的真實(shí)路徑，得到的結(jié)果越精確，但網(wǎng)格劃分的越細(xì)，節(jié)點(diǎn)數(shù)越多，計(jì)算量越大，從而導(dǎo)致程序運(yùn)行的時(shí)間過(guò)長(zhǎng)。實(shí)際應(yīng)用中可以根據(jù)場(chǎng)地實(shí)測(cè)資料和工程精度要求來(lái)選擇合適的網(wǎng)格大小。另外，還可以增加激發(fā)點(diǎn)和接收點(diǎn)的數(shù)目，使更多的節(jié)點(diǎn)有射線通過(guò)，以此來(lái)提高波速計(jì)算的精度。

（4）該方法考慮了地下介質(zhì)的非均勻性，不僅能夠得到地下介質(zhì)波在地下介質(zhì)中的傳播路徑。

參考文獻(xiàn)

［1］Julian B R， Gubbins D. Three-dimensional seismic ray tracing［J］. J Geophys， 1977，43：95-114.

［2］Thurber C H， Ellsworth W L. Rapid solution of ray tracing problems in heterogeneous media［J］. Bulletin of the Seismological Society of America， 1980，70（4）：1137-1148.

［3］Vidale J. Finite-difference travel time calculation［J］. Bulletin of the Seismological Society of American， 1988，78：2062-2076.

［4］趙改善，郝守玲，楊爾皓，等. 基于旅行時(shí)線性插值的地震射線追蹤算法［J］. 石油物探，1998，37（2）：14-24.

［5］Paige C， Saunders M. LSQR： An algorithm for sparse linear equations and sparse least squares［J］. Association of Computational Mechanics Transactions and Mathematical Software. 1982，8：43-71.

［6］李志輝，劉爭(zhēng)平. 最短路徑法射線追蹤的MATLAB實(shí)現(xiàn)［J］. 工程地質(zhì)計(jì)算機(jī)應(yīng)用， 2004，36（4）：16-21.

［7］魏亦文. 基于復(fù)雜介質(zhì)最短路徑射線追蹤層析成像正演研究［D］. 湖南：中南大學(xué). 2007.

［8］陸基孟，王永剛. 地震勘探原理［M］. 北京：中國(guó)石油大學(xué)出版社. 2011.

［9］張建中，陳世軍，余大祥. 最短路徑射線追蹤方法及其改進(jìn)［J］. 地球物理學(xué)進(jìn)展， 2003，18（1）：146-150.

［10］喻維秋. 初至波表層模型層析反演及其應(yīng)用［D］. 成都：成都理工大學(xué). 2007.

A Forward and Inversion Algorithm for Wave Velocity Test

DING Jianhua，WAN Xinlin，HE Xiaohua
（School of Civil engineering ， Anhui Jianzhu University ， Anhui Hefei，230601）

Abstract：Forward and inversion algorithm of wave velocity test to calculate the underground or medium of crosshole velocity structure and seismic wave propagation path through traveling time of the earth's surface or crosshole seismic earthquake wave.The steps are as follows:（1）Data collection.（2）The underground medium must be divided into the same size of the grid and initializes the slowness （the reciprocal of speed），Ray tracing technology is used to derive the ray path and minimum travel time of the grid nodes.（3）Establishing linear equations（4） solving linear equations by using inversion techniques， then correcting the initial model of the slowness.（5）Repeating steps （2） to （4），Continuously revising slowness， until meet the precision requirement.This paper uses the shortest path ray tracing algorithm and the improved simultaneous iterative reconstruction techniques（SIRT） to fnd the result in 2 d plane ，and compares the results to- that got with the traditional single hole method so as to understand the reliability，practicality and accuracy of wave velocity test method in this article.

Keywords：shortest path ray tracing；SIRT algorithm；puckering method

收稿日期：2015-07-15

DOI:10.11921/j.issn.2095-8382.20160206

中圖分類號(hào)：TP319

文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：2095-8382（2016）02-024-06