CN1139805C - 一种利用声波探测物体内部结构的方法 - Google Patents

Abstract

本发明属于超声检测技术领域，特别涉及一种利用声波探测物体内部结构的方法。在阻尼最小二乘法的基础上，通过加入与标准阻尼最小二乘法不同的阶梯式修正法和结构重整步骤，利用声波波动方程的全反演算法，多次重复整体迭代直至将被探测物体的内部声学参数计算出来。应用它可比现有的基于单向或双向波场理论的反演算的算法探测出更复杂和更深层的内部结构。本发明的方法在接收到的声波信号有一定误差时仍能得到稳定的正确结果。

Description

一种利用声波探测物体内部结构的方法

                      技术领域

本发明属于超声检测技术领域，特别涉及一种利用声波探测物体内部结构的方法。

                      背景技术

利用发声源使声波在被探测物体中传播，并通过声接收器记录下在被探测物体中传播后的声波信号，再经过反演算方法确定出被探测物体的内部结构的方法现在已经有了一些应用。在这种探测方法中反演算的算法是整个方法的关键。

我们知道声波在物体中的传播一般是用弹性波动方程来描述的。而当切向波可以忽略或传播距离较小时，声信号可用一个标度波动方程 $\frac{\mathrm{\ρ}\left(x\right)}{\mathrm{\κ}\left(x\right)}\frac{{\∂}^2{P}^s(x,t)}{{\∂t}^2}=\mathrm{\ρ}\left(x\right)\▿\·[\frac{1}{\mathrm{\ρ}\left(x\right)}\▿P^s(x,t)]+\mathrm{\δ}(x-x_s)f\left(t\right),----\left(1\right)$ 近似地描述。式中：x为空间的位置矢量；t为时间；P^s(x，t)为由第s发声源所产生的声压；x_s为第s发声源的位置矢量；f(t)为源所产生的脉冲波形；ρ(x)为物质内部的密度，κ(x)为物质内部的压缩模量；式中表示对空间求梯度，而·表示对空间求散度。这里我们将近一步假定物质内部的密度ρ(x)＝ρ为常数。这样方程(1)中将只有一个参量 $K\left(x\right)=\frac{\mathrm{\κ}\left(x\right)}{\mathrm{\ρ}}={c\left(x\right)}^2,$ 其中c(x)为物体内的压缩声波的传播速度。如果我们通过运算能够求出K(x)或c(x)，我们就可以在某种程度上确定被探测物体的内部结构。

目前在大部份的应用中，如：在探测地质结构的应用中所用的反演算的算法主要是基于单向波场理论，正如C.P.A.Wapenaar和A.J.Berkhout，《弹性波场的演化》(Elastic wave fi eld extrapolation)，(Elsevier，Amsterdam，1989)所指出的，这一理论首先将声场分解成正向传播P^s ₊(x，t)和负向传播P^s _-(x，t)的两部分，并将它们写成两个相互有偶合的一阶偏微分方程组。这个有相互偶合的一阶偏微分方程组与方程(1)是等价的。但在实际处理中，为了使实际运算得以进行，在计算中需要忽略这个相互偶合项。忽略这个相互偶合项在物理上意味着忽略了多次反射效应。而这一点从根本上限制了基于单向波场理论的反演算方法的发展，使这种方法难以探测较多的岩层。这是因为经过多层后传出的声波信号总是与多次反射的信号混在一起的。这种方法可探测的标准层数一般为3-5层。

为了能够探测出更多的层数，严格按照声波波动方程(1)进行的全反演算的算法有着更光明的前景。严格按照声波波动方程(1)进行的全反演算的基本手段之一是非线性最小二乘法。它的计算目的是使平方和 $\underset{r,s,t}{\mathrm{\Σ}}{\left[{\mathrm{\δd}}^s(x_r,t)\right]}^2=\underset{r,s,t}{\mathrm{\Σ}}{[P^s(x_r,t)-P_m^s(x_r,t)]}^2----\left(2\right)$ 最小。其中：δd^s(x_r，t)≡P^s(x_r，t)-P^s _m(x_r，t)；

