基于彎曲射線的跨孔層析成像算法研究

第一章序言PAGEPAGE51內(nèi)容提要地學(xué)CT技術(shù)是80年代從醫(yī)學(xué)引入地學(xué)的一門新興地球物理探測技術(shù)，由于其具有探測范圍廣，分辨率高，得到的圖像結(jié)果直觀，便于工作者分析地下地層結(jié)構(gòu)的特點(diǎn)，因此這種技術(shù)已被廣泛的應(yīng)用于工程物探領(lǐng)域以及煤炭、石油、和金屬等資源勘探開發(fā)方面，因此研究此方法技術(shù)是具有理論和實(shí)際意義的。本文從CT技術(shù)原理出發(fā)，介紹了CT技術(shù)的正、反演的幾種方法。其中正演方法即射線追蹤包括平直射線的射線追蹤和彎曲射線追蹤，本文在matlab平臺上編制程序，實(shí)現(xiàn)了兩種方法，并對比了兩種方法，得到線性插值射線追蹤法比最短路徑法更精確，速度更快的結(jié)論。在層析成像反演方法中，重點(diǎn)介紹了LSQR算法及SIRT算法，利用給定的模型在平直射線追蹤時運(yùn)用LSQR算法進(jìn)行了層析反演，在彎曲射線追蹤時運(yùn)用SIRT算法進(jìn)行了層析反演，檢驗(yàn)了兩種算法。關(guān)鍵詞：層析成像、線性插值射線追蹤法、最短路徑法、最小二乘QR分解法、聯(lián)合迭代重建算法。目錄TOC\o"1-2"\h\z\u第一章序言 11.1、研究背景及意義 11.2、國內(nèi)外研究現(xiàn)狀及存在的問題 21.3、研究內(nèi)容 4第二章層析成像的基本理論 62.1、Radon變換 62.2、速度層析成像的基本原理 82.3、衰減層析成像的原理 10第三章層析成像中的正演理論 133.1、平直射線的射線追蹤 133.2、彎曲射線的射線追蹤 16第四章層析成像中的反演方法 354.1、代數(shù)重建算法（ART算法） 354.2、最小二乘QR分解法（LSQR算法） 374.3、聯(lián)合迭代重建算法（SIRT算法） 40第五章結(jié)論 49參考文獻(xiàn) 50發(fā)表論文 53致謝 54中文摘要 IAbstract III基于彎曲射線的跨孔層析成像的算法研究

第一章序言1.1、研究背景及意義隨著國民經(jīng)濟(jì)飛速發(fā)展，各種占地面積大、設(shè)計精巧的高層建筑物大量出現(xiàn)，如高架橋梁，摩天大廈等。科學(xué)技術(shù)的發(fā)展也促使地基基礎(chǔ)的設(shè)計、施工采用了許多新的工藝與方法。從而對地基勘察提出了許多更新、更高的要求。對于地下不明隱蔽物或喀斯特溶洞基巖構(gòu)造破碎帶等復(fù)雜地質(zhì)問題的工程評價更是地基勘探中的重要課題。目前除了采用超淺層高分辯人工地震和精密磁測方法外，可以采用高精度跨孔層析成像(CT)這一新技術(shù)進(jìn)行勘探[1]。CT（ComputerizedTomography）技術(shù)，即計算機(jī)層析成像，是由投影重建圖像的應(yīng)用技術(shù)之一。圖像重建是圖像處理的一個很重要的分支，目前已廣泛應(yīng)用于各領(lǐng)域中，其理論基礎(chǔ)源于奧地利數(shù)學(xué)家Radon于1917年發(fā)表的“論如何根據(jù)函數(shù)在流形上的積分來確定函數(shù)”一文。在該文中，Radon系統(tǒng)地論證了由積分值(即投影函數(shù)值projection)確定被積函數(shù)(即重建圖像函數(shù))的整套理論方法，CT技術(shù)的形成和發(fā)展提供了可靠的理論依據(jù)。美國物理學(xué)家A.M.Cormack在Radon理論基礎(chǔ)上做了進(jìn)一步研究工作，他于1963年在“JournalofAppliedPhysics”上發(fā)表文章“用線積分表示函數(shù)的方法及其在放射醫(yī)學(xué)上的應(yīng)用”，并且完成了仿真與實(shí)驗(yàn)研究，這標(biāo)志著由X射線投影重建圖像的解析數(shù)學(xué)方法的確立。1971年，由英國EMI公司工程師G.N.Hounsfield研制出的第一臺CT掃描裝置在英國一家醫(yī)院成功安裝。至此，CT技術(shù)在醫(yī)學(xué)上的成功應(yīng)用已轟動全球，其在成像方面無可爭辯的優(yōu)越性引起世人矚目。為此，1979年的諾貝爾生理和醫(yī)學(xué)獎授予了CT的發(fā)明者和物理學(xué)家A.M.Cormack和工程師G.N.Hounsfiled。此舉推動了CT技術(shù)的進(jìn)一步發(fā)展，越來越多的人進(jìn)行CT的理論研究和實(shí)驗(yàn)工作，不斷地對算法和技術(shù)手段進(jìn)行改進(jìn)，使得CT技術(shù)迅速地應(yīng)用于各個領(lǐng)域[2]。80年代末CT技術(shù)被引入地學(xué)中，國內(nèi)簡稱“CT”探測。目前國際上存在聲波（AT）、電磁波（EM）和電阻率（RT）三種地學(xué)“CT”途徑，是地球物理勘探技術(shù)的一項(xiàng)重點(diǎn)創(chuàng)新。其中，電磁波CT探測技術(shù)成像精度高、對鉆孔壁無破壞性、成本低，很適于在地基工程方面普遍推廣。對于復(fù)雜地質(zhì)條件，如地下溶洞、地下不明隱蔽物、不同風(fēng)化程度的基巖層、斷裂或斷層分布特征等。電磁波CT層析技術(shù)除了能揭示常規(guī)勘探折射法可望推斷的構(gòu)造特征外，還能進(jìn)一步直觀而精確地查明那些局部“異常病變體”的位置、形狀與性質(zhì)，這是傳統(tǒng)勘探方法難以解決而電磁波CT層析技術(shù)卻擅長的突出的優(yōu)點(diǎn)之所在。例如，對于地下巖溶洞的勘探，以往慣常優(yōu)先采用地質(zhì)鉆孔方法。當(dāng)然對鉆孔直接揭露的溶洞特征可以了解清楚，但由于巖溶洞在空間發(fā)育常作串珠狀或不規(guī)則狀，因此單純依靠鉆探來展示巖溶溶洞在整個平面上的發(fā)育特征，不僅成本昂貴，而且根本解決不了問題。若采用高分辨超淺層人工地震方法，往往也只能查明垂深向第一個溶洞的頂面埋深情況，而溶洞在不同深度也常呈串珠狀。因此也無法圓滿解決巖溶洞勘探這一復(fù)雜問題。采用電磁波CT層析技術(shù)，只要在相鄰兩個鉆孔中進(jìn)行探測，便可取得鉆孔之間的剖面圖像，結(jié)合已有的鉆孔地質(zhì)資料。便能精確判斷整個剖面的地質(zhì)構(gòu)造特征、溶洞形態(tài)及其它一些細(xì)部結(jié)構(gòu)。若采用足夠多的井間CT剖面，即能控制住整個場地地下巖溶洞的三維分布特征。對于重大工程和水利樞紐、大橋及高層、超高層建筑等，常要求搞清地下持力層(堅硬土層或基巖)是否存在斷裂構(gòu)造問題，這同樣可采用電磁波CT層析技術(shù)予以解決。目前，對于地基中基巖風(fēng)化的判斷由于受鉆孔數(shù)量的限制，幾乎無法完整地構(gòu)畫出持力層的平面分布圖。由于基巖風(fēng)化程度的不同，其孔隙度和飽和度有很大的差異，必然使其吸收系數(shù)值有很大的變化，從而在電磁波CT層析圖像顯示上有明顯差異，據(jù)此可區(qū)別風(fēng)化帶界限，為設(shè)計提供更為詳盡可靠的地質(zhì)資料。不過，電磁波CT層析成像技術(shù)在淺層水飽和軟土中的使用效果不明顯，因?yàn)橛倌噘|(zhì)土的水飽和度很大，而其吸收系數(shù)β值差異相對要小，使CT成像不太清楚，但對不同飽和度的老粘土、砂土、淤質(zhì)粘土、風(fēng)化層、基巖等會有不同的反應(yīng)而得到清晰的CT圖像[1]。1.2、國內(nèi)外研究現(xiàn)狀及存在的問題層析成像涉及到通過參數(shù)的觀測結(jié)果來求解參數(shù)的空間分布這一問題。在資源勘探及工程勘查領(lǐng)域，所調(diào)查的對象不是人體，而是地下的構(gòu)造情況，這一點(diǎn)雖然從基本原理及其理論背景上講有許多類似之處，卻和醫(yī)學(xué)界的層析成像有許多不同，如在輻射手段方面，醫(yī)學(xué)界用的是X射線或伽瑪射線，而地球物理學(xué)界用的是地震波、聲波、電磁波等；在醫(yī)學(xué)界射線的輻射路徑都是直線傳播，而在地球物理學(xué)界，因介質(zhì)變化波的傳播路徑會出現(xiàn)反射、折射、衍射等；醫(yī)學(xué)界的層析成像技術(shù)測定的參數(shù)是接受輻射強(qiáng)度，而地球物理學(xué)層析曾向技術(shù)測定的參數(shù)包括速度、吸收系數(shù)等。將地球物理層析技術(shù)應(yīng)用于工程中，是目前國際上的一項(xiàng)前沿性研究課題，涉及到理論方法、軟件系統(tǒng)、模擬試驗(yàn)、巖性測試、采集儀器、觀測技術(shù)以及資料分析等一系列重要環(huán)節(jié)，并且只有完整地解決這些問題，實(shí)際的工程應(yīng)用才成為可能。目前國外己研究和應(yīng)用較多的方法主要有電阻率層析成像、地震層析成像和電磁波層析成像三種，以達(dá)到直觀地透視地下構(gòu)造的目的(馮銳，林宣明，1996；楊文采，杜劍淵，1994；高星，2000；全亞榮、陳周平，2000;劉昌銼、孫振國，2000；趙大鵬,2001)。總的來說，又以電磁波層析成像(ElectromagneticWaveComputerizedTomography，簡稱EWCT)研究和應(yīng)用較多[3]。在我國，吳以仁(1982)等對EM法做了系統(tǒng)性的總結(jié)，易永森等(1990)還采用陰影交匯法取得了許多工程物探應(yīng)用成果，楊文采等對CT成像算法作了許多研究。到二十世紀(jì)九十年代初期，國家地震局地球物理研究所馮銳等對該技術(shù)在我國工程中的應(yīng)用做了大量開創(chuàng)性的工作，在國內(nèi)首先建立起來了巖石高頻電磁參數(shù)成套測試設(shè)備和技術(shù)，解決了一系列計算機(jī)軟硬件接口問題，并將多種理論方法實(shí)用化，提高了圖像重建算法的速度與精度，同時取得了大量工程應(yīng)用成果，其技術(shù)處于國際領(lǐng)先水平。就整體而言，國內(nèi)國際上EWCT技術(shù)處于急速發(fā)展的試驗(yàn)應(yīng)用研究階段。雖有許多應(yīng)用實(shí)例，但深度和廣度依然不夠，大量的空白點(diǎn)有待解決。盡管層析技術(shù)在醫(yī)學(xué)領(lǐng)域己經(jīng)取得了巨大的成功，但在地球物理學(xué)中仍存在著一些困難，諸如波長偏大、射線彎曲、非完全投影和成像目標(biāo)復(fù)雜等。為此，在觀測系統(tǒng)與數(shù)據(jù)處理中必須采取一些特殊的措施來彌補(bǔ)。1.3、研究內(nèi)容本文的主要研究內(nèi)容包括層析成像的正、反演。在層析成像技術(shù)中，當(dāng)?shù)叵陆橘|(zhì)速度變化不大時用直射線代替真實(shí)射線，可以取得較好的效果。在速度變化較大的情況下，實(shí)際傳播路徑與直線相差很大，用直射線追蹤已失去意義，必須研究彎曲射線的層析成像方法，即進(jìn)行射線追蹤。射線追蹤的方法有很多，如：試射法，基于旅行時的線性插值射線追蹤法和基于費(fèi)馬原理的最短路徑法等。試射法是最早提出和使用的射線追蹤方法，在數(shù)學(xué)上描述為初值問題。它是己知射線的初始點(diǎn)(激發(fā)點(diǎn))和初始出射的方向，求地震波傳播的射線路徑，具體做法是從震源點(diǎn)出發(fā)，給定一系列的出射角，按snell定律逐段追蹤計算扇形區(qū)射線束各射線的路徑，把滿足一定誤差條件的最靠近接收點(diǎn)位置的射線作為實(shí)際射線路徑。圖1.1線性插值示意圖線性旅行時插值方法[4]是把模型離散成均勻的正方形單元，旅行時和射線路徑的確定只與單元邊界的點(diǎn)有關(guān)。假設(shè)單元邊界上任一點(diǎn)的旅行時可由該邊界上相鄰兩個離散點(diǎn)的旅行時線性插值得到。如圖所示為一勻速正方形單元，單元邊界平行于坐標(biāo)軸，A、B是二個旅行時已知的點(diǎn)，要求取穿過A、B邊界到達(dá)D點(diǎn)的最小旅行時及射線路徑。設(shè)射線從C點(diǎn)通過，C點(diǎn)旅行時刻用線性內(nèi)插公式由A、B二點(diǎn)的值表示。D點(diǎn)的旅行時為C點(diǎn)旅行時與波在C、D間直線傳播時間之和，然后根據(jù)Fermat原理，就可以求出C點(diǎn)的位置（即射線路徑）和D點(diǎn)的最小旅行時，該方法分向前處理和向后處理，向前處理只計算各單元邊界上的節(jié)點(diǎn)的旅行時；向后處理根據(jù)Fermat原理追蹤射線路徑。確定的射線路徑不是單元邊界上離散節(jié)點(diǎn)的連線，而是穿過單元邊界上正好滿足最小旅行時條件的那一點(diǎn)。這樣，在一個單元內(nèi)射線路徑為一直線，同時邊界上的折射角隨入射角連續(xù)變化。Moser(1991)提出了基于Huygens原理和網(wǎng)絡(luò)理論的最短路徑射線追蹤方法(SPR)[5]。該方法將模型分為由弧線連接的節(jié)點(diǎn)構(gòu)成的網(wǎng)格，每個節(jié)點(diǎn)與相鄰接點(diǎn)相聯(lián)系。從激發(fā)點(diǎn)到所有節(jié)點(diǎn)的最短路徑構(gòu)成一最短路徑樹，每一射線節(jié)點(diǎn)即為繞射點(diǎn)，使能量不斷向前傳播。計算時從激發(fā)點(diǎn)開始，逐步向四周的單元擴(kuò)展，求出每個單元內(nèi)任意兩射線節(jié)點(diǎn)間的走時，并按Fermat原理確定出其最小走時和最短路徑，最后求出激發(fā)點(diǎn)到接收點(diǎn)的最短射線路徑和最小走時。在該方法中，波至能量從一個節(jié)點(diǎn)傳到另一個節(jié)點(diǎn)，能夠避免Vidale方法中失去一些潛在的最短射線路徑問題，但是，波至從一個節(jié)點(diǎn)傳至另一個節(jié)點(diǎn)的限制，勢必影響該方法的精度。如果節(jié)點(diǎn)較密，又需消耗大量的計算機(jī)時間和內(nèi)存。Fischer和Lees(1993)在單元邊界上用Snell定理改變射線方向，提出了在低速度區(qū)有效獲得正確射線的方法，使各單元上所需計算的節(jié)點(diǎn)數(shù)大大減少，從而提高了該方法的計算效率。本人是在前人的研究基礎(chǔ)上，重點(diǎn)討論了線性插值射線追蹤法和最短路徑射線追蹤方法，并編制程序?qū)崿F(xiàn)了這兩種方法。層析反演是跨孔層析成像的另一關(guān)鍵技術(shù)。成像方法包括對介質(zhì)波速度或波慢度成像稱為走時反演和對介質(zhì)視吸收系數(shù)成像稱為振幅反演。從算法上來講，振幅反演和走時反演都屬于求解積分方程的問題，可以化為相似的代數(shù)方程組求解。所用的算法是基本相同的。常見的方法有反投影技術(shù)（BPT），代數(shù)重建技術(shù)（ART），聯(lián)合迭代重建技術(shù)（SIRT），共軛梯度最小二乘法（CGLS），正交分解最小二乘法（LSQR）等[6]。本文主要運(yùn)用的是聯(lián)合迭代重建技術(shù)和正交分解最小二乘法進(jìn)行層析反演。

