專利名稱：一種利用聲波探測物體內(nèi)部結(jié)構(gòu)的方法
技術(shù)領(lǐng)域：
本發(fā)明屬于超聲檢測技術(shù)領(lǐng)域，特別涉及一種利用聲波探測物體內(nèi)部結(jié)構(gòu)的方法。
利用發(fā)聲源使聲波在被探測物體中傳播，并通過聲接收器記錄下在被探測物體中傳播后的聲波信號，再經(jīng)過反演算方法確定出被探測物體的內(nèi)部結(jié)構(gòu)的方法現(xiàn)在已經(jīng)有了一些應(yīng)用。在這種探測方法中反演算的算法是整個方法的關(guān)鍵。
我們知道聲波在物體中的傳播一般是用彈性波動方程來描述的。而當切向波可以忽略或傳播距離較小時，聲信號可用一個標度波動方程ρ(x)κ(x)∂2Ps(x,t)∂t2=ρ(x)▿•[1ρ(x)▿Ps(x,t)]+δ(x-xs)f(t),----(1)]]>近似地描述。式中x為空間的位置矢量；t為時間；Ps(x，t)為由第s發(fā)聲源所產(chǎn)生的聲壓；xs為第s發(fā)聲源的位置矢量；f(t)為源所產(chǎn)生的脈沖波形；ρ(x)為物質(zhì)內(nèi)部的密度，κ(x)為物質(zhì)內(nèi)部的壓縮模量；式中表示對空間求散度，而·表示對空間求梯度。這里我們將近一步假定物質(zhì)內(nèi)部的密度ρ(x)＝ρ為常數(shù)。這樣方程(1)中將只有一個參量K(x)=κ(x)ρ=c(x)2]]>，其中c(x)為物體內(nèi)的壓縮聲波的傳播速度。如果我們通過運算能夠求出K(x)或c(x)，我們就可以在某種程度上確定被探測物體的內(nèi)部結(jié)構(gòu)。
目前在大部份的應(yīng)用中，如在探測地質(zhì)結(jié)構(gòu)的應(yīng)用中所用的反演算的算法主要是基于單向波場理論，正如C.P.A.Wapenaar和A.J.Berkhout，《彈性波場的演化》(Elastic wave fi eld extrapolation)，(Elsevier，Amsterdam，1989)所指出的，這一理論首先將聲場分解成正向傳播 和負向傳播 的兩部分，并將它們寫成兩個相互有偶合的一階偏微分方程組。這個有相互偶合的一階偏微分方程組與方程(1)是等價的。但在實際處理中，為了使實際運算得以進行，在計算中需要忽略這個相互偶合項。忽略這個相互偶合項在物理上意味著忽略了多次反射效應(yīng)。而這一點從根本上限制了基于單向波場理論的反演算方法的發(fā)展，使這種方法難以探測較多的巖層。這是因為經(jīng)過多層后傳出的聲波信號總是與多次反射的信號混在一起的。這種方法可探測的標準層數(shù)一般為3-5層。
為了能夠探測出更多的層數(shù)，嚴格按照聲波波動方程(1)進行的全反演算的算法有著更光明的前景。嚴格按照聲波波動方程(1)進行的全反演算的基本手段之一是非線性最小二乘法。它的計算目的是使平方和Σr,s,t[δds(xr,t)]2=Σr,s,t[Ps(xr,t)-Pms(xr,t)]2----(2)]]>最小。其中δds(xr,t)≡Ps(xr,t)-Pms(xr,t);Σr,s,t]]>表示對所有發(fā)聲源，接受器及時間求和，xr為第r接收器的位置矢量，Ps(xr，t)表示在假定模型為K(x)的條件下按(1)式計算出的壓力，而 表示實際測量的壓力。
非線性最小二乘法與線性最小二乘法的最大區(qū)別為線性最小二乘法可一次求得滿足公式(2)的條件，這里即K(x)的值，且結(jié)果是穩(wěn)定的；而非線性最小二乘法則需要在某一初始值進行多次的迭代，且最終并不一定能得到正確的結(jié)果。在我們的問題中假定初始模型為K(x)，非線性最小二乘法用來修改K(x)的公式為ATA(δK)＝AT(δd)， (3)這里A是偏微商矩陣，階數(shù)為L×Nx。其中Nx為描寫K(x)所需的參數(shù)的個數(shù)，如在完整的有限元空間描寫K(x)所需的參數(shù)為將K(x)分立化時所產(chǎn)生的點數(shù)；而L＝nr×ns×nt，其中nr和ns分別為在反演算中所用的接收器和發(fā)聲源的個數(shù)，nt為在反演算中所使用的分立的時間的步數(shù)。AT為矩陣A的轉(zhuǎn)置。δd為L階矢量，它是公式(2)中的δds(xr，t)的矢量描述形式。δK＝δK(x)為一個Nx階的矢量，它是模型參數(shù)K(x)的修正量。這里K(x)也為Nx階矢量，它是模型參數(shù)的矢量描述形式。修正后的模型參數(shù)為K(x)+δK(x)。在作完此修正后，將把修正后的模型參數(shù)做為新的初始模型K(x)，并重復(fù)上述修正過程。
在實際的反演算運算過程中，最小二乘法的修正式(3)對于δK經(jīng)常是病態(tài)的，即解不唯一。因而，在實際計算中需將原始的最小二乘法作一些改進使它能有效的處理病態(tài)問題。阻尼最小二乘法是這類改進中的一種，應(yīng)用它可在一定程度上改善方程的病態(tài)性質(zhì)。阻尼最小二乘法的主要原理是在原始最小二乘法的基礎(chǔ)上加入一個能使模型K(x)在一定程度上保持原有數(shù)值的衰減系數(shù)γ。在阻尼最小二乘法中，改進模型K(x)的公式為(ATA+γI)(δK)＝AT(δd). (4)其中，I為單位矩陣。在計算中我們需適當調(diào)整衰減系數(shù)γ的數(shù)值，使每次模型K(x)的改進向好的方向進行。
利用這種非線性阻尼最小二乘法，在某一初始結(jié)構(gòu)下(如將內(nèi)部結(jié)構(gòu)設(shè)置成與表面相等的均勻值)，經(jīng)過M0次的迭代，逐漸改變模型，最終得到與被探測物體的內(nèi)部性質(zhì)一致的K(x)的完整流程圖，如圖1所示，其中，K(x)＝K(x)+δK(x)代表阻尼最小二乘法。但是，在實際的全反演算的模擬計算中，這種簡單的應(yīng)用最小二乘法的反演算算法沒能達到人們所期望的效果，如F.Jurado，M.Cuer和V.Richard，在《地球物理》(Geophysics)，60，1857(1995)；W.W.Symes和Carazzone，在《地球物理》(Geophysics)，56，654(1991)中分析失敗原因所指出的，人們認為主要是由于在反演算的過程中，方程(3)或(4)通常具有病態(tài)的性質(zhì)，以及殘差的平方和，即(2)式具有眾多的局域極小值。這一系列問題嚴重地妨害了反演算的正確進程。而如何解決這些問題至今還沒有找到有效的答案。為此，利用全反演算的算法來探測物體內(nèi)部結(jié)構(gòu)，至今還處于探索階段。
本發(fā)明的目的在于提供一種利用聲波探測物體內(nèi)部結(jié)構(gòu)的方法，其是利用發(fā)聲源在被探測物體表面發(fā)聲，使聲波在被探測物體中傳播，并在被探測物體的同一表面用聲接收器接收從物體內(nèi)部傳出的聲信號，將發(fā)聲源所發(fā)的聲波的特性，及聲接收器接收到的在被探測物體中傳播后的聲信號，通過基于聲波波動方程進行的全反演算的算法將被探測物體的內(nèi)部聲學參數(shù)(這里為物體內(nèi)的壓縮波的聲速或等價的參數(shù))計算出來。在計算中加入了與傳統(tǒng)修正方式不同的步驟，有效地防止了(2)式的殘差平方和落入局域極小值，從而可準確地確定出被探測物體的內(nèi)部結(jié)構(gòu)。
本發(fā)明利用發(fā)聲源、聲接收器、以及通過聲波波動方程的全反演算法對物體內(nèi)部結(jié)構(gòu)進行探測，可用來對地質(zhì)結(jié)構(gòu)、海洋中的聲速水平分布、大氣分層等進行探測。
本發(fā)明包括計算機、發(fā)聲源、聲接收器和被探測物體，其特征在于(1).將聲接收器均勻分布放置在被探測物體的某一表面上，記錄下它們的位置；在放置聲接收器的同一表面上選定擬放發(fā)聲源的位置，并記錄下它們的位置；將在每一處接收到的由發(fā)聲源所發(fā)出的聲壓存入計算機，當對所有擬放發(fā)聲源的位置進行發(fā)聲和接收記錄后，結(jié)束數(shù)據(jù)采集工作；(2).在基于阻尼最小二乘法，利用聲波波動方程的全反演算法的基礎(chǔ)上，通過加入與標準阻尼最小二乘法不同的階梯式修正法和結(jié)構(gòu)重整步驟，且每次的階梯式修正法，阻尼最小二乘法和結(jié)構(gòu)重整過程都對當前模型進行修正，并將每次修正后所得的模型參數(shù)K(x)設(shè)置為新的當前模型，經(jīng)多次重復(fù)整體迭代，使Σr,s,t[Ps(xr,t)-Pms(xr,t)]2]]>的值更有效地趨于整體的最小值，從而確定被探測物體的內(nèi)部聲學參數(shù)。
所述的整體迭代，首先將聲學參數(shù)K(x)的初始結(jié)構(gòu)設(shè)置為內(nèi)部結(jié)構(gòu)與表面值相等的均勻值，整體迭代的指標J設(shè)置為1；初始化完成后，整體迭代部分首先設(shè)置對階梯式修正法的迭代的指標J1為0，這時進入一個J1是否小于M1的判斷，當J1確實小于M1時，進行階梯式修正法的修正；修正完畢后，將迭代的指標J1在現(xiàn)有值上加1，并再次進入J1是否小于M1的判斷；如果J1小于M1，重復(fù)進行上面的過程直至J1等于M1；完成了整體迭代中的第一部分后，則進入整體迭代中的第二部分，在此設(shè)置對阻尼最小二乘法的迭代指標J0為0，然后進入一個對J0是否小于M0的判斷，當J0確實小于M0時，進行阻尼最小二乘法的修正；修正完畢后，將迭代指標J0在現(xiàn)有值上加1，并再次進入J0是否小于M0的判斷；如果J0小于M0，重復(fù)進行上面的過程直至J0等于M0；完成了整體迭代中的前兩部分后，將進入對整體迭代的指標J是否小于M的判斷，當J確實小于M時，進行結(jié)構(gòu)重整；結(jié)構(gòu)重整后，將整體迭代指標J在現(xiàn)有值上加1，并回到下一次的整體迭代；多次重復(fù)整體迭代直至J等于M，當J＝M后，將進行最后一次對整體迭代中的前兩部分的運算，而不再作結(jié)構(gòu)重整，然后將當前模型參數(shù)K(x)的結(jié)果輸出，并同時將Σr,s,t[Ps(xr,t)-Pms(xr,t)]2]]>的值輸出，用以檢驗最后結(jié)果的正確性，結(jié)束利用聲波探測物體內(nèi)部結(jié)構(gòu)的方法。
所述的階梯式修正法是通過對由阻尼最小二乘法所求得的修正量δK(x)進行一個相應(yīng)的變換來求得當前模型的修正量δK′(x)，其首先確定由通常的阻尼最小二乘法所求得的修正量δK(x)的極值的位置及極值的性質(zhì)，當兩極值的距離足夠接近時，確定它們的中點為起始修正點，它們的深層的極值和淺層的極值的差值為該起始修正點之后的修正值；找出滿足上述條件的所有修正值中的最大值δK′max，設(shè)定一個小于1的相對域值τ，對所有修正值大于τδK′max的地方進行修正，其為深于起始修正點的所有聲學參數(shù)值加上它相應(yīng)的修正值，在每一次計算中可同時修正復(fù)數(shù)個這樣的起始修正點和修正值，這些修正的總和構(gòu)成δK′(x)。
所述的結(jié)構(gòu)重整是首先將現(xiàn)有的模型參數(shù)做調(diào)整，設(shè)定一個小的域值和一個新模型參數(shù)，其新模型參數(shù)的表面值與現(xiàn)有模型參數(shù)的表面值相等，從表面開始逐次做參數(shù)重整，直至整個新模型參數(shù)被確定；在某一當前位置，當現(xiàn)有模型參數(shù)值與淺層鄰近的現(xiàn)有模型參數(shù)值之差的絕對值小于或等于這個域值時，則設(shè)定這一位置的新模型參數(shù)值與它淺層的新模型參數(shù)值相等，而當這個差值的絕對值大于這個域值時，則設(shè)定這一位置的新模型參數(shù)值為淺層鄰近的新模型參數(shù)值加上這個差值；最后，將所得的新模型參數(shù)值設(shè)置為新的當前模型以進行新的迭代運算。
本發(fā)明由4部分組成1.計算機；2.發(fā)聲源；3.聲接收器；4.被探測物體。其中，每一部分具有如下特征。
計算機有一定計算能力的計算機。
發(fā)聲源可發(fā)出某種特殊脈沖波形的發(fā)聲源。如電聲換能器、爆炸等。
聲接收器可在某特定點實時記錄該點的聲壓的聲接收器。如電聲換能器、聲強計、實時聲強分析儀等。
被探測物體被探測物體由具有不同聲學參數(shù)的物質(zhì)構(gòu)成。如由具有不同聲速的物質(zhì)構(gòu)成的物體。(注本發(fā)明是根據(jù)聲學參數(shù)的分布來確定物質(zhì)內(nèi)部結(jié)構(gòu)的，相同聲學參數(shù)的不同物質(zhì)在本發(fā)明中是無法區(qū)分的)。
本發(fā)明的工作步驟如下將聲接收器盡可能均勻分布地放置在被探測物體的某一表面上，并記錄下它們的位置，做好接收聲信號的準備。在放置聲接收器的同一表面上選定擬放發(fā)聲源的位置，并記錄下它們的位置。擬放發(fā)聲源的位置也應(yīng)盡可能均勻地分布在這一表面。擬放發(fā)聲源的位置、放置聲接收器的位置和被探測物體的關(guān)系如圖2所示。
逐一在擬放發(fā)聲源的位置xs放置發(fā)聲源，并使其發(fā)放一脈沖聲波f(t)。從發(fā)聲源開始發(fā)聲時起直至聲波在物體內(nèi)部消失的時段里，記錄所有聲接收器接收到的實時聲壓波形。將在每一xr處接收到的由第s發(fā)聲源所發(fā)出的聲壓記為 ，并存入計算機。當對所有擬放發(fā)聲源的位置進行過前述的發(fā)聲和接收記錄后，數(shù)據(jù)的采集工作完結(jié)。
在數(shù)據(jù)采集完結(jié)后，本發(fā)明將由計算機根據(jù)采集到的數(shù)據(jù)，通過聲波波動方程的全反演算法將被探測物體的內(nèi)部聲學參數(shù)計算出來。我們所采用的算法是在實時空間的阻尼最小二乘法(圖1)的基礎(chǔ)上，通過加入一個與標準阻尼最小二乘法不同的稱做階梯式修正法的算法，和稱做結(jié)構(gòu)重整的算法，使公式(2)的殘差平方和更有效地趨于整體的最小值。具體計算過程如圖3所示。
圖3中的初始結(jié)構(gòu)的作用為將當前模型的聲學參數(shù)K(x)設(shè)置為內(nèi)部結(jié)構(gòu)與表面值相等的均勻值。
整個計算過程包括M次的整體迭代(對圖3中的指標J)。在每次整體迭代又包括M1次的對階梯式修正法的迭代(對圖3中的指標J1)，M0次的對阻尼最小二乘法的迭代(對圖3的中指標J0)和一次結(jié)構(gòu)重整。每次的階梯式修正法，阻尼最小二乘法和結(jié)構(gòu)重整過程都按照下面介紹的特定方式對當前模型進行修正，并將修正后所得的模型參數(shù)K(x)設(shè)置為新的當前模型。
圖3中的阻尼最小二乘法對應(yīng)于通常的阻尼最小二乘法。它的聲學參數(shù)修正量δK(x)可由通常的阻尼最小二乘法線性方程組，即公式(4)求得。
圖3中的階梯式修正法是本發(fā)明所特有的修正法。它對當前模型的修正量δK＇(x)是通過對由通常的阻尼最小二乘法所求得的修正量δK(x)進行一個相應(yīng)的變換來求得的。具體變換的規(guī)則如下。首先確定由通常的阻尼最小二乘法所求得的修正量δK(x)的極值的位置及極值的性質(zhì)(極大值或極小值)。當兩極值的距離足夠接近時(如小于等于3)，我們確定它們的中點為起始修正點，它們的差值(深層的極值減去淺層的極值)為該起始修正點之后的修正值。找出滿足上述條件的所有修正值中的最大值δK＇max。設(shè)定一個相對域值τ(如等于0.8)。對所有修正值大于τδK′max的地方做如下修正，即深于起始修正點的所有聲學參數(shù)值將加上它相應(yīng)的修正值。在每一次計算中可同時修正復(fù)數(shù)個這樣的起始修正點和修正值，這些修正的總和構(gòu)成δK＇(x)。
圖3中的結(jié)構(gòu)重整將現(xiàn)有的模型參數(shù)做如下的調(diào)整。設(shè)定一個小的域值(如等于0.01)和一個新模型參數(shù)。新模型參數(shù)的表面值與現(xiàn)有模型參數(shù)的表面值相等。從表面開始逐次做下述重整過程直至整個新模型參數(shù)被確定。在某一當前位置，當現(xiàn)有模型參數(shù)值與淺層鄰近的現(xiàn)有模型參數(shù)值之差的絕對值小于或等于這個域值時、則設(shè)定這一位置的新模型參數(shù)值與它淺層的新模型參數(shù)值相等，而當這個差值(當前位置的現(xiàn)有模型參數(shù)值減去淺層鄰近的現(xiàn)有模型參數(shù)值)的絕對值大于這個域值時，則設(shè)定這一位置的新模型參數(shù)值為淺層鄰近的新模型參數(shù)值加上這個差值。最后，將所得的新模型參數(shù)值設(shè)置為新的當前模型以進行新的迭代運算。
在進行階梯式修正法和阻尼最小二乘法的修正中都將用到阻尼最小二乘法的衰減系數(shù)γ，它的值在運算中被調(diào)整至使修正后的殘差平方和，即(2)式比修正前的殘差平方和要小或相等。由于當衰減系數(shù)趨于無窮大時意味著沒有修正，這時殘差的平方和不變。因此，逐漸加大衰減系數(shù)使殘差平方和變小或不變顯然總是可以做到的。
在實際計算中可根據(jù)計算機的計算能力確定一組合適的M、M1和M0。典型的M1和M0的值約為10左右。M的值越大結(jié)果越準確可靠。
按照圖3的流程圖以及前面的說明，模型參數(shù)將按下述的方式逐漸逼近物體所具有的真實值。首先、將聲學參數(shù)K(x)的初始結(jié)構(gòu)設(shè)置為內(nèi)部結(jié)構(gòu)與表面值相等的均勻值(如圖4)。同時將整體迭代次數(shù)、階梯式修正法的迭代次數(shù)和阻尼最小二乘法的迭代次數(shù)分別確定為M、M1和M0，并將整體迭代的指標J設(shè)置為1。初始化完成后，計算進入整體迭代部分。在整體迭代部分首先設(shè)置對階梯式修正法的迭代的指標J1為0。這時進入一個J1是否小于M1的判斷。當J1確實小于M1時，進行階梯式修正法的修正。我們稱這個修正法為階梯式修正法是因為這個修正法每次給出的修正值均為一種階梯形式，這一點在從最早的初始模型(即與表面值相等的均勻值)進行第一次的階梯式修正時最為明顯(如圖5)。從圖5中可以看出修正后的模型參數(shù)從原來的均勻值變?yōu)橐粋€階梯形式。修正完畢后，將迭代的指標J1在現(xiàn)有值上加1，并再次進入J1是否小于M1的判斷。如果J1小于M1，重復(fù)進行上面的過程直至J1等于M1。當J1等于M1時，我們應(yīng)已經(jīng)進行了M1次的階梯式修正法的修正。這時模型參數(shù)已為M1次的階梯式修正的疊加，因此模型參數(shù)的階梯形式已不如單一次的修正時(圖5)明顯。但淺層的大模樣已接近物體所具有的真實值，但在淺層仍存在細小的差別(圖6)。完成了整體迭代中的第一部分后，我們將進入整體迭代中的第二部分。在此我們設(shè)置對阻尼最小二乘法的迭代指標J0為0。然后進入一個對J0是否小于M0的判斷。當J0確實小于M0時，進行阻尼最小二乘法的修正。阻尼最小二乘法的修正的作用是在前面的計算基礎(chǔ)上使淺層仍存在的細小差別逐漸消失。修正完畢后，將迭代指標J0在現(xiàn)有值上加1，并再次進入J0是否小于M0的判斷。如果J0小于M0，重復(fù)進行上面的過程直至J0等于M0。當J0等于M0時，我們應(yīng)已經(jīng)進行了M0次的阻尼最小二乘法的修正。這時模型參數(shù)在淺層的細小差別已明顯減少，而深層則變?yōu)橐环N連續(xù)變化與階梯變化的混合形式(圖7)。從圖7中可以看到在深層的階梯變化形式與物體所具有的真實值相似，而連續(xù)變化部分與真實值不符。完成了整體迭代中的前兩部分后，我們將進入對整體迭代的指標J是否小于M的判斷。當J確實小于M時，進行結(jié)構(gòu)重整。結(jié)構(gòu)重整的作用是將與真實值不符的連續(xù)變化部分濾掉，使模型參數(shù)仍為一種階梯變化的形式(圖8)。結(jié)構(gòu)重整后，將整體迭代指標J在現(xiàn)有值上加1，并回到下一次的整體迭代。每次整體迭代會使更深層的模型參數(shù)與物體所具有的真實值一致。圖9顯示了經(jīng)過4次整體迭代后的模型參數(shù)。該圖顯示深層的模型參數(shù)已經(jīng)與真實值很接近，它比經(jīng)過1次整體迭代后的模型參數(shù)(圖8)改進了許多。多次重復(fù)整體迭代直至J等于M。當J＝M后，我們將進行最后一次對整體迭代中的前兩部分的運算(即對階梯式修正法的迭代和對阻尼最小二乘法的迭代)，而不再作結(jié)構(gòu)重整，然后將當前模型參數(shù)做為最后結(jié)果輸出。圖10為模型參數(shù)的最后結(jié)果的示意圖。這個結(jié)果與物體所具有的真實值完全一致。
本發(fā)明的最大特點是在全反演算的算法中加入了不同于現(xiàn)有最小二乘法的階梯式修正法和結(jié)構(gòu)重整。使基于聲波波動方程進行的全反演算的算法有了實際的效果。應(yīng)用它可比現(xiàn)有的基于單向或雙向波場理論的反演算的算法探測出更復(fù)雜和更深層的內(nèi)部結(jié)構(gòu)。本發(fā)明的方法在接收到的聲波信號有一定誤差時仍能得到穩(wěn)定的正確結(jié)果，這是其它的利用全反演算的算法探測物體內(nèi)部結(jié)構(gòu)的方法所作不到的。
下面結(jié)合附圖和實施例對本發(fā)明的技術(shù)方案作進一步的描述。