表示对所有发声源，接受器及时间求和，x_r为第r接收器的位置矢量，P^s(x_r，t)表示在假定模型的声学参数为K(x)的条件下按(1)式计算出的压力，而P^s _m(x_r，t)表示实际测量的压力。

非线性最小二乘法与线性最小二乘法的最大区别为：线性最小二乘法可一次求得满足公式(2)的条件，这里即：K(x)的值，且结果是稳定的；而非线性最小二乘法则需要在某一初始值进行多次的迭代，且最终并不一定能得到正确的结果。在我们的问题中假定初始的模型的声学参数为K(x)，非线性最小二乘法用来修改K(x)的公式为A^TA(δK)＝A^T(δd)，                         (3)这里A是偏微商矩阵，阶数为L×N_x。其中：N_x为描写K(x)所需的参数的个数，如：在完整的有限元空间描写K(x)所需的参数为将K(x)分立化时所产生的点数；而L＝n_r×n_s×n_t，其中：n_r和n_s分别为在反演算中所用的接收器和发声源的个数，n_t为在反演算中所使用的分立的时间的步数。A^T为矩阵A的转置。δd为L阶矢量，它是公式(2)中的δd^s(x_r，t)的矢量描述形式。δK＝δK(x)为一个N_x阶的矢量，它是声学参数K(x)的修正量。这里K(x)也为N_x阶矢量，它是声学参数的矢量描述形式。修正后的声学参数为K(x)δK(x)。在作完此修正后，将把修正后的声学参数做为新的初始模型K(x)，并重复上述修正过程。

在实际的反演算运算过程中，最小二乘法的修正式(3)对于δK经常是病态的，即：解不唯一。因而，在实际计算中需将原始的最小二乘法作一些改进使它能有效的处理病态问题。阻尼最小二乘法是这类改进中的一种，应用它可在一定程度上改善方程的病态性质。阻尼最小二乘法的主要原理是在原始最小二乘法的基础上加入一个能使声学参数K(x)在一定程度上保持原有数值的衰减系数γ。在阻尼最小二乘法中，改进声学参数K(x)的公式为(A^TA+γI)(δK)＝A^T(δd).                      (4)其中，I为单位矩阵。在计算中我们需适当调整衰减系数γ的数值，使每次K(x)的改进向好的方向进行。

利用这种非线性阻尼最小二乘法，在某一初始的声学参数下(如：将内部结构设置成与表面相等的均匀值)，经过M₀次的迭代，逐渐改变模型的声学参数，最终得到与被探测物体的内部声学参数一致的K(x)的完整流程图，如图1所示，其中，K(x)＝K(x)δK(x)代表阻尼最小二乘法。但是，在实际的全反演算的模拟计算中，这种简单的应用最小二乘法的反演算算法没能达到人们所期望的效果，如F.Jurado，M.Cuer和V.Richard，在《地球物理》(Geophysics)，60，1857(1995)；W.W.Symes和Carazzone，在《地球物理》(Geophysics)，56，654(1991)中分析失败原因所指出的，人们认为主要是由于在反演算的过程中，方程(3)或(4)通常具有病态的性质，以及残差的平方和，即(2)式具有众多的局域极小值。这一系列问题严重地妨害了反演算的正确进程。而如何解决这些问题至今还没有找到有效的答案。为此，利用全反演算的算法来探测物体内部结构，至今还处于探索阶段。

                    发明目的

本发明的目的在于提供一种利用声波探测物体内部结构的方法，其是利用发声源在被探测物体表面发声，使声波在被探测物体中传播，并在被探测物体的同一表面用声接收器接收从物体内部传出的声信号，将发声源所发的声波的特性，及声接收器接收到的在被探测物体中传播后的声信号，通过基于声波波动方程进行的全反演算的算法将被探测物体的内部声学参数(这里为物体内的压缩波的声速或等价的参数)计算出来。在计算中加入了与传统修正方式不同的步骤，有效地防止了(2)式的残差平方和落入局域极小值，从而可准确地确定出被探测物体的内部结构。

本发明利用发声源、声接收器、以及通过声波波动方程的全反演算法对物体内部结构进行探测，可用来对地质结构、海洋中的声速水平分布、大气分层等进行探测。