第二章層析成像的基本理論在地球物理勘探中，人們通過觀測包含地下地質(zhì)體物理特征信號的數(shù)據(jù)，推斷地下地質(zhì)體的形態(tài)特征。這種根據(jù)觀測數(shù)據(jù)推斷地下地質(zhì)體特性的工作過程就是所謂的“反演問題”。相反，“正演問題”就是在給定地下地質(zhì)體特征和特定的物理定律成立的前提下，確定所能記錄到的數(shù)據(jù)，層析成像就是正反演問題的求解過程?？缈讓游龀上袼惴煞譃閮纱箢悾夯诓▌臃匠痰牟ㄐ螌游龀上窈突谶\(yùn)動方程的射線層析成像。其中，射線層析成像僅利用了電磁波初至旅行時和波場衰減參數(shù)，方法原理簡單，干擾因素小，只要能充分利用可觀測空間和介質(zhì)的先驗(yàn)信息，采用誤差較小的反演算法，就可以獲得滿意的效果。因此，目前實(shí)際應(yīng)用中廣泛使用的仍然是射線層析成像法。2.1、Radon變換[7][8][9]層析成像的基本思想是利用物體外部的投影來重建物體內(nèi)部的參數(shù)分布。對不同的研究對象，可采用不同形式的發(fā)射源，利用層析成像技術(shù)獲得物體內(nèi)部不同的物性參數(shù)分布情況。在電磁波層析成像中，以電磁波作為發(fā)射源，求解介質(zhì)的速度分布。圖2.1直線L在兩坐標(biāo)系中的幾何關(guān)系從物體內(nèi)部圖像重建的角度看，一張物體切片圖像是兩個空間變量的圖像函數(shù)。從不同角度觀測目標(biāo)體，觀測到的波場信息應(yīng)是入射波方向和觀測點(diǎn)位置兩個變量的函數(shù)，稱其為投影函數(shù)。若L為平面上的任一直線，如圖2.1所示，則沿直線L的線積分：（2-1）式中：dl是直線L的線元素增量。對于直線L，其方程可表示為：或（2-2）于是(2-1)式的線積分又可以寫成關(guān)于t和的函數(shù)形式，即為的Radon變換，記作：（2-3）Radon變換與Fourier變化的關(guān)系記函數(shù)的Fourier變換為（2-4）令（2-5）在（2-4）中作代換并令即得（2-6）反演公式利用Radon變換與Fourier變換之間的關(guān)系式（2-6）及Fourier變換的反演公式（2-7）即可求得Radon變換的反演公式（2-8）上面的反演公式還可以進(jìn)步化為（2-9）由此可見，Radon變換的反演公式包括：求導(dǎo)，Hilbert變換及對求平均三種運(yùn)算。因此當(dāng)滿足適當(dāng)條件時（這在實(shí)際問題中一般是滿足的），由可以唯一確定，也就是說Radon變換的反演問題是存在且唯一的。但是要成為數(shù)學(xué)上適定的問題，還必須滿足穩(wěn)定性的要求；另外，反演公式（2-8）或（2-9）只有理論上的意義，并不適宜具體數(shù)值計算。2.2、速度層析成像的基本原理我們下面以速度層析為例介紹跨孔層析成像的原理[10][11][12]。圖2.2層析成像原理示意圖設(shè)圖2.2為某一地段的兩個井孔，假定在左端井孔放置發(fā)射天線右端放置接收天線且接收的是電磁波的穿透直達(dá)波，則有：=（2-10）其中：為介質(zhì)的波速，為波從激發(fā)點(diǎn)到接收點(diǎn)的傳播走時，為收發(fā)距，又稱為射線。(2-10)式可改寫為：=Ls(2-11)其中：，為波速度的倒數(shù)，常稱之為波慢度。對斷面作均等剖分，橫向依?x長度分割為M個網(wǎng)格單元，縱向依?y長度分割成N個網(wǎng)格單元。設(shè)：激發(fā)點(diǎn)的序號為i=1，2，……，N；接收點(diǎn)的序號為j=1，2，……，N；每個激發(fā)點(diǎn)與接收點(diǎn)均在網(wǎng)格邊線的中點(diǎn)上；網(wǎng)格單元序號為k=1，2，……K，（K=M*N），設(shè)每一格網(wǎng)格中的波的速度(或慢度)為一平均值.則可得到（2-12）為第i點(diǎn)激發(fā)時，在第j點(diǎn)接收到的波的走時；為第i點(diǎn)激發(fā)，第j點(diǎn)接收時射線通過網(wǎng)格單元k的路徑長度；為第k個網(wǎng)格單元的慢度。當(dāng)i為某一定值時，j可取遍1~N，即得N個（2-12）式的方程，構(gòu)成一個方程組：（2-13）即：（2-14）上式可寫成矩陣形式：（2-15）也可寫成其中：表示第i點(diǎn)激發(fā)時,得到的j個走時形成的矩陣，是觀測值；表示第i點(diǎn)激發(fā)時所有射線的所有段元形成的矩陣；表示所有網(wǎng)格中波傳播速度形成的矩陣，是待求量。當(dāng)i從1取到N時，則可得到一個個方程的方程組，可表示為：（2-16）因此在速度層析中是要對（2-17）求解，并將其轉(zhuǎn)化為一個的矩陣并以灰度或彩色圖像顯示出來。2.3、衰減層析成像的原理在做衰減層析成像時，幅度不能從線性積分直接獲得，但通過取對數(shù)，問題就可變成線性層析成像反演。探地雷達(dá)發(fā)射的一般是球面波，因此電場幅度隨距離的衰減可表示為：（2-18）式中表示天線輻射功率的一個常數(shù)；為天線的指向性；為介質(zhì)衰減常數(shù)。為電磁波的傳播距離，實(shí)測的射線振幅應(yīng)取決于發(fā)射天線和接收天線的指向性增益，如接收天線和發(fā)射天線輻射模式相同，則滿足：（2-19）式中：為沿射線方向電磁波振幅的衰減總量；是與系統(tǒng)有關(guān)的歸一化常數(shù)[4]；分別為發(fā)射天線和接收天線的方向性增益；為電磁波的傳播距離。假設(shè)傳播為遠(yuǎn)場傳播并且用的是偶極天線則有：（2-20）對上式取對數(shù)，整理得（2-21）其中：是偶極發(fā)射天線的輻射模式；是接收天線靈敏度函數(shù)；和分別是接收天線與發(fā)射天線和射線路徑之間的夾角。與速度層析相同將研究斷面分成若干網(wǎng)格單元，線性積分通過求和的方法來計算，則有：（2-22）其中為第條射線在第個單元中的長度；為第個單元的衰減系數(shù)。式（2-22）也可寫成矩陣形式：（2-23）解出此方程組即可得到各單元內(nèi)的衰減系數(shù)，從而了解地下介質(zhì)的情況。層析成像的關(guān)鍵就是從方程組（2-16）（2-23）求得地下介質(zhì)的速度分布或衰減系數(shù)分布。在直射線層析成像中，方程組（2-16）（2-23）的系數(shù)矩陣是容易求出的，且是不變的，使得該方程組成為一個線性方程組，可用迭代法求解次方程組得到慢度向量。在彎曲射線層析成像中，慢度向量（或衰減系數(shù)向量）為未知的，射線路徑也是未知的，即方程組（2-16）（2-23）的系數(shù)矩陣中的每個元素都是未知的，因此從方程組中解出慢度向量的問題是一個非線性問題。解決這個非線性問題一般做法是，現(xiàn)給定慢度向量的初始值（或衰減系數(shù)向量的初始值，通過射線追蹤（下一章將具體介紹）求出方程組的系數(shù)矩陣的近似矩陣，然后，僅將慢度向量（或衰減系數(shù)向量）看成未知向量，方程組就成為了線性方程組，用層析成像基本算法（第四章具體介紹）解此線性方程組可得到新的慢度向量（或衰減系數(shù)向量），在通過射線追蹤求出在新的慢度向量（或衰減系數(shù)）下的初至波旅行時（或振幅），根據(jù)此時的初至波旅行時（或振幅）與已知的初至波旅行時（或振幅）的差向量來修正慢度向量（或衰減系數(shù)向量），并將修正后的慢度向量（或衰減系數(shù)向量）作為新的初始值，在重復(fù)以上的過程，知道初至波旅行時（或振幅）誤差減到足夠小為之。第三章層析成像中的正演理論目前，人們對地球內(nèi)部的物理性質(zhì)（包括速度、密度、電導(dǎo)率、溫度等）以及礦產(chǎn)資源的分布已經(jīng)有了不同程度的了解。這些知識多數(shù)來自于地表地質(zhì)和地球物理、地球化學(xué)資料的反演和解釋，而不是來自于鉆井資料。對于物體內(nèi)部的詳細(xì)了解，一個有效的途徑就是依據(jù)對表層的探測采集獲取觀測數(shù)據(jù)，加強(qiáng)資料的處理和反演解釋，以獲取更接近真實(shí)的地球物理和地球化學(xué)模型。層析成像技術(shù)的反演理論的目的是根據(jù)觀測數(shù)據(jù)求取響應(yīng)的地球物理模型。因此首先必須確定觀測數(shù)據(jù)和地球物理模型參數(shù)之間的函數(shù)關(guān)系，使地球物理工作者既可以根據(jù)給定的模型參數(shù)計算相應(yīng)的觀測數(shù)據(jù)（即實(shí)現(xiàn)正演計算），也可以根據(jù)觀測數(shù)據(jù)求取地球物理模型的參數(shù)，實(shí)現(xiàn)反演映射。顯然正演是反演的前提和條件，只有解決了正演計算，不管是靠解析的方法還是數(shù)值的方法，才有可能實(shí)現(xiàn)反演映射。層析成像的正演要解決的問題就是找到接收點(diǎn)的旅行時及激發(fā)點(diǎn)與接收點(diǎn)之間的射線路徑[5][14]。當(dāng)?shù)叵陆橘|(zhì)較均勻時，電磁波可看做沿直線傳播，而當(dāng)?shù)叵麓嬖诓ㄋ佼惓］^大的非均勻體時，電磁波在地層中傳播時的射線彎曲現(xiàn)象必須加以考慮。我們要獲得地下構(gòu)造的清晰圖像，其關(guān)鍵環(huán)節(jié)是實(shí)現(xiàn)激發(fā)點(diǎn)和發(fā)射點(diǎn)之間的定位，即射線追蹤。3.1、平直射線的射線追蹤當(dāng)射線路徑為平直射線時，射線路徑應(yīng)為激發(fā)點(diǎn)與接收點(diǎn)之間的直線段，此時系數(shù)子矩陣的求取關(guān)鍵在于判斷此條射線與兩相鄰縱向線的交點(diǎn)的位置關(guān)系，并以此來判定網(wǎng)格內(nèi)是否有射線經(jīng)過以及射線經(jīng)過網(wǎng)格時，網(wǎng)格中線段元的長度。具體實(shí)現(xiàn)過程如下：假設(shè)從第i點(diǎn)發(fā)射，共掃描到N條射線，每條射線都與網(wǎng)格的縱向線有交點(diǎn)，根據(jù)這些射線的方程可得到這些交點(diǎn)的縱坐標(biāo)并形成一個交點(diǎn)矩陣其中j代表第j條射線，k代表從左邊算起的第k條縱向線，k從1取到M+1，即整個網(wǎng)格共有M+1條縱向線，i代表第i個發(fā)射點(diǎn)。令可得到每條射線與水平線的交角的正弦值組成的矩陣如圖3.1所示：設(shè)為小于等于且與最接近的橫向網(wǎng)格線的值。用if判斷語句判斷與的位置關(guān)系，其關(guān)系如下：圖3.1一條射線與一個網(wǎng)格單元兩邊界交點(diǎn)的位置關(guān)系（1）、此時網(wǎng)格中有值，且值（2）、且≠此時仍然只有網(wǎng)格中有值，且其值為（3）、此時網(wǎng)格和中都有值（4）此時網(wǎng)格和內(nèi)有值，且值為：（5）此時網(wǎng)格和以及兩網(wǎng)格之間的縱向網(wǎng)格內(nèi)都有值，且值為:在此需要用一個for循環(huán)令a從到則網(wǎng)格（6）此時同情況（5）相同但值有所變化在此也需要用一個for循環(huán)令a從到則網(wǎng)格將以上判斷i從1循環(huán)到N，j從1循環(huán)到N，k從1循環(huán)到M。在Matlab平臺上可得到N個矩陣，將這N個矩陣合到一起即是所要求得系數(shù)矩陣L。之后可利用層析反演來判定介質(zhì)中異常體的位置。圖3.215*30網(wǎng)格內(nèi)的平直射線的射線路徑3.2、彎曲射線的射線追蹤上述直射線追蹤層析成像適合于比較均勻的介質(zhì)，而在實(shí)際的地質(zhì)勘查中，地質(zhì)體一般為波速差異較大的非均勻介質(zhì)，這時用直射線追蹤反演誤差較大，所以需要采用彎曲射線追蹤。現(xiàn)有的表示復(fù)雜介質(zhì)的方法基本上有兩種：eq\o\ac(○,1)把復(fù)雜介質(zhì)離散成矩形或三角形，每個單元的速度（或慢度）為常數(shù)[15]；eq\o\ac(○,2)用矩形網(wǎng)格劃分復(fù)雜介質(zhì)，給定網(wǎng)格節(jié)點(diǎn)的速度（或慢度），其它地方的速度（或慢度）有這些節(jié)點(diǎn)的速度（或慢度）的線性插值表示。基于這些模型，形成了許多射線追蹤方法。傳統(tǒng)的射線追蹤方法，通常指試射法（Shootingmethod）和彎曲法（Bengingmethod）。試射法[16]是最早提出和使用的射線追蹤方法，在數(shù)學(xué)上描述為初值問題。它己知射線的初始點(diǎn)(激發(fā)點(diǎn))和初始出射方向，求波傳播的射線路徑。具體作法是給定一系列的出射角，從激發(fā)點(diǎn)發(fā)出一組到接收點(diǎn)附近的射線，把滿足一定誤差條件的最靠近接收點(diǎn)位置的射線作為實(shí)際的射線路徑或由靠近接收點(diǎn)的兩條射線的旅行時內(nèi)插出接收點(diǎn)的旅行時。但是，當(dāng)介質(zhì)速度結(jié)構(gòu)較復(fù)雜時，即使從激發(fā)點(diǎn)發(fā)出的射線很密，也很難確定激發(fā)點(diǎn)到所有接收點(diǎn)的射線路徑，而且計算費(fèi)時，又會形成射線屏蔽區(qū)。彎曲法[17]基于Fermat時間穩(wěn)定原理，屬數(shù)學(xué)上的二點(diǎn)邊值問題。從激發(fā)點(diǎn)到接收點(diǎn)傳播的路徑，是在真實(shí)射線路徑附近變動的所有路徑中能使旅行時最小或穩(wěn)定的路徑。若激發(fā)點(diǎn)和接收點(diǎn)給定，先給出射線路徑的初始猜測值，再用變分原理求泛函極值的有關(guān)算法，逐次迭代修正射線路徑，直至收斂于真實(shí)射線路徑。但彎曲法有時會陷入局部收斂，得不到全局極小解，難于處理速度的較大變化，而且計算效率仍然很低。傳統(tǒng)的射線追蹤算法往往會收斂到一個局部最小旅行時的路徑，而且計算成本太高，但如果收斂到全局極小，便會有較高的精度，點(diǎn)效率較低。為了克服傳統(tǒng)射線追蹤方法的上述缺點(diǎn)，從1990年代初以來，許多地球物理工作者致力于發(fā)展精度較高、效率較高而且實(shí)用的計算射線路徑和初至旅行時的波前射線追蹤方法。這些方法歸納起來主要有：旅行時線性插值方法、最短路徑法等。3.2.1、基于旅行時的線性插值射線追蹤方法1993年，Asawaka和Kawanaka等人針對井孔間的初至波射線追蹤，提出了基于Fermat原理的旅行時線性插值射線追蹤方法（LTI）。這種方法基于網(wǎng)格單元模型，適用于任意變速介質(zhì)，可以追蹤包括直達(dá)波、折射波、透射波等多種射線路徑。用該算法計算旅行時和追蹤射線路徑比其他常規(guī)方法更為精確，更加適應(yīng)速度變化大等復(fù)雜條件，這是一種很有發(fā)展前途的地震射線追蹤方法。旅行時插值射線追蹤方法把模型離散成均勻的正方形單元，在單元邊界上設(shè)置一些計算旅行時的節(jié)點(diǎn)，假設(shè)單元邊界上任一點(diǎn)的旅行時可由該邊界上相鄰兩個離散點(diǎn)的旅行時線性插值得到。計算時分兩步，第一步，從激發(fā)點(diǎn)開始，用插值公式計算各單元邊界上的節(jié)點(diǎn)的旅行時；第二步，從接收點(diǎn)開始，根據(jù)Fermat原理和插值公式確定射線路徑。該方法確定的射線路徑不是單元邊界上離散節(jié)點(diǎn)的連線，而是穿過單.元邊界上正好滿足最小旅行時條件的那一點(diǎn)，邊界上的折射角隨入射角連續(xù)變化。3.2.1.1、旅行時線性插值射線追蹤法的基本原理[18]。在非均勻介質(zhì)中進(jìn)行進(jìn)行射線追蹤時，需把介質(zhì)剖分成若干個均勻單元。在一個小單元內(nèi)，若假設(shè)地震波旅行時隨距離線性變化，那么，任意一點(diǎn)的旅行時可用其兩側(cè)2個已知點(diǎn)旅行時的線性插值來估算，再利用Fermat原理便可求出從這2個已知點(diǎn)之間穿過到達(dá)另一點(diǎn)的最小旅行時。圖3.3由二個已知旅行時點(diǎn)插值求令一點(diǎn)旅行時的關(guān)系幾何圖如圖3.3所示，A、B為2個已知點(diǎn)，旅行時分別為、，射線從A、B之間的C點(diǎn)通過到達(dá)D點(diǎn)。求C點(diǎn)的位置和D點(diǎn)的最小旅行時。設(shè)AB與AD線段之間的夾角為，AC距離為r，在A和B之間進(jìn)行旅行時線性插值，可得到C點(diǎn)的旅行時為：（3-1）這時D點(diǎn)的旅行時可表示為：（3-2）其中：、、均可通過A、B、D點(diǎn)的坐標(biāo)求取。根據(jù)費(fèi)馬原理，將上式對r求導(dǎo)，并令其等于零，可得：（3-3）（3-4）其中：此時C點(diǎn)坐標(biāo)可由下式確定：（3-5）由于射線從A、B之間通過，那么必須滿足約束條件。該式是在最小旅行時條件下射線穿過A、B之間的正確位置到達(dá)D點(diǎn)的約束條件。只有在此條件滿足的情況下，才能得到射線在A、B間通過的準(zhǔn)確位置C點(diǎn)的坐標(biāo)和到達(dá)D點(diǎn)最小旅行時。當(dāng)A、B連線平行于Z軸時如圖3.4，式（3-4）可簡化為：圖3.4由平行于z軸的連線上二個已知點(diǎn)插值求另一個點(diǎn)旅行時的幾何關(guān)系圖（3-6）此時C點(diǎn)坐標(biāo)為：（3-7）3.2.1.2、旅行時線性插值射線追蹤法的具體實(shí)現(xiàn)過程線性插值射線追蹤法是以線性旅行時假設(shè)為前提的，首先把要計算的模型劃分成等正方形的像元，在每個像元邊界上確定若干個計算點(diǎn)，假設(shè)邊界上任一點(diǎn)的旅行時是該邊界上相鄰兩離散點(diǎn)旅行時的插值然后進(jìn)行向前和向后處理。向前處理是指從激發(fā)點(diǎn)開始逐步計算出各個節(jié)點(diǎn)上的最小旅行時，具體步驟如下：圖3.5線性插值射線追蹤法向前處理示意圖（1）計算激發(fā)點(diǎn)所在單元邊界上各計算點(diǎn)的旅行時，如圖a所示(圖中S表示激發(fā)點(diǎn))；（2）計算激發(fā)點(diǎn)單元所在列各單元邊界上各計算點(diǎn)的旅行時，這時只考慮來自該單元下邊界，即激發(fā)點(diǎn)所在單元的上邊界的已知點(diǎn)，并利用公式求出計算點(diǎn)上的最小旅行時。再以該單元為激發(fā)單元，重復(fù)上述步驟，繼續(xù)向上計算，直到計算出激發(fā)點(diǎn)上方所有節(jié)點(diǎn)的最小旅行時為止。然后再以同樣的方法計算出激發(fā)點(diǎn)下方所有節(jié)點(diǎn)上的最小旅行時，如圖(b)所示。這樣，我們就得到了激發(fā)點(diǎn)所在列的所有單元邊界上節(jié)點(diǎn)的最小旅行時。計算激發(fā)點(diǎn)所在單元右側(cè)一列所有節(jié)點(diǎn)的最小旅行時。首先計算各單元水平邊界上各節(jié)點(diǎn)的最小旅行時，如圖(c)所示。這時只考慮來自左邊界的射線。對于任一水平單元來說，它既是上部單元的下邊界，又是下部單元的上邊界，這需要在二個單元中分別計算旅行時，再取二者中旅行時的最小者。在①中只考慮了來自單元左邊界的射線的情況，我們還應(yīng)該考慮來自單元上邊界或下邊界的射線的情況。如圖(c)所示，利用單元下邊界上的已知點(diǎn)計算上邊界上節(jié)點(diǎn)的最小旅行時;用單元上邊界上的已知點(diǎn)計算下邊界上節(jié)點(diǎn)的最小旅行時。通過①、②步就得到了該列所有單元水平邊界上各節(jié)點(diǎn)的最小旅行時。計算該列各單元右側(cè)邊界上的最小旅行時。這時要考慮來自單元另外三個邊界上所有的己知點(diǎn)，計算右邊界上各節(jié)點(diǎn)的最小旅行時，如圖(e)所示。重復(fù)第(3)步，直到計算出激發(fā)點(diǎn)右側(cè)各列所有節(jié)點(diǎn)的最小旅行時為止。向后處理是根據(jù)費(fèi)馬原理對全部激發(fā)點(diǎn)-接收點(diǎn)追蹤射線路徑。得到滿足最小旅行時條件下與像元邊界交點(diǎn)的射線路徑。在像元內(nèi)的射線路徑總是直線，同時像元邊界處的射線折射角隨入射角連續(xù)變化，其具體步驟如下：（1）首先在接收點(diǎn)所在單元，找出至接收點(diǎn)的旅行時為最小的邊界節(jié)點(diǎn)，如圖(a)所示。其中未知點(diǎn)為接收點(diǎn)，已知點(diǎn)為該單元邊界上所有節(jié)點(diǎn)，節(jié)點(diǎn)的旅行時已由前述的向前處理得到。然后根據(jù)互換原理式確定出該單元上節(jié)點(diǎn)至接收點(diǎn)的旅行時為最小的節(jié)點(diǎn)。（2）確定與上步找出的最小旅行時節(jié)點(diǎn)對應(yīng)的插值段。對這些插值段應(yīng)用，求出穿過各段至接收點(diǎn)的旅行時，其中的最小者就是該接收點(diǎn)對應(yīng)的初至旅行時，與該初至旅行時對應(yīng)的插值點(diǎn)便是所要求的初至射線路徑與單元邊界的交點(diǎn)，如圖(b)所示。圖3.6線性插值射線追蹤法向后處理示意圖（3）將求得的射線交點(diǎn)作為新的接收點(diǎn)，重復(fù)上述第(1)步和第(2)步，直至激發(fā)單元為止，如圖c、圖d所示。（4）將激發(fā)點(diǎn)和找出的左右交點(diǎn)及接收點(diǎn)順次相連就得到了從激發(fā)點(diǎn)到接收點(diǎn)的初至射線路徑，如圖e所示。圖3.7、圖3.8、圖3.9分別為不同模型下利用線性插值射線追蹤法得到的射線路徑?？梢钥闯觯涸谀Ｐ蜑榫鶆蚪橘|(zhì)的情況下，射線路徑接近直射線，當(dāng)存在高速異常的情況下，多條射線路徑向異常區(qū)域彎曲，而當(dāng)存在低速異常的情況下，本應(yīng)穿過異常區(qū)域的射線繞過此低速區(qū)域行進(jìn)。是基本符合Fermat原理與實(shí)際情況的。圖3.7均勻模型下的射線路徑圖3.8存在高速異常的模型下的射線路徑圖3.9存在低速異常的模型下的射線路徑（3）方法驗(yàn)證：圖3.10驗(yàn)證線性插值射線追蹤法示意圖圖3.10是一個的網(wǎng)格模型，背景速度為，第二行第四列和第三行第四列存在一個高速異常區(qū)，其速度為。通過線性插值射線追蹤我們可以得到一條穿過異常區(qū)的射線路徑，我們也可得到A、B、C點(diǎn)的坐標(biāo)，分別為A（2，2.686），B（3，2.3802），C（3.264，2），其中B點(diǎn)為射線穿過異常區(qū)的入射點(diǎn)，通過A、B、C三點(diǎn)的坐標(biāo)我們可以得到射線穿過B點(diǎn)時的入射角的正弦值為，折射角的正弦值為。則有而可以看出，是符合斯奈爾定律的。由此驗(yàn)證可知線性插值射線追蹤法是一種精度較高的射線追蹤法。3.2.2、最短路徑射線追蹤法最短路徑方法起源于網(wǎng)絡(luò)理論，首次由Nakanishi和Yamaguchi應(yīng)用于地震射線追蹤中。Moser以及Klimes和Kvasnicha對最短路徑方法進(jìn)行了詳細(xì)研究。通過科技人員的不斷研究，最短路徑方法目前已發(fā)展較為成熟，其基本算法的計算程序也較為固定。最短路徑射線追蹤方法基于Fermat最小旅行時原理和網(wǎng)絡(luò)理論中的最短路徑算法來實(shí)現(xiàn)。把地下介質(zhì)離散成若干小單元體。并在各單元邊界上設(shè)置一些節(jié)點(diǎn)。地下速度模型就可表示成由這些節(jié)點(diǎn)以及它們之間的連線所形成的網(wǎng)絡(luò)。網(wǎng)絡(luò)中的每一個節(jié)點(diǎn)只能與彼此相鄰的節(jié)點(diǎn)連接。相鄰節(jié)點(diǎn)之間的連接權(quán)等于電磁波沿其連線傳播的旅行時。一條路徑是由相互連接的節(jié)點(diǎn)序列組成的，沿著該路徑的旅行時為該路徑上所有連接權(quán)之和。從一個節(jié)點(diǎn)到另一個可能有無數(shù)條路徑，按照Fermat原理，把旅行時最小的路徑近似為電磁波傳播通過的射線。圖3.11最短路徑法的原理示意圖3.2.2.1速度解決復(fù)雜介質(zhì)的地球物理正、反演數(shù)值計算的有效途徑，是把復(fù)雜介質(zhì)離散成一系列小單元，使每個單元的幾何和物理參數(shù)變得簡單，以至于能進(jìn)行求解。以二維地質(zhì)模型為例，實(shí)際近地表模型是復(fù)雜的，需要對其進(jìn)行離散話處理。可采用四邊形單元、三角形單元、多邊形單元對地質(zhì)模型剖分，本文討論矩形單元剖分。如圖3.11，用由M條橫向直線和N條豎向直線組成的網(wǎng)絡(luò)把研究區(qū)劃分成個矩形單元。網(wǎng)格節(jié)點(diǎn)坐標(biāo)為、速度為，其中；。3.2.2.2本文在二維模型中采用矩形參數(shù)化方式，介質(zhì)模型都貝離散化為個大小相同的矩形網(wǎng)格。此時存在兩種節(jié)點(diǎn)設(shè)置的方法：模型一：將節(jié)點(diǎn)分布在矩形網(wǎng)格的四個角點(diǎn)上，同時波速也分配在節(jié)點(diǎn)上，如圖3.12所示，每個節(jié)點(diǎn)都可以與相鄰的節(jié)點(diǎn)連接，相連節(jié)點(diǎn)之間的走時為它們之間的歐式距離除兩個節(jié)點(diǎn)的平均慢度之積。一個節(jié)點(diǎn)與其周圍16個節(jié)點(diǎn)直接相連。這樣，一個激發(fā)點(diǎn)或下一級激發(fā)點(diǎn)可以有16個出射方向和16個與之對應(yīng)的計算節(jié)點(diǎn)；這種設(shè)置便于繪制旅行時和速度的等值線圖。更重要的是可以引入速度界面進(jìn)行反射波的射線追蹤。在實(shí)際應(yīng)用中，Moser等以及Klimes和Kvasnicha均采用了這種模型。(a)(b)圖3.12最短路徑法兩種節(jié)點(diǎn)設(shè)置中的模型一（2）模型二：網(wǎng)絡(luò)節(jié)點(diǎn)按一定規(guī)則分布在矩形單元的邊界上（不含單元角點(diǎn)），每一個矩形單元內(nèi)，速度或慢度為常數(shù)。只有在兩個節(jié)點(diǎn)之間不存在單元邊界時，這兩個節(jié)點(diǎn)才可以連接。如圖3.13所示：相連節(jié)點(diǎn)之間的旅行時由它們之間的歐式距離乘以它們所在的單元的慢度得到。這種模型由常速度單元組成，是兩個相連節(jié)點(diǎn)之間的旅行時的計算簡單而快捷并且在一個單元內(nèi)的射線追蹤是精確的。圖3.13最短路徑法兩種節(jié)點(diǎn)設(shè)置中的模型二3.2.2.3最短路徑射線追蹤法的具體實(shí)現(xiàn)過程計算從發(fā)射點(diǎn)出發(fā)到模型中所有節(jié)點(diǎn)i的最小走時，然后根據(jù)Fermat原理尋找射線路徑。從發(fā)射點(diǎn)到模型中任一個節(jié)點(diǎn)i的最小走時應(yīng)滿足Bellman方程式中表示從發(fā)射點(diǎn)到節(jié)點(diǎn)i的最小走時，表示相鄰節(jié)點(diǎn)j的最小走時表示i和j節(jié)點(diǎn)間的距離，c表示i和j節(jié)點(diǎn)間的慢度平均值。上式應(yīng)滿足的初始條件為，即發(fā)射點(diǎn)所在節(jié)點(diǎn)處走時為零。需要指出的是：1）對于復(fù)雜構(gòu)造的地質(zhì)模型必須精細(xì)的劃分，否則，誤差將較大；2）射線追蹤的精度主要取決于節(jié)點(diǎn)數(shù)M，精度要求越高，所需要的節(jié)點(diǎn)數(shù)越多。最短路徑射線追蹤法的具體實(shí)現(xiàn)步驟：把網(wǎng)格上所有節(jié)點(diǎn)分成集合P和Q，P為已知最小旅行時的節(jié)點(diǎn)總數(shù)集合，Q為未知最小旅行時的節(jié)點(diǎn)的集合。若節(jié)點(diǎn)總數(shù)為N，經(jīng)過N次迭代后可為求出所有節(jié)點(diǎn)的最小旅行時。（1）初始時Q集合包含所有節(jié)點(diǎn)，除激發(fā)點(diǎn)的旅行時為外，其余所有節(jié)點(diǎn)的旅行時均為：。圖3.14最短路徑法的實(shí)現(xiàn)過程示意圖對來說有；（2）在D中找一個旅行時最小的節(jié)點(diǎn)i，定它的旅行時為；（3）確定與節(jié)點(diǎn)i相連的所有節(jié)點(diǎn)的集合V；（4）求節(jié)點(diǎn)與節(jié)點(diǎn)i連線的旅行時；（5）求節(jié)點(diǎn)的新旅行時（取原有旅行時與的最小值）；（6）將i點(diǎn)從Q集合轉(zhuǎn)到P集合；（7）若P集合中的節(jié)點(diǎn)個數(shù)小于總節(jié)點(diǎn)數(shù)N，轉(zhuǎn)到步驟(2），否則結(jié)束追蹤；（8）從接收點(diǎn)開始倒推出各個發(fā)射點(diǎn)到接收點(diǎn)的射線路徑，只要每個節(jié)點(diǎn)記下使它形成最小旅行時的前一個節(jié)點(diǎn)，就很容易推出射線路徑。圖3.15-圖3.23為利用最短路徑法的到的射線路徑。由圖中可見射線全部都是通過計算點(diǎn)向前行進(jìn)的，其中在圖3.15-圖3.17中，節(jié)點(diǎn)設(shè)置是應(yīng)用了模型一，即節(jié)點(diǎn)設(shè)置在單元個邊界的角點(diǎn)上，圖3.15是均勻介質(zhì)情況下得到的射線路徑，圖3.16、圖3.17分別是從存在高速、低速異常區(qū)域的介質(zhì)模型中的到的射線路徑；圖3.18、圖3.19、圖3.20中節(jié)點(diǎn)設(shè)置都是運(yùn)用了模型二，即節(jié)點(diǎn)設(shè)置在單元邊界上，圖3.18是均勻情況下的射線路徑，圖3.19和圖3.20分別為存在高速異常和存在低速異常情況下得到的射線路徑，由圖可以看出，在最短路徑法中得到的射線路徑同實(shí)際射線路徑是有一定差異的，這與節(jié)點(diǎn)的設(shè)置有關(guān)，當(dāng)節(jié)點(diǎn)設(shè)置的更加密集時，是可以改善這種誤差的。圖3.15采用模型一節(jié)點(diǎn)設(shè)置時得到的均勻介質(zhì)下的射線路徑圖3.16采用模型一節(jié)點(diǎn)設(shè)置時得到的低速非均勻介質(zhì)下的射線路徑圖3.17采用模型一節(jié)點(diǎn)設(shè)置時得到的低速非均勻介質(zhì)下的射線路徑圖3.18采用模型二節(jié)點(diǎn)設(shè)置時得到的均勻介質(zhì)下的射線路徑圖3.19采用模型二節(jié)點(diǎn)設(shè)置時得到的存在低速異常介質(zhì)下的射線路徑圖3.20采用模型二節(jié)點(diǎn)設(shè)置時得到的存在高速異常介質(zhì)下的射線路徑圖3.21單元邊界上設(shè)置三個節(jié)點(diǎn)的情況（均勻介質(zhì)）圖3.22單元邊界上設(shè)置三個節(jié)點(diǎn)的情況（存在低速異常的介質(zhì)）圖3.23單元邊界上設(shè)置三個節(jié)點(diǎn)的情況（存在高速異常的介質(zhì)）圖3.21-圖3.23是利用模型二最短路徑法得到的射線，即計算點(diǎn)設(shè)置在單元邊界上，其中每個單元邊界上設(shè)置三個計算點(diǎn)。與圖3.18-圖3.20對比可以看出當(dāng)單元邊界上的計算節(jié)點(diǎn)加多時，射線路徑更接近實(shí)際情況。即前文提的精度要求越高所需的計算點(diǎn)數(shù)越多。對比最短路徑射線追蹤法和線性插值射線追蹤法可以看出：1、最短路徑射線追蹤法中的射線路徑與單元邊界的交點(diǎn)都是之前按照一定規(guī)律設(shè)定好的節(jié)點(diǎn)，而線性插值射線追蹤法的射線路徑與單元邊界的交點(diǎn)是通過線性插值得到的，因此精度要高于最短路徑法。2、最短路徑射線追蹤法在要求精度較高的時候，要求在網(wǎng)格邊界上設(shè)置的節(jié)點(diǎn)要多，這時運(yùn)算速度就會下降，而線性插值射線追蹤法只在網(wǎng)格邊界上設(shè)置一個節(jié)點(diǎn)就可以得到較高精度的射線路徑，其運(yùn)算速度遠(yuǎn)快于最短路徑射線追蹤法。

第四章層析成像中的反演方法層析成像技術(shù)的關(guān)鍵在于求出方程組（2-16）的合理解，通過上一章介紹的射線追蹤法可以得到該方程組的系數(shù)矩陣，此方程組通常是大型、稀疏、強(qiáng)超定、欠定甚至不相容的方程組，所以要求求解此方程組的算法具有穩(wěn)定、收斂、節(jié)省內(nèi)存、效率高等特點(diǎn)。目前，適用的算法有很多，其中最為基本的算法有代數(shù)重建方法、共軛梯度方法、濾波反投影方法、SVD方法等。這里主要介紹代數(shù)重建方法（ART方法）和聯(lián)合迭代法（SIRT方法）以及結(jié)合Lanczos投影法、最小二乘法和矩陣的QR分解法而等到的LSQR算法。本文在層析反演過程中采用了后兩種方法，其中當(dāng)射線為平直射線時運(yùn)用了LSQR算法，而當(dāng)射線彎曲時采用了SIRT算法。4.1、代數(shù)重建算法（ART算法）ART（AigebraicReconstructionTechnique）方法是按射線依次修改有關(guān)象元的圖像向量的一類迭代算法。其迭代過程為：首先，給定慢度向量的初值，然后，循環(huán)地按照方程組的第一個方程到最后一個方程，依次對慢度向量進(jìn)行修正，知道修正后的慢度向量滿足預(yù)定誤差的要求為止。在方程中令圖像向量產(chǎn)生一個增量，有：（4-1）作為迭代算法要根據(jù)第j條射線的走時差求慢度的修改增量。由于方程可能是欠定的或病態(tài)的，可以用它作為約束的模的極小解。由拉格朗日乘子法令目標(biāo)函數(shù)為：(其中為拉格朗日乘子)，（4-2）由：得：?；卮胧接校杭匆虼?，再由上式便可寫出第j條射線及第i個像元求波慢度修改增量的公式：（4-3）上式就是所謂的加法修正的公式。當(dāng)然我們不一定非取的模極小，也可以取任意階的模，如模極小來求慢度的修改增量，其中。由向量范數(shù)定義式可知，此時目標(biāo)函數(shù)可取為：（4-4）令及最后可得：（4-5）此式雖然可以用來作迭代修改，但是當(dāng)時涉及開方運(yùn)算，速度太慢，一般很少采用。只是在令時，可導(dǎo)出ART迭代的最簡單修正公式：（4-6）它說明走時差平均地分配給每一條射線j通過的單元，而不考慮象元內(nèi)射線的長短。綜上所述，ART方法的具體步驟為：eq\o\ac(○,1)選定初值；eq\o\ac(○,2)計算第j個觀測值與第j個方程的估算值之差；eq\o\ac(○,3)計算慢度向量。其中k從1開始遞增。對于每個k都作eq\o\ac(○,2)、eq\o\ac(○,3)兩步。隨著k的遞增，其對應(yīng)的方程從第一個方程到最后一個方程逐輪循環(huán)。每完成一輪循環(huán)，要判定迭代結(jié)果滿足預(yù)定誤差的要求，滿足則停止，不滿足則進(jìn)入下一輪循環(huán)。每次完成第eq\o\ac(○,3)步后，需要對加以約束。4.2、最小二乘QR分解法（LSQR算法）在求出系數(shù)矩陣后即可用LSQR算法解線性方程組。LSQR方法是Paige和Sanders在1982年提出的[6]，它是利用Lanczos迭代法求解最小二乘問題的一種方法。LSQR方法具有計算量小的優(yōu)點(diǎn),并且能很容易地利用矩陣的稀疏性簡化計算,因而適合求解大型稀疏問題。對于方程，其最小二乘問題可以通過雙對角化來求解假定和是正交陣，且為如下的的下雙角陣（4-7）用下列迭代方法可實(shí)現(xiàn)矩陣的雙對角分解：（4-8）其中，。使（4-8）式又可寫成如下形式（4-9）其中表示階單位矩陣的第行，再設(shè)可以確定：（4-10）在滿足給定精度時停止迭代。由于我們希望盡量小，且理論上是正交陣，取使最小。解最小二乘問題這就構(gòu)成LSQR算法的基礎(chǔ)。LSQR主要步驟總結(jié)如下[8]：(1)初始化，，，，，其中、為m維向量，、為n維向量，、、、為實(shí)數(shù)。(2)對i=1,2,3,?作以下各步(3)雙對角化矩陣(a)(b)(4)修改參數(shù)(a)(b)(c)(d)(e)(f)(g)(5)迭代求解(a)(b)(6)收斂判別最簡單的測試就是，當(dāng)?shù)螖?shù)增加時，所求得的解沒有明顯變化就可以停止迭代。當(dāng)用LSQR法解出方程組得解后就可以將結(jié)果以灰度圖象顯示出來，從而可以直觀的辨別出研究斷面的異常體。圖4.1平直射線模型及反演成果圖圖4.1為運(yùn)用平直射線追蹤及LSQR算法得到的層析成像反演結(jié)果，可以看出，結(jié)果是比較理想的。圖4.2平直射線模型及反演成果圖圖4.2是一個存在高速，低速兩個異常區(qū)域的模型，背景速度為0.13m/ns,高速異常區(qū)為0.3m/ns,低速異常區(qū)為0.033m/ns。該圖同樣是利用平直射線的射線追蹤，再運(yùn)用LSQR算法進(jìn)行反演得到層析成像反演結(jié)果。兩個異常區(qū)域都很好的反演，呈現(xiàn)出來。4.3、聯(lián)合迭代重建算法（SIRT算法）最優(yōu)化準(zhǔn)則和聯(lián)合迭代重建技術(shù)(SIRT)也稱逐點(diǎn)重建法，最初是由Gilbert(1972)提出[41]-[42]。SIRT和ART類似，其區(qū)別在于：ART每次修正值只考慮一條射線，SIRT則是利用一個像素內(nèi)通過的所有射線的修正值來確定這一像素的平均修正值。SIRT是所有射線通過計算以后，才完成一次迭代。取平均修正值可以壓制一些干擾因素，而且SIRT的計算結(jié)果與觀測數(shù)據(jù)使用順序無關(guān)。SIRT算法原理是：首先給出一組速度初值，計算估算值，得到觀測值和估算值的差,其中。則在像素內(nèi)的平均修正值（4-11）式中：表示穿過第j個網(wǎng)格的射線數(shù)，n表示總的射線數(shù)，m表示模型離散后的網(wǎng)格數(shù)。用上式對第j個像素的慢度值進(jìn)行修正，得到當(dāng)估算值與觀測值之間的差小于給定誤差時停止迭代，否則進(jìn)行下一輪迭代。其具體步驟如下：首先給未知數(shù)一組初值，令，。計算第一到第M個方程的估算值，（4-12）式中的是由射線追蹤算法求取到的系數(shù)矩陣。求出觀測值與計算值之差（4-13）設(shè)第i條射線在第j個象元內(nèi)的第q次慢度修正值為，則（4-14）設(shè)修正值正比于第j個象元內(nèi)射線所通過的路徑與該射線長度之比，即（4-15）式中，是第i條射線的比例常數(shù)，是第i條射線的全長（4-16）將式（4-16）帶入式（4-14）與式（4-15）中，整理簡化，可得：（4-17）將式（4-17）帶入式式（4-15）中得：（4-18）（5）為求出第j號象元內(nèi)的平均修正值，假設(shè)j象元內(nèi)共有條射線通過，則（4-19）（6）用平均修正值對第j個象元的進(jìn)行修正。在修正時，值需受下列物理?xiàng)l件的約束，即若，則取若，則取和為介質(zhì)慢度的范圍，不同的介質(zhì)，慢度值的范圍不同。（7）對求得的值，用下式判斷其收斂程度（4-20）如果條件成立，則認(rèn)為值達(dá)到預(yù)定的收斂要求，否則還要做下一輪迭代，直至滿足條件時停止。從上面第（4）、（5）步計算過程可以看出不是用第i條射線的修正值來求慢度而是將所有射線得到的修改值保存下來，在本輪對射線迭代結(jié)束后求所有射線在象元內(nèi)的平均值，然后再由（4-21）對每個象元的波慢作修改，再判斷是否達(dá)到收斂要求。給定初始模型給定初始模型通過射線追蹤得到系數(shù)矩陣L是否滿足精度滿足精度輸出S是修改慢度向量S否圖4.3層析成像流程圖以一個簡單模型為例，將一個存在高速異常的介質(zhì)模型離散成的網(wǎng)格，第二行第二列和第三行第二列存在高速異常，運(yùn)用線性插值射線追蹤法得到射線路徑，再利用SIRT算法得到的結(jié)果如圖4.4所示。圖4.4網(wǎng)格模型及反演結(jié)果圖4.5層狀介質(zhì)層析成像結(jié)果圖圖4.5為一層狀介質(zhì)，在3-5m處為一低速層，20-24m處有一高速層，圖為迭代10次的結(jié)果圖，兩個速度異常的層區(qū)域的速度、位置、大小都以呈現(xiàn)出來，即層析成像的到了比較理想的結(jié)果。（a）存在高速異常的介質(zhì)（b）存在低速異常的介質(zhì)圖4.6存在方形異常介質(zhì)的層析成像結(jié)果圖4.6分別為存在低、高速方形異常的介質(zhì)，兩個速度異常的方形區(qū)域的速度位置、大小都呈現(xiàn)出來，即層析成像的到了比較理想的結(jié)果。（a）存在高速異常區(qū)的介質(zhì)（b）存在低速異常區(qū)的介質(zhì)圖4.7存在橫向異常介質(zhì)的層析成像結(jié)果圖4.7為從存在橫向異常介質(zhì)模型的到的層析成像結(jié)果，從圖中可以看出當(dāng)異常為橫向時無論速度值、異常區(qū)的位置、大小都達(dá)到較好的吻合。（a）存在低速異常區(qū)的介質(zhì)，介質(zhì)模型被劃分為的網(wǎng)格（b）存在低速異常區(qū)的介質(zhì)，介質(zhì)模型被劃分為的網(wǎng)格（c）存在高速異常區(qū)的介質(zhì)，介質(zhì)模型被劃分為的（d）存在低速異常區(qū)的介質(zhì)，介質(zhì)模型被劃分為的網(wǎng)格圖4.8存在垂向異常介質(zhì)的層析成像結(jié)果圖4.8中的四幅圖為存在垂向異常介質(zhì)模型的到的層析成像的結(jié)果。對比圖4.7、4.8可以看出當(dāng)異常區(qū)域?yàn)榇瓜驎r，結(jié)果不如異常為橫向時理想。