圖1.基于簡單最小二程法的算法流程圖，其中J0為循環(huán)指標，M0為循環(huán)次數(shù)；K(x)＝K(x)+δK(x)為阻尼最小二乘法；圖2.本發(fā)明所對應(yīng)的被探測物體K0、發(fā)聲源●和聲接收器○的相對位置示意圖；圖3.本發(fā)明的算法流程圖，其中J、J1和J0為循環(huán)指標，M、M1和M0為循環(huán)次數(shù)；K(x)＝K(x)+δK(x)為階梯式修正法，K(x)＝K(x)+δK(x)為阻尼最小二乘法；圖4.最早的初始結(jié)構(gòu)的示意圖，圖中實線為初始結(jié)構(gòu)，點線為模型的真實值；圖5.從最早的初始模型進行第一次階梯式修正后的計算結(jié)果，圖中實線為計算結(jié)果，點線為模型的真實值；圖6.從最早的初始模型進行M1次的階梯式修正后的計算結(jié)果，圖中實線為計算結(jié)果，點線為模型的真實值；圖7.在圖6的基礎(chǔ)上進行M0次的阻尼最小二乘法的修正后的計算結(jié)果，圖中實線為計算結(jié)果，點線為模型的真實值；圖8.從最早的初始模型完成一次整體迭代后的計算結(jié)果，圖中實線為計算結(jié)果，點線為模型的真實值；圖9.從最早的初始模型完成4次整體迭代后的計算結(jié)果，圖中實線為計算結(jié)果，點線為模型的真實值；圖10.最終的計算結(jié)果，圖中實線為計算結(jié)果，點線為模型的真實值；圖11.實施例1中的被探測柱狀物體的截面圖，Km為膠凍物體，其尺寸為40厘米×40厘米，Am為吸收層。
實施例1將本發(fā)明的方法具體應(yīng)用到探測一個柱狀膠凍物體，它的截面圖如圖11所示。對于一個柱狀物體，如果它的長度遠大于它的寬度，則在柱中心附近的截面內(nèi)其波動方程可用二維波動方程近似地描述。對于膠凍物體切向波可忽略，因而聲波僅為壓縮波。為了簡化問題，被探測物體由密度近似相等的膠凍物體制作，且壓縮模量僅隨深度變化。這樣本被探測物體只有一個聲學參數(shù)K(x)，且它僅為深度的函數(shù)。適當選取時間單位使表面的模型參數(shù)K(x)為1.0。在這一時間單位下，準備一個可產(chǎn)生f(t)=(π2(t-10)282-12)exp(-π2(t-10)282)]]>的脈沖聲波的電聲換能器作為發(fā)聲源。選擇這樣的脈沖聲波的理由是這種脈沖形式與爆炸所產(chǎn)生的脈沖形式類似，因而使本實施例的結(jié)果可類比用爆炸來探測地質(zhì)結(jié)構(gòu)的情況。另外準備10個電聲換能器作為聲接收器，并將其均勻放置在被探測膠凍物體的中心橫截面的表面上，其具體位置為從膠凍物體的中心橫截面的表面的左端開始4cm、8cm、12cm、16cm、20cm、24cm、28cm、32cm、36cm和40cm處。類似地、在膠凍物體的中心橫截面的表面上均勻選取擬放發(fā)聲源的位置，其具體位置為從膠凍物體的中心橫截面的表面的左端開始1cm、5cm、9cm、13cm、17cm、21cm、25cm、29cm、33cm和37cm處。逐一在擬放發(fā)聲源的位置放置作為發(fā)聲源的電聲換能器，并使其產(chǎn)生脈沖為f(t)的聲波。從t＝0時刻開始至t＝512的時段里(t＝512時膠凍物體內(nèi)部的聲波已基本消失)，用作為聲接收器的各個電聲換能器接收其所在位置的實時聲壓波形。將在xr處的電聲換能器接收到的由第s發(fā)聲源所發(fā)出的聲壓波形記為 ，并存入計算機。當對所有擬放發(fā)聲源的位置進行過前述的發(fā)聲和接收記錄后，數(shù)據(jù)的采集工作完結(jié)。電聲換能器的相對誤差小于5％。
在數(shù)據(jù)采集完結(jié)后，將由計算機根據(jù)圖3中全反演算的算法，利用現(xiàn)已采集到的數(shù)據(jù)將被探測物體的內(nèi)部聲學參數(shù)計算出來。選取圖3中的M＝8、M1＝M0＝10。初始結(jié)構(gòu)設(shè)置為內(nèi)部結(jié)構(gòu)與表面值相等的均勻值，在此即K(x)＝1.0(如圖4)。初始模型確定后將整體迭代的指標J設(shè)置為1，并進入整體迭代部分。在整體迭代部分首先設(shè)置對階梯式修正法的迭代指標J1為0。這時進入一個J1是否小于M1的判斷。當J1確實小于M1時，進行階梯式修正法的修正。從最早的初始模型經(jīng)過第一次的階梯式修正后，階梯式修正法給出了一個階梯形式模型參數(shù)(如圖5)。修正完畢后，將迭代指標J1在現(xiàn)有值上加1，并再次進入J1是否小于M1的判斷。如果J1小于M1，重復(fù)進行上面的過程直至J1等于M1。當J1等于M1時，我們已經(jīng)進行了M1次的階梯式修正法的修正。這時模型參數(shù)已為M1次的階梯式修正的疊加，因此模型參數(shù)的階梯形式已不如單一次的修正時(圖5)明顯。此時在淺層，模型參數(shù)的大模樣已接近物體所具有的真實值，但仍存在細小的差別，而在深層則存在較大的差別(圖6)。完成了整體迭代中的第一部分后，我們進入整體迭代中的第二部分。在此我們設(shè)置對阻尼最小二乘法的迭代指標J0為0。然后進入一個對J0是否小于M0的判斷。當J0確實小于M0時，進行阻尼最小二乘法的修正。修正完畢后，將迭代指標J0在現(xiàn)有值上加1，并再次進入J0是否小于M0的判斷。如果J0小于M0，重復(fù)進行上面的過程直至J0等于M0。當J0等于M0時，我們已經(jīng)進行了M0次的阻尼最小二乘法的修正。這時模型參數(shù)在淺層的細小差別已明顯減少，而深層則變?yōu)橐环N連續(xù)變化與階梯變化的混合形式(圖7)。完成了整體迭代中的前兩部分后，我們進入對整體迭代的指標J是否小于M的判斷。當J確實小于M時，進行結(jié)構(gòu)重整。結(jié)構(gòu)重整后，與真實值不符的連續(xù)變化部分被部分地濾掉，使模型參數(shù)仍變回到一種階梯變化的形式(圖8)。完成結(jié)構(gòu)重整后，將整體迭代指標J在現(xiàn)有值上加1，并回到下一次的整體迭代。每次整體迭代會使更深層的模型參數(shù)與物體所具有的真實值更趨于一致。經(jīng)過4次整體迭代后的模型參數(shù)如圖9所示。該圖顯示深層的模型參數(shù)已經(jīng)與真實值很接近，它比經(jīng)過1次整體迭代后的模型參數(shù)(圖8)改進了許多。多次重復(fù)整體迭代直至J等于M。當J＝M后，即J＝8，我們進行最后一次對整體迭代中的前兩部分的運算(即對階梯式修正法的迭代和對阻尼最小二乘法的迭代)，而不再作結(jié)構(gòu)重整，然后將當前模型參數(shù)做為最后結(jié)果輸出。圖10顯示了本實施例的最終結(jié)果。從圖10可以看到經(jīng)8次整體迭代后，實線和點線完全吻合，表明8層不同聲學參數(shù)的區(qū)域均被準確探測到。
在本實施例中，膠凍物體共有8層不同參數(shù)的區(qū)域，大大超出了應(yīng)用基于單向波場理論所能探測的標準層數(shù)(3-5層)。此外、由于電聲換能器發(fā)出和接收到的聲波信號有一定誤差，而本發(fā)明的方法在此情況下仍能得到穩(wěn)定的正確結(jié)果。這是其它的利用全反演算的算法探測物體內(nèi)部結(jié)構(gòu)的方法所作不到的。
權(quán)利要求
1.一種利用聲波探測物體內(nèi)部結(jié)構(gòu)的方法，包括計算機、發(fā)聲源、聲接收器和被探測物體，其特征在于(1).將聲接收器均勻分布放置在被探測物體的某一表面上，記錄下它們的位置；在放置聲接收器的同一表面上選定擬放發(fā)聲源的位置，并記錄下它們的位置；將在每一處接收到的由發(fā)聲源所發(fā)出的聲壓存入計算機，當對所有擬放發(fā)聲源的位置進行發(fā)聲和接收記錄后，結(jié)束數(shù)據(jù)采集工作；(2).在基于阻尼最小二乘法，利用聲波波動方程的全反演算法的基礎(chǔ)上，通過加入與標準阻尼最小二乘法不同的階梯式修正法和結(jié)構(gòu)重整步驟，且每次的階梯式修正法，阻尼最小二乘法和結(jié)構(gòu)重整過程都對當前模型進行修正，并將每次修正后所得的模型參數(shù)K(x)設(shè)置為新的當前模型，經(jīng)多次重復(fù)整體迭代，使Σr,s,t[Ps(xr,t)-Pms(xr,t)]2]]>的值更有效地趨于整體的最小值，從而確定被探測物體的內(nèi)部聲學參數(shù)。
2.如權(quán)利要求1所述的利用聲波探測物體內(nèi)部結(jié)構(gòu)的方法，其特征在于所述的整體迭代，首先將聲學參數(shù)K(x)的初始結(jié)構(gòu)設(shè)置為內(nèi)部結(jié)構(gòu)與表面值相等的均勻值，整體迭代的指標J設(shè)置為1；初始化完成后，整體迭代部分首先設(shè)置對階梯式修正法的迭代的指標J1為0，這時進入一個J1是否小于M1的判斷，當J1確實小于M1時，進行階梯式修正法的修正；修正完畢后，將迭代的指標J1在現(xiàn)有值上加1，并再次進入J1是否小于M1的判斷；如果J1小于M1，重復(fù)進行上面的過程直至J1等于M1；完成了整體迭代中的第一部分后，則進入整體迭代中的第二部分，在此設(shè)置對阻尼最小二乘法的迭代指標J0為0，然后進入一個對J0是否小于M0的判斷，當J0確實小于M0時，進行阻尼最小二乘法的修正；修正完畢后，將迭代指標J0在現(xiàn)有值上加1，并再次進入J0是否小于M0的判斷；如果J0小于M0，重復(fù)進行上面的過程直至J0等于M0；完成了整體迭代中的前兩部分后，將進入對整體迭代的指標J是否小于M的判斷，當J確實小于M時，進行結(jié)構(gòu)重整；結(jié)構(gòu)重整后，將整體迭代指標J在現(xiàn)有值上加1，并回到下一次的整體迭代；多次重復(fù)整體迭代直至J等于M，當J＝M后，將進行最后一次對整體迭代中的前兩部分的運算，而不再作結(jié)構(gòu)重整，然后將當前模型參數(shù)K(x)的結(jié)果輸出，并同時將Σr,s,t[Ps(xr,t)-Pms(xr,t)]2]]>的值輸出，用以檢驗最后結(jié)果的正確性，結(jié)束利用聲波探測物體內(nèi)部結(jié)構(gòu)的方法。
3.如權(quán)利要求1或2所述的利用聲波探測物體內(nèi)部結(jié)構(gòu)的方法，其特征在于所述的階梯式修正法是通過對由阻尼最小二乘法所求得的修正量δK(x)進行一個相應(yīng)的變換來求得當前模型的修正量δK′(x)，其首先確定由通常的阻尼最小二乘法所求得的修正量δK(x)的極值的位置及極值的性質(zhì)，當兩極值的距離足夠接近時，確定它們的中點為起始修正點，它們的深層的極值和淺層的極值的差值為該起始修正點之后的修正值；找出滿足上述條件的所有修正值中的最大值δK′max，設(shè)定一個小于1的相對域值τ，對所有修正值大于τδK＇max的地方進行修正，其為深于起始修正點的所有聲學參數(shù)值加上它相應(yīng)的修正值，在每一次計算中可同時修正復(fù)數(shù)個這樣的起始修正點和修正值，這些修正的總和構(gòu)成δK＇(x)。
4.如權(quán)利要求1或2所述的利用聲波探測物體內(nèi)部結(jié)構(gòu)的方法，其特征在于所述的結(jié)構(gòu)重整是首先將現(xiàn)有的模型參數(shù)做調(diào)整，設(shè)定一個小的域值和一個新模型參數(shù)，其新模型參數(shù)的表面值與現(xiàn)有模型參數(shù)的表面值相等，從表面開始逐次做參數(shù)重整，直至整個新模型參數(shù)被確定；在某一當前位置，當現(xiàn)有模型參數(shù)值與淺層鄰近的現(xiàn)有模型參數(shù)值之差的絕對值小于或等于這個域值時，則設(shè)定這一位置的新模型參數(shù)值與它淺層的新模型參數(shù)值相等，而當這個差值的絕對值大于這個域值時，則設(shè)定這一位置的新模型參數(shù)值為淺層鄰近的新模型參數(shù)值加上這個差值；最后，將所得的新模型參數(shù)值設(shè)置為新的當前模型以進行新的迭代運算。
全文摘要
本發(fā)明屬于超聲檢測技術(shù)領(lǐng)域,特別涉及一種利用聲波探測物體內(nèi)部結(jié)構(gòu)的方法。在阻尼最小二乘法的基礎(chǔ)上,通過加入與標準阻尼最小二乘法不同的階梯式修正法和結(jié)構(gòu)重整步驟,利用聲波波動方程的全反演算法,多次重復(fù)整體迭代直至將被探測物體的內(nèi)部聲學參數(shù)計算出來。應(yīng)用它可比現(xiàn)有的基于單向或雙向波場理論的反演算的算法探測出更復(fù)雜和更深層的內(nèi)部結(jié)構(gòu)。本發(fā)明的方法在接收到的聲波信號有一定誤差時仍能得到穩(wěn)定的正確結(jié)果。
文檔編號G01N29/00GK1285510SQ9911151
公開日2001年2月28日 申請日期1999年8月18日 優(yōu)先權(quán)日1999年8月18日
發(fā)明者孫剛, 常謙順, 沈平 申請人:中國科學院物理研究所, 中國科學院應(yīng)用數(shù)學研究所, 香港科技大學