本发明包括计算机、发声源、声接收器和被探测物体，其特征在于：

(1).将声接收器均匀分布放置在被探测物体的某一表面上，记录下它们的位置；在放置声接收器的同一表面上选定拟放发声源的位置，并记录下它们的位置；将在每一处接收到的由发声源所发出的声压存入计算机，当对所有拟放发声源的位置进行发声和接收记录后，结束数据采集工作；

(2).在基于阻尼最小二乘法，利用声波波动方程的全反演算法的基础上，通过加入与标准阻尼最小二乘法不同的阶梯式修正法和结构重整步骤，且每次的阶梯式修正法，阻尼最小二乘法和结构重整过程都对当前模型的声学参数进行修正，并将每次修正后所得的K(x)设置为新的当前模型的声学参数，经多次重复整体迭代，使 $\underset{r,s,t}{\mathrm{\Σ}}{[P^s(x_r,t)-P_m^s(x_r,t)]}^2$ 的值更有效地趋于整体的最小值，从而确定被探测物体的内部声学参数。

所述的整体迭代，首先将声学参数K(x)的初始结构设置为内部结构与表面值相等的均匀值，整体迭代的指标J设置为1；初始化完成后，整体迭代部分首先设置对阶梯式修正法的迭代的指标J₁为0，这时进入一个J₁是否小于M₁的判断，当J₁确实小于M₁时，进行阶梯式修正法的修正；修正完毕后，将迭代的指标J₁在现有值上加1，并再次进入J₁是否小于M₁的判断；如果J₁小于M₁，重复进行上面的过程直至J₁等于M₁；完成了整体迭代中的第一部分后，则进入整体迭代中的第二部分，在此设置对阻尼最小二乘法的迭代指标J₀为0，然后进入一个对J₀是否小于M₀的判断，当J₀确实小于M₀时，进行阻尼最小二乘法的修正；修正完毕后，将迭代指标J₀在现有值上加1，并再次进入J₀是否小于M₀的判断；如果J₀小于M₀，重复进行上面的过程直至J₀等于M₀；完成了整体迭代中的前两部分后，将进入对整体迭代的指标J是否小于M的判断，当J确实小于M时，进行结构重整；

结构重整后，将整体迭代指标J在现有值上加1，并回到下一次的整体迭代；

多次重复整体迭代直至J等于M，当J＝M后，将进行最后一次对整体迭代中的前两部分的运算，而不再作结构重整，然后将当前声学参数K(x)的结果输出，并同时将 $\underset{r,s,t}{\mathrm{\Σ}}{[P^s(x_r,t)-P_m^s(x_r,t)]}^2$ 的值输出，用以检验最后结果的正确性，结束利用声波探测物体内部结构的方法。

所述的阶梯式修正法是通过对由阻尼最小二乘法所求得的修正量δK(x)进行一个相应的变换来求得当前模型的声学参数的修正量δK′(x)，其首先确定由标准的阻尼最小二乘法所求得的修正量δK(x)的极值的位置及极值的性质，当两极值的距离足够接近时，确定它们的中点为起始修正点，它们的深层的极值和浅层的极值的差值为该起始修正点之后的修正值；找出满足上述条件的所有修正值中的最大值δK′_max，设定一个小于1的相对域值τ，对所有修正值大于τδK′_max的地方进行修正，其为深于起始修正点的所有声学参数值加上它相应的修正值，在每一次计算中可同时修正多个这样的起始修正点和修正值，这些修正的总和构成δK′(x)。

所述的结构重整是首先将现有的声学参数做调整，设定一个小的域值和一个新声学参数，其新声学参数的表面值与现有声学参数的表面值相等，从表面开始逐次做参数重整，直至整个新声学参数被确定；在某一当前位置，当现有声学参数值与浅层邻近的现有声学参数值之差的绝对值小于或等于这个域值时，则设定这一位置的新声学参数值与它浅层的新声学参数值相等，而当这个差值的绝对值大于这个域值时，则设定这一位置的新声学参数值为浅层邻近的新声学参数值加上这个差值；最后，将所得的新声学参数值设置为新的当前模型以进行新的迭代运算。