第五章結(jié)論本文從層析成像的基本理論出發(fā)，研究論述了層析成像正、反演兩個方面的算法，在正演方面，即射線追蹤方法，從算法的介紹和得到的圖像可以看出：基于旅行時線性插值射線追蹤法具有算法簡單、精度高等優(yōu)點(diǎn)，它適合于速度突變的情況，無須對單元附近的速度進(jìn)行平滑處理。因此這是一中很有應(yīng)用前途的射線追蹤算法。最短路徑射線追蹤法具有在復(fù)雜速度結(jié)構(gòu)下無條件穩(wěn)定的優(yōu)點(diǎn)，射線追蹤的精度主要取決于節(jié)點(diǎn)數(shù)M，精度要求越高，所需要的節(jié)點(diǎn)數(shù)越多。但此時大大增加了計算量。在反演方面，本文主要介紹了三種反演算法，分別是:ART算法、LSQR算法和SIRT算法，其中重點(diǎn)運(yùn)用了LSQR算法和SIRT算法，從算法的原理及公式的推導(dǎo)可以看出：LSQR法是利用Lanczos方法求解最小二乘問題的一種投影法，在迭代過程中，它只涉及到非零元素，減少了存貯空間，提高了運(yùn)算速度，所以特別適于求解系數(shù)為大型稀疏矩陣的方程組。SIRT算法利用一個像素內(nèi)通過的所有射線的修正值來確定這一像素的平均修正值。它是所有射線通過計算以后，才完成一次迭代。取平均修正值可以壓制一些干擾因素，而且SIRT的計算結(jié)果與觀測數(shù)據(jù)使用順序無關(guān)。從實(shí)例計算的結(jié)果看：上述模型在若干次迭代后，當(dāng)速度特征體為橫向時，其形狀、大小、位置與所給定的理論速度模型都能達(dá)到較好的吻合；但當(dāng)速度特征體為垂向時，得到的結(jié)果就不夠理想總體來說用跨孔雷達(dá)層析成像的方法進(jìn)行地下地質(zhì)地層結(jié)構(gòu)勘測是方便可行的，因?yàn)槔眠@種方法可以直觀的反應(yīng)地下介質(zhì)的某些參數(shù)，如：速度，慢度，以及介質(zhì)的吸收系數(shù)等等。從而，可以判斷地下異常體的位置，形狀，大小等。因此，此種方法是工程，環(huán)境，水文等領(lǐng)域中的一種十分有效的勘測方法。參考文獻(xiàn)宋文榮、吳儀芳、劉建達(dá)、許漢剛等，電磁波層析成像技術(shù)（CT）在地基勘探中的應(yīng)用[J]，地震學(xué)刊，1994年第二期，62-64。華衛(wèi)兵，基于代數(shù)重建算法的CT圖像仿真[D]，北京科技大學(xué)，2005。李才名，電磁波層析成像（EWCT）技術(shù)及其在混凝土樁基礎(chǔ)質(zhì)量檢測中的應(yīng)用[D]，南京大學(xué)，2003。張建中、陳世軍，復(fù)雜介質(zhì)地震初至波數(shù)值模擬[J]，計算物理，2003年第20卷第5期，429-433。李志輝，CT技術(shù)在路基病害勘查治理中的研究及應(yīng)用[D]，西南交通大學(xué)碩士論文，2006。張賽民，跨空地震層析成像研究[D]，中南大學(xué)碩士學(xué)位論文，2003。袁志亮，井間聲波電磁波層析成像技術(shù)應(yīng)用研究與軟件開發(fā)[D]，中國地質(zhì)大學(xué)博士學(xué)位論文，2007。楊文采，地球物理反演和地震層析成像[M]，地質(zhì)出版社，1989。赫爾曼，由投影重建圖像—CT的理論基礎(chǔ)[M]，科學(xué)出版社，1985。KiHaLeeandGanquanXie,Anewapproachtoimagingwithlow-frequencyelectromagneticfields[J],Geophyscs.VOL.58,No.6(June1993),780-796。S.HanafyandS.A.alHagrey,Ground-penetratingradartomographyforsoil-moistureheterogeneity,Geophysics[J],VOL.71,No.1(January-February2006),P.K9-K18。曾昭發(fā)、劉四新、王者江、薛建，探地雷達(dá)方法原理及應(yīng)用[M]，科學(xué)出版社，2006。黃家會、宋雷、崔廣心、楊維好，應(yīng)用跨孔雷達(dá)層析成像技術(shù)研究深部巖層特性[J]，中國礦業(yè)大學(xué)學(xué)報，1999年11月，第28卷第6期，578-581。王家映，地球物理反演理論[M]，高等教育出版社（第二版），2002。CoultripRL,High-accuracywavefronttracingtraveltimecalulation[J],Geophysics,58(2):284-292。JulianBR,GubbinsD,Three-dimensionalseismicraytracing[J],Geophy,1977,43:95-114。UmJ，ThurberC，Afastalgorithmfortwo-pointsseismicraytracing[J]，Bull.Seis.Soc.Am，1987,79:972-986。趙改善、郝守玲、楊爾皓等，基于旅行時線性插值的地震射線追蹤算法[J]，石油物探，1998，37（2）：14-24。EduardoL.FariaandPaulL.Stoffa,Traveltimecomputationintransverselyisotropicmedia[J],Geophysics,1994,59(12):272-281。Asawaka,E.andKawanaka,T.Seismicraytracingusinglineartraveltimeinterpolation[J],Geophysics,1993,57(2):326-333。張建中，復(fù)雜地區(qū)地震勘探靜校正方法技術(shù)及應(yīng)用研究[M]，廈門大學(xué)、勝利石油管理局博士后科研工作報告，2004。黃靚、黃政宇,線性插值射線追蹤的改進(jìn)方法[J]，湘潭大學(xué)自然科學(xué)學(xué)報，第24卷第4期，105-108。李志輝、劉爭平，最短路徑法射線追蹤的MATLAB實(shí)現(xiàn)[J]，工程地質(zhì)計算機(jī)應(yīng)用，2004，第4期，16-21。段心標(biāo)、金維浚，井間層析成像混合最短射線追蹤旅行時反演[J]，工程地球物理學(xué)報，2007年6月，第4卷第3期，201-206。王汪根、劉盛東、張平松，彎曲射線追蹤中的DIJKSTRA算法的改進(jìn)與實(shí)現(xiàn)[J]，地球物理學(xué)進(jìn)展，2006年12月，第21卷第4期，1120-1126。嚴(yán)蔚敏、吳偉民，數(shù)據(jù)結(jié)構(gòu)[M]，清華大學(xué)出版社，1997。項(xiàng)榮武、劉艷杰、胡忠盛，圖論中最短路徑問題的解法[J]，沈陽航空工業(yè)學(xué)院學(xué)報，2004年4月，第21卷第2期，86-88。張建中、陳世軍、余大祥，最短路徑射線追蹤方法及其改進(jìn)[J]，地球物理學(xué)進(jìn)展，2003年3月，第18卷第1期，146-150。張云姝、顧漢明、師學(xué)明，基于Moser曲線射線追蹤的SIRT聲波層析成像[J]，工程地球物理學(xué)報，2005年6月，第2卷第3期，167-176。陳國金、曹輝、吳永栓、唐金良，最短路徑層析成像技術(shù)在井間地震中的應(yīng)用[J]，石油物探，2004年7月，第43卷第4期，327-330。徐初偉、張建中、李海洋，基于薄層模型的有限頻率地震波射線追蹤[J]，勘探地球物理進(jìn)展，2004年6月，第27卷第3期，182-185。徐鳳生、李天志，所有最短路徑的求解算法[J]，計算機(jī)工程與科學(xué)，2006年第28卷第12期，83-84。李才明、張善法，電磁波層析成像阻尼因子引入與應(yīng)用[J]，地球物理學(xué)進(jìn)展，2005年3月，第20卷第1期，221-224。楊慧珠、聶建新、譚桂華，反射波地震層析成像算法研究[J]，遼寧工程技術(shù)大學(xué)學(xué)報，2004年4月，第23卷第2期，159-161。張秀軍、徐安農(nóng)、李安坤、蔣利華，改進(jìn)的共軛梯度法及其收斂性[J]，桂林電子工業(yè)學(xué)院學(xué)報，2005年12月，第25卷第6期，64-67。張順利、王建武、劉清，基于代數(shù)重建法的計算機(jī)斷層成像[J]，咸陽師范學(xué)院學(xué)報，2004年8月，第19卷第4期，34-37。杜正春、牛振勇、方萬良，基于分塊QR分解的一種狀態(tài)估計算法[J]，中國電機(jī)工程學(xué)報，2003年8月，第23卷第8期，50-55。王振宇、劉國華、梁國錢，基于廣義逆的層析成像反演方法研究[J]，浙江大學(xué)學(xué)報，2005年1月，第39卷第1期，1-5。何永斌、范嘯濤、安紅巖、何果，線性最小二乘問題解法的理論分析[J]，成都理工大學(xué)學(xué)報，2003年10月，第30卷第5期，529-533。宛新林、席道瑛、高爾根、沈兆武，用改進(jìn)的光滑約束最小二乘正交分解法實(shí)現(xiàn)電阻率三維反演[J]，地球物理學(xué)報，2005年3月，第48卷第2期，439-444。宿淑春、王守東、吳律，井間層析成像的平滑SIRT算法[J]，石油大學(xué)學(xué)報（自然科學(xué)版），2001年12月，第25卷第6期，29-31。趙連鋒，井間地震波速與衰減聯(lián)合層析成像方法研究[D]，成都理工大學(xué)博士論文2003年。