本发明由4部分组成：1.计算机；2.发声源；3.声接收器；4.被探测物体。其中，每一部分具有如下特征。

计算机：有一定计算能力的计算机。

发声源：可发出某种特殊脉冲波形的发声源。如：电声换能器、爆炸等。

声接收器：可在某特定点实时记录该点的声压的声接收器。如：电声换能器、声强计、实时声强分析仪等。

被探测物体：被探测物体由具有不同声学参数的物质构成。如：由具有不同声速的物质构成的物体。(注：本发明是根据声学参数的分布来确定物质内部结构的，相同声学参数的不同物质在本发明中是无法区分的)。

本发明的工作步骤如下：

将声接收器尽可能均匀分布地放置在被探测物体的一侧表面上，并记录下它们的位置，做好接收声信号的准备。在放置声接收器的同一表面上选定拟放发声源的位置，并记录下它们的位置。拟放发声源的位置也应尽可能均匀地分布在这一表面。拟放发声源的位置、放置声接收器的位置和被探测物体的关系如图2所示。

逐一在拟放发声源的位置x_s放置发声源，并使其发放一脉冲声波f(t)。从发声源开始发声时起直至声波在物体内部消失的时段里，记录所有声接收器接收到的实时声压波形。将在每一x_r处接收到的由第s发声源所发出的声压记为P_m ^s(x_r，t)，并存入计算机。当对所有拟放发声源的位置进行过前述的发声和接收记录后，数据的采集工作完结。

在数据采集完结后，本发明将由计算机根据采集到的数据，通过声波波动方程的全反演算法将被探测物体的内部声学参数计算出来。我们所采用的算法是在实时空间的阻尼最小二乘法(图1)的基础上，通过加入一个与标准阻尼最小二乘法不同的称做阶梯式修正法的算法，和称做结构重整的算法，使公式(2)的残差平方和更有效地趋于整体的最小值。具体计算过程如图3所示。

图3中的初始结构的作用为将当前模型的声学参数K(x)设置为内部结构与表面值相等的均匀值。

整个计算过程包括M次的整体迭代(对图3中的指标J)。在每次整体迭代又包括M₁次的对阶梯式修正法的迭代(对图3中的指标J₁)，M₀次的对阻尼最小二乘法的迭代(对图3的中指标J₀)和一次结构重整。每次的阶梯式修正法，阻尼最小二乘法和结构重整过程都按照下面介绍的特定方式对当前模型的声学参数进行修正，并将修正后所得的K(x)设置为新的当前模型的声学参数。

图3中的阻尼最小二乘法对应于标准的阻尼最小二乘法。它的声学参数修正量δK(x)可由标准的阻尼最小二乘法线性方程组，即公式(4)求得。

图3中的阶梯式修正法是本发明所特有的修正法。它对当前模型的声学参数的修正量δK′(x)是通过对由标准的阻尼最小二乘法所求得的修正量δK(x)进行一个相应的变换来求得的。具体变换的规则如下。首先确定由标准的阻尼最小二乘法所求得的修正量δK(x)的极值的位置及极值的性质(极大值或极小值)。当两极值的距离足够接近时(如小于等于3)，我们确定它们的中点为起始修正点，它们的差值(深层的极值减去浅层的极值)为该起始修正点之后的修正值。找出满足上述条件的所有修正值中的最大值δK′_max。设定一个相对域值τ(如等于0.8)。对所有修正值大于τδK′_max的地方做如下修正，即深于起始修正点的所有声学参数值将加上它相应的修正值。在每一次计算中可同时修正多个这样的起始修正点和修正值，这些修正的总和构成δK′(x)。

图3中的结构重整将现有的声学参数做如下的调整。设定一个小的域值(如等于0.01)和一个新声学参数。新声学参数的表面值与现有声学参数的表面值相等。从表面开始逐次做下述重整过程直至整个新声学参数被确定。在某一当前位置，当现有声学参数值与浅层邻近的现有声学参数值之差的绝对值小于或等于这个域值时、则设定这一位置的新声学参数值与它浅层的新声学参数值相等，而当这个差值(当前位置的现有声学参数值减去浅层邻近的现有声学参数值)的绝对值大于这个域值时，则设定这一位置的新声学参数值为浅层邻近的新声学参数值加上这个差值。最后，将所得的新声学参数值设置为新的当前模型的声学参数以进行新的迭代运算。

在进行阶梯式修正法和阻尼最小二乘法的修正中都将用到阻尼最小二乘法的衰减系数γ，它的值在运算中被调整至使修正后的残差平方和，即(2)式比修正前的残差平方和要小或相等。由于当衰减系数趋于无穷大时意味着没有修正，这时残差的平方和不变。因此，逐渐加大衰减系数使残差平方和变小或不变显然总是可以做到的。

在实际计算中可根据计算机的计算能力确定一组合适的M、M₁和M₀。典型的M₁和M₀的值约为10左右。M的值越大结果越准确可靠。

按照图3的流程图以及前面的说明，声学参数将按下述的方式逐渐逼近物体所具有的真实值。首先、将声学参数K(x)的初始结构设置为内部结构与表面值相等的均匀值(如图4)。同时将整体迭代次数、阶梯式修正法的迭代次数和阻尼最小二乘法的迭代次数分别确定为M、M₁和M₀，并将整体迭代的指标J设置为1。初始化完成后，计算进入整体迭代部分。在整体迭代部分首先设置对阶梯式修正法的迭代的指标J₁为0。这时进入一个J₁是否小于M₁的判断。当J₁确实小于M₁时，进行阶梯式修正法的修正。我们称这个修正法为阶梯式修正法是因为这个修正法每次给出的修正值均为一种阶梯形式，这一点在从最早的初始模型(即：与表面值相等的均匀值)进行第一次的阶梯式修正时最为明显(如图5)。从图5中可以看出修正后的模型的声学参数从原来的均匀值变为一个阶梯形式。修正完毕后，将迭代的指标J₁在现有值上加1，并再次进入J₁是否小于M₁的判断。如果J₁小于M₁，重复进行上面的过程直至J₁等于M₁。当J₁等于M₁时，我们应已经进行了M₁次的阶梯式修正法的修正。这时模型的声学参数已为M₁次的阶梯式修正的叠加，因此阶梯形式已不如单一次的修正时(图5)明显。但浅层的大模样已接近物体所具有的真实值，但在浅层仍存在细小的差别(图6)。完成了整体迭代中的第一部分后，我们将进入整体迭代中的第二部分。在此我们设置对阻尼最小二乘法的迭代指标J₀为0。然后进入一个对J₀是否小于M₀的判断。当J₀确实小于M₀时，进行阻尼最小二乘法的修正。阻尼最小二乘法的修正的作用是在前面的计算基础上使浅层仍存在的细小差别逐渐消失。修正完毕后，将迭代指标J₀在现有值上加1，并再次进入J₀是否小于M₀的判断。如果J₀小于M₀，重复进行上面的过程直至J₀等于M₀。当J₀等于M₀时，我们应已经进行了M₀次的阻尼最小二乘法的修正。这时模型的声学参数在浅层的细小差别已明显减少，而深层则变为一种连续变化与阶梯变化的混合形式(图7)。从图7中可以看到在深层的阶梯变化形式与物体所具有的真实值相似，而连续变化部分与真实值不符。完成了整体迭代中的前两部分后，我们将进入对整体迭代的指标J是否小于M的判断。当J确实小于M时，进行结构重整。结构重整的作用是将与真实值不符的连续变化部分滤掉，使模型的声学参数仍为一种阶梯变化的形式(图8)。结构重整后，将整体迭代指标J在现有值上加1，并回到下一次的整体迭代。每次整体迭代会使更深层的声学参数与物体所具有的真实值一致。图9显示了经过4次整体迭代后的模型的声学参数。该图显示深层的声学参数已经与真实值很接近，它比经过1次整体迭代后的模型的声学参数(图8)改进了许多。多次重复整体迭代直至J等于M。当J＝M后，我们将进行最后一次对整体迭代中的前两部分的运算(即：对阶梯式修正法的迭代和对阻尼最小二乘法的迭代)，而不再作结构重整，然后将当前模型的声学参数做为最后结果输出。图10为模型的声学参数的最后结果的示意图。这个结果与物体所具有的真实值完全一致。