發(fā)表論文（1）楊薇，劉四新，馮彥謙，跨孔層析成像LSQR算法研究，物探與化探，2008年4月，第32卷第2期，199-202。

致謝首先,衷心感謝我的導(dǎo)師劉四新教授。兩年來，無論是在學(xué)習(xí)上還是在生活上，他都給了我無微不至的關(guān)懷和幫助。論文從選題到完成的每一個環(huán)節(jié)，無不滲透著他辛勤的汗水。他治學(xué)嚴(yán)謹(jǐn)、精益求精，他忘我工作、誨人不倦，他言傳身教、關(guān)愛有佳，深深地銘刻在我的心中，將永遠(yuǎn)激勵著我在科學(xué)的海洋里不懈拼博。感謝曾昭發(fā)教授，鹿琪副教授，王春暉，仝傳雪等師兄師姐和同窗好友的鼓勵和幫助。最后，真誠感謝本文所援引參考文獻(xiàn)的作者，是他們的成果給予了我知識和啟迪，讓我順利完成了碩士畢業(yè)論文。中文摘要作者姓名：楊薇專業(yè)：地球探測與信息技術(shù)導(dǎo)師姓名及職稱：劉四新教授探地雷達(dá)（GroundPenetratingRadar）是一種高效的淺層地球物理探測技術(shù)。它通過發(fā)射高頻電磁脈沖波，利用地下介質(zhì)電性參數(shù)的差異，根據(jù)回波的振幅、波形和頻率等運(yùn)動學(xué)和動力學(xué)特征來分析和推斷介質(zhì)結(jié)構(gòu)和物性特征。但探地雷達(dá)的探測范圍十分有限。鉆孔探地雷達(dá)是利用電磁波在地層中的傳播特性來獲取地層信息，從而解釋地下構(gòu)造的一種井中地球物理方法。它解決了地面探地雷達(dá)的探測范圍小的問題。跨孔層析成像模式是鉆孔探地雷達(dá)得主要工作模式的一種。這種技術(shù)是利用井間電磁波走時數(shù)據(jù)或振幅數(shù)據(jù)反演出研究區(qū)域的速度分布或衰減系數(shù)的分布的一種新的地球物理勘探方法。該技術(shù)主要依據(jù)電磁波在不同介質(zhì)中的速度不同或吸收系數(shù)不同進(jìn)行電磁波透射得到最大數(shù)據(jù)量，由計算機(jī)處理得到勘探剖面上的速度或吸收系數(shù)并以灰度或彩色成像，形象直觀地研究復(fù)雜地質(zhì)現(xiàn)象的成因及局部結(jié)構(gòu)，對于基巖中的斷層、軟弱夾層、地下溶洞及其發(fā)育特征、地下不明隱蔽物等勘探具有明顯的優(yōu)勢。本文從層析成像的基本理論出發(fā)，研究了層析成像的正、反演方法。本文在層析成像的正演方法方面主要論述了線性插值射線追蹤法和最短路徑法的原理與實(shí)現(xiàn)過程。線性插值法是一種由Asakawa提出基于旅行時的射線追蹤法，這種方法基于網(wǎng)格單元模型，適用于任意變速介質(zhì)，可以追蹤包括直達(dá)波、折射波、透射波等多種射線路徑。通過對假定的物理模型的計算可以看出，用該算法是一種快速、精確的計算旅行時及追蹤射線路徑的方法，該方法更加適應(yīng)速度變化大等復(fù)雜條件，這是一種很有發(fā)展前途的射線追蹤方法。最短路徑射線追蹤方法基于Fermat最小旅行時原理和網(wǎng)絡(luò)理論中的最短路徑算法來實(shí)現(xiàn)。最短路徑方法起源于網(wǎng)絡(luò)理論，首次由Nakanishi和Yamaguchi應(yīng)用于地震射線追蹤中。Moser以及Klimes和Kvasnicha對最短路徑方法進(jìn)行了詳細(xì)研究。通過科技人員的不斷研究，最短路徑方法目前已發(fā)展較為成熟，其基本算法的計算程序也較為固定。在層析成像反演方面，本文主要介紹了代數(shù)重建算法（ART算法）、最小二乘QR分解法（LSQR算法）、聯(lián)合迭代重建算法（SIRT算法）三種算法，推導(dǎo)出迭代公式，并在matlab平臺上編制程序，并利用給定模型檢驗(yàn)之。關(guān)鍵詞：層析成像，射線追蹤，線性插值法，最短路徑法，LSQR算法，SIRT算法