本发明的最大特点是在全反演算的算法中加入了不同于现有最小二乘法的阶梯式修正法和结构重整。使基于声波波动方程进行的全反演算的算法有了实际的效果。应用它可比现有的基于单向或双向波场理论的反演算的算法探测出更复杂和更深层的内部结构。本发明的方法在接收到的声波信号有一定误差时仍能得到稳定的正确结果，这是其它的利用全反演算的算法探测物体内部结构的方法所作不到的。下面结合附图和实施例对本发明的技术方案作进一步的描述。

                     附图说明

图1.基于简单最小二程法的算法流程图，其中J₀为循环指标，M₀为循环次数；K(x)＝K(x)+δK(x)为阻尼最小二乘法；

图2.本发明所对应的被探测物体K₀、发声源和声接收器的相对位置示意图；

图3.本发明的算法流程图，其中J、J₁和J₀为循环指标，M、M₁和M₀为循环次数；K(x)＝K(x)+δK(x)为阶梯式修正法，K(x)＝K(x)+δK(x)为阻尼最小二乘法；

图4.初始模型的声学参数的示意图，图中实线为初始结构，点线为真实值的待测物体的声学参数；

图5.从最早的初始模型的声学参数进行第一次阶梯式修正后的计算结果，图中实线为计算结果，点线为真实值的待测物体的声学参数；

图6.从最早的初始模型的声学参数进行M₁次的阶梯式修正后的计算结果，图中实线为计算结果，点线为真实值的待测物体的声学参数；

图7.在图6的基础上进行M₀次的阻尼最小二乘法的修正后的计算结果，图中实线为计算结果，点线为真实值的待测物体的声学参数；

图8.从最早的初始模型的声学参数完成一次整体迭代后的计算结果，图中实线为计算结果，点线为真实值的待测物体的声学参数；

图9.从最早的初始模型的声学参数完成4次整体迭代后的计算结果，图中实线为计算结果，点线为真实值的待测物体的声学参数；

图10.最终的计算结果，图中实线为计算结果，点线为真实值的待测物体的声学参数；

图11.实施例1中的被探测柱状物体的截面图，K_m为胶冻物体，其尺寸为40厘米×40厘米，A_m为吸收层。

            具体实施方式

实施例1

将本发明的方法具体应用到探测一个柱状胶冻物体，它的截面图如图11所示。对于一个柱状物体，如果它的长度远大于它的宽度，则在柱中心附近的截面内其波动方程可用二维波动方程近似地描述。对于胶冻物体切向波可忽略，因而声波仅为压缩波。为了简化问题，被探测物体由密度近似相等的胶冻物体制作，且压缩模量仅随深度变化。这样本被探测物体只有一个声学参数K(x)，且它仅为深度的函数。适当选取时间单位使表面的声学参数K(x)为1.0。在这一时间单位下，准备一个可产生 $f\left(t\right)=(\frac{{\mathrm{\π}}^2{(t-10)}^2}{8^2}-\frac{1}{2})\exp (-\frac{{\mathrm{\π}}^2{(t-10)}^2}{8^2})$ 的脉冲声波的电声换能器作为发声源。选择这样的脉冲声波的理由是这种脉冲形式与爆炸所产生的脉冲形式类似，因而使本实施例的结果可类比用爆炸来探测地质结构的情况。另外准备10个电声换能器作为声接收器，并将其均匀放置在被探测胶冻物体的中心横截面的表面上，其具体位置为从胶冻物体的中心横截面的表面的左端开始4cm、8cm、12cm、16cm、20cm、24cm、28cm、32cm、36cm和40cm处。类似地、在胶冻物体的中心横截面的表面上均匀选取拟放发声源的位置，其具体位置为从胶冻物体的中心横截面的表面的左端开始1cm、5cm、9cm、13cm、17cm、21cm、25cm、29cm、33cm和37cm处。逐一在拟放发声源的位置放置作为发声源的电声换能器，并使其产生脉冲为f(t)的声波。从t＝0时刻开始至t＝512的时段里(t＝512时胶冻物体内部的声波已基本消失)，用作为声接收器的各个电声换能器接收其所在位置的实时声压波形。将在x_r处的电声换能器接收到的由第s发声源所发出的声压波形记为P_m ^s(x_r，t)，并存入计算机。当对所有拟放发声源的位置进行过前述的发声和接收记录后，数据的采集工作完结。电声换能器的相对误差小于5％。

在数据采集完结后，将由计算机根据图3中全反演算的算法，利用现已采集到的数据将被探测物体的内部声学参数计算出来。选取图3中的M＝8、M₁＝M₀＝10。初始模型的声学参数设置为内部结构与表面值相等的均匀值，在此即K(x)＝1.0(如图4)。初始结构确定后将整体迭代的指标J设置为1，并进入整体迭代部分。在整体迭代部分首先设置对阶梯式修正法的迭代指标J₁为0。这时进入一个J₁是否小于M₁的判断。当J₁确实小于M₁时，进行阶梯式修正法的修正。从最早的初始结构经过第一次的阶梯式修正后，阶梯式修正法给出了一个阶梯形式声学参数(如图5)。修正完毕后，将迭代指标J₁在现有值上加1，并再次进入J₁是否小于M₁的判断。如果J₁小于M₁，重复进行上面的过程直至J₁等于M₁。当J₁等于M₁时，我们已经进行了M₁次的阶梯式修正法的修正。这时模型的声学参数已为M₁次的阶梯式修正的叠加，因此阶梯形式已不如单一次的修正时(图5)明显。此时在浅层，模型的声学参数的大模样已接近物体所具有的真实值，但仍存在细小的差别，而在深层则存在较大的差别(图6)。完成了整体迭代中的第一部分后，我们进入整体迭代中的第二部分。在此我们设置对阻尼最小二乘法的迭代指标J₀为0。然后进入一个对J₀是否小于M₀的判断。当J₀确实小于M₀时，进行阻尼最小二乘法的修正。修正完毕后，将迭代指标J₀在现有值上加1，并再次进入J₀是否小于M₀的判断。如果J₀小于M₀，重复进行上面的过程直至J₀等于M₀。当J₀等于M₀时，我们已经进行了M₀次的阻尼最小二乘法的修正。这时模型的声学参数在浅层的细小差别已明显减少，而深层则变为一种连续变化与阶梯变化的混合形式(图7)。完成了整体迭代中的前两部分后，我们进入对整体迭代的指标J是否小于M的判断。当J确实小于M时，进行结构重整。结构重整后，与真实值不符的连续变化部分被部分地滤掉，使模型的声学参数仍变回到一种阶梯变化的形式(图8)。完成结构重整后，将整体迭代指标J在现有值上加1，并回到下一次的整体迭代。每次整体迭代会使更深层的声学参数与物体所具有的真实值更趋于一致。经过4次整体迭代后的模型的声学参数如图9所示。该图显示深层的声学参数已经与真实值很接近，它比经过1次整体迭代后的声学参数(图8)改进了许多。多次重复整体迭代直至J等于M。当J＝M后，即：J＝8，我们进行最后一次对整体迭代中的前两部分的运算(即：对阶梯式修正法的迭代和对阻尼最小二乘法的迭代)，而不再作结构重整，然后将当前模型的声学参数做为最后结果输出。图10显示了本实施例的最终结果。从图10可以看到经8次整体迭代后，实线和点线完全吻合，表明被测胶冻物体中8层不同声学参数的区域均被准确探测到。

在本实施例中，胶冻物体共有8层不同参数的区域，大大超出了应用基于单向波场理论所能探测的标准层数(3-5层)。此外、由于电声换能器发出和接收到的声波信号有一定误差，而本发明的方法在此情况下仍能得到稳定的正确结果。这是其它的利用全反演算的算法探测物体内部结构的方法所作不到的。

Claims (4)

1.一种利用声波探测物体内部结构的方法，包括计算机、发声源、声接收器和被探测物体，其特征在于：

(1).将声接收器均匀分布放置在被探测物体的一侧表面上，记录下它们的位置；在放置声接收器的同一表面上选定拟放发声源的位置，并记录下它们的位置；将在每一处接收到的由发声源所发出的声压存入计算机，当对所有拟放发声源的位置进行发声和接收记录后，结束数据采集工作；