AbstractAutor:YangWeiMajor:GeoexplorationAndInformationTecknologyTutor’snameandtitle:LiuSixinProfessorGroundpenetratingradarisakindofefficientlytechniqueforshallowlayergeophysicalexplore.Itsendshighfrequencyelectromagnetismimpulsewaveandmakesuseofthedifferenceoftheelectricalpropertyofundergroundmediumtoanalyzeanddeducethestructureandcharacterofthemedium,basedonkinematicsanddynamicscharacterofthewave,suchasswing,waveformandfrequency.GPRisapplyingonmanyfieldssuchasengineering,hydrologyandenvironmentetc,andnowit’sapplicationscaleisenlarging.ThedetectionscaleoftheGPRislimitedbecauseofeffectoftectumandbeltofweathering.BoreholeantennasarevaluetoolstoprobethesubsurfacewhenthetopmostsoilshowhighconductivityandpreventpenetrationofEMenergyfromthesurface,orwhenthetargetisbeyondtherangeatagivenfrequency.Ithashighresolutionandlargedetectionscale.Ithastwoworkingmodes:oneissingle-holereflectionmode,anotheriscross-holetomographymode.Thistextmainlyintroducedthebasicprincipleandthecalculationmethodofcross-holetomographymode(engineeringCTtechnique),includingthemethodofgetthemodulusmatrixbycontrastthepointofradialandgridding,thecourseofsynthesizetimematrixwiththeslownessmatrixgiven,theprincipleofinversioncalculationwithLSQRcalculationmethod,andthecourseofconversiontheformofthematrixandmakingpicture.Inthispaperwestudyforwardmethodsandinversionmethodfromtheprincipleofcross-holetomography.Thispaperdiscussestheprinciplesandimplementationprocessofthelineartraveltimeinterpolationmethodandtheshortestpathsr