(2).在基于阻尼最小二乘法，利用声波波动方程的全反演算法的基础上，通过加入与标准阻尼最小二乘法不同的阶梯式修正法和结构重整步骤，且每次的阶梯式修正法，阻尼最小二乘法和结构重整过程都对当前模型的声学参数进行修正，并将每次修正后所得的声学参数K(x)设置为新的当前模型的声学参数，经多次重复整体迭代，使 $\underset{r,s,t}{\mathrm{\Σ}}{[P^s(x_r,t)-P_m^s(x_r,t)]}^2$ 的值更有效地趋于整体的最小值，从而确定被探测物体的内部声学参数。

2.如权利要求1所述的利用声波探测物体内部结构的方法，其特征在于：所述的整体迭代，首先将声学参数K(x)的初始结构设置为内部结构与表面值相等的均匀值，整体迭代的指标J设置为1；初始化完成后，整体迭代部分首先设置对阶梯式修正法的迭代的指标J₁为0，这时进入一个J₁是否小于M₁的判断，当J₁确实小于M₁时，进行阶梯式修正法的修正；修正完毕后，将迭代的指标J₁在现有值上加1，并再次进入J₁是否小于M₁的判断；如果J₁小于M₁，重复进行上面的过程直至J₁等于M₁；完成了整体迭代中的第一部分后，则进入整体迭代中的第二部分，在此设置对阻尼最小二乘法的迭代指标J₀为0，然后进入一个对J₀是否小于M₀的判断，当J₀确实小于M₀时，进行阻尼最小二乘法的修正；修正完毕后，将迭代指标J₀在现有值上加1，并再次进入J₀是否小于M₀的判断；如果J₀小于M₀，重复进行上面的过程直至J₀等于M₀；完成了整体迭代中的前两部分后，将进入对整体迭代的指标J是否小于M的判断，当J确实小于M时，进行结构重整；

结构重整后，将整体迭代指标J在现有值上加1，并回到下一次的整体迭代；

多次重复整体迭代直至J等于M，当J＝M后，将进行最后一次对整体迭代中的前两部分的运算，而不再作结构重整，然后将当前模型的声学参数K(x)的结果输出，并同时将 $\underset{r,s,t}{\mathrm{\Σ}}{[P^s(x_r,t)-P_m^s(x_r,t)]}^2$ 的值输出，用以检验最后结果的正确性，结束利用声波探测物体内部结构的方法。

3.如权利要求1或2所述的利用声波探测物体内部结构的方法，其特征在于所述的阶梯式修正法是通过对由阻尼最小二乘法所求得的修正量δK(x)进行一个相应的变换来求得当前模型的修正量δK′(x)，其首先确定由标准的阻尼最小二乘法所求得的修正量δK(x)的极值的位置及极值的性质，当两极值的距离足够接近时，确定它们的中点为起始修正点，它们的深层的极值和浅层的极值的差值为该起始修正点之后的修正值；找出满足上述条件的所有修正值中的最大值δK′_max，设定一个小于1的相对域值τ，对所有修正值大于τδK′_max的地方进行修正，其为深于起始修正点的所有声学参数值加上它相应的修正值，在每一次计算中可同时修正多个这样的起始修正点和修正值，这些修正的总和构成δK′(x)。

4.如权利要求1或2所述的利用声波探测物体内部结构的方法，其特征在于所述的结构重整是首先将现有模型的声学参数做调整，设定一个小的域值和一个新声学参数，其新声学参数的表面值与现有声学参数的表面值相等，从表面开始逐次做参数重整，直至整个新声学参数被确定；在某一当前位置，当现有声学参数值与浅层邻近的现有声学参数值之差的绝对值小于或等于这个域值时，则设定这一位置的新声学参数值与它浅层的新声学参数值相等，而当这个差值的绝对值大于这个域值时，则设定这一位置的新声学参数值为浅层邻近的新声学参数值加上这个差值；最后，将所得的新声学参数值设置为新的当前模型的声学参数以进行新的迭代运算。

Images