分布式目标方位估计具体讲解

1、西北工业大学硕士论文 第二章 分布式目标回波分析方法 第二章 分布式目标回波分析方法 当声波在传播中遇到目标时， 会在目标表面激发起次级声源， 它们向周围介 质中辐射次级声波，我们把这些次级声波统称为散射波。在收发合置的情况下， 返回声源方向的那部分散射波被称为目标回波。 目标回波与目标的散射特征密切 相关，它是入射波与目标相互作用后才产生的，在此过程中， 有关目标本身的某 些特征信息会包含在回波中， 人们通过对回波的分析处理， 可以将目标的特征信 息提取出来，再辅以某些先验知识，从而实现对目标的探测。例如，通过接收基 阵采集目标对某种入射激励信号的回波， 再对阵列的快拍数据进行分析处理， 就

2、 可以估计出目标的个数、方位、距离、速度等参数。由此可见，研究目标回波的 特性，在工程上是很有意义的。 目标回波不仅与目标本身的各种特征有关，如形状、大小、材料、结构等， 还与入射波的形式有关，入射波时宽、带宽将直接决定回波的结构。 宽带入射信 号将产生宽带回波，窄带入射信号引起窄带回波，时宽入射信号产生稳态响应， 窄脉冲信号产生瞬态响应。 本文主要就高频窄带信号在稳态和瞬态两个方面进行 介绍。此外，目标回波一般来说包括几何散射波和弹性散射波， 几何散射波是由 目标本身的外形特征决定的， 与目标的材料无关； 而弹性散射波是由于目标的弹 性声阻抗与介质不同而引起的， 与目标本身的材料有关系。 本

3、文所关心的主要是 目标的几何特征，因此，以后将对几何回波作着重研究。 2.1 离散散射型目标回波 所谓离散散射型目标是把目标当作若干个散射很强的亮点组成的的散射体 , 目标的回波是由所有亮点的子回波的线性叠加。式 (1.1 )、(1.2 )分别给出了单 亮点和多亮点目标的传递函数， 要得到目标的传递函数， 必须先确定式中的三组 参数，即幅度因子 Am、 时间延迟 m、 相位跳变 m( m=1，2，,N )，人们通过 研究发现了一些几何体的参数。 1凸光滑曲面的镜反射 西北工业大学硕士论文 第二章 分布式目标回波分析方法 设 Pi( ) 为镜反射亮点处入射波的频谱， 则由该亮点产生的镜反射波的频

4、谱 为： Pb(r, ) 0.5| R1 | R2 |e .Pi( ) (2.1) b (1 |R1|/r)(1 |R2 |/r) r i 其中 ,R1、R2为镜反射亮点处的目标表面的两个主曲率半径， r 为空间散射场中某 一点与亮点之间的距离， k 为入射波波数， ejkr /r 用于修正信道传输的影响，因 为无论目标多么复杂，在远场中散射波都以 ejkr /r 的规律扩展。三个亮点参数分 别为： A(r, ) 0.5| R1 | R2 |， =0， =0(2.2 )若 (1 |R1 | / r )(1 |R2 |/r) 满足 rR1及 rR2，则 A | R1 |R2 | 2(2.3) 2

5、有限长圆柱体 如图 2-1 所示 , 为声波相对于圆柱 体的入射角 , L 为圆柱体长度， 2为直径， 取棱角 1 为参考点，则被照射到的三个亮点 的散射参数为： /4 2 3 /4 (2.4a ) (2.4b ) sin1/2 cos ， 3 2Lcos /c ， 3 5 /4 (2.4c ) 西北工业大学硕士论文 第二章 分布式目标回波分析方法 可见，幅度因子和时间延迟是波数 k、半径 a、圆柱体长度 L以及入射角 的函 数，而相位跳变是一个常数。 3无限长圆柱体 根据 Sommerfeld-Watson 变换分析和实验结果，无限长圆柱回波包括镜反 射波和各型表面环绕波组成， 在远场高频条

6、件下， 回波的主要成分为镜反射回波。 jkr e ) a / 2 M r (2.5 ) 镜反射回波可以表示为： P sp ( r , 镜反射亮点反射参数为： Aspa/2M , sp 0 , sp 0 (2.6) 这里参考点取在镜反射点上 ,M 为镜反射系数 , 当高频 (大 ka) 时可以用平面波正 入射到柱体表面的反射系数近似。 一旦确定的目标的亮点散射参数， 就可以根据式 (1.2 )模拟出目标的散射回 波。 2.2 连续散射型目标回波 所谓连续散射型目标是指把分布式目标看作表面各点的散射特性在空间位 置上是连续变化的几何体 , 目标表面各点都对回波有贡献。这种观点更有利于准 确地描述目

7、标的散射特性以及精确探测目标的回波结构。由于表面各点的位置、 法线向量、声阻抗等参数不同，它们对回波贡献也就不同。 对于连续散射型目标 回波是通过一些基于面元法的思想来描述和分析的， 将目标表面的面积微元对回 波的贡献在整个表面进行面积分，就可以获得回波响应。若表面形状很复杂， 对 表面积分就很繁琐，甚至根本不能得出结果，在这种情况下， 可以将表面按某种 方式网格化， 再运用数值计算的方法就可以得出模拟的回波结果， 人们在实际中 发现，这些数值方法得出的结果具有足够的精确度，可以满足工程应用的要求。 西北工业大学硕士论文 第二章 分布式目标回波分析方法 A 设有一目标置于均匀介质空间中，目标的

8、外表面为封闭曲面 S，其外法线向 量为 n, 有一点声源位于点，下面计算介质空间中另一点处的散射声场。 如图 2-2 所示， B 点散射声场的赫姆霍兹积分公式为： 41 s r12 ejkr2. n( s) s n(r12e 1 jkr2 )dS (2.7 ) 其中， k 为波数， s 是散射声场的速度势函数，n 表示在 S 面上的外法线方 向导数, r 1、r2分别是 A、 B点到表面上面元 dS的距离。由于被积函数中的 s及 其法线方向导数n 是未知量，所以，上式的积分一般不能直接求出。但是， 应用边界条件，可将被积函数中的未知量用已知量表示。 假设目标表面 S 是刚性的，则在 S上有 i

9、 s 0 (2.8 ) n n s 其中， i 是入射波势函数，它可表示为 i Aejkr1(2.9 )这 r1 里 A是一个表示振幅的常数(式中略去了时间因子 ej t )。由式 (2.8 )及(2.9 )， 并考虑远场高频条件 kr 11,我们可以得到 西北工业大学硕士论文 第二章 分布式目标回波分析方法 cos(n, r1 ) (2.10) 其中, (n,r1) 表示矢量 r1与法线向量 n之间的夹角。同样地，在 kr 21时有 1 ejkr2jk ejkr2 cos(n,r2)(2.11 ) n r2r2 作为一种足够精确的近似，可以认为刚性目标表面上的散射声场等于入射声场， 于是 s

10、 s A ejkr1(2.12 ) r1 将式(2.10 )、(2.11 )、(2.12 )代入式(2.7 )，就可以得到 s jkA 1 ejk(r1 r2) cos(n,r1) cos(n,r2) dS(2.13 ) 4 s r1r2 在收发合置的情况下， r 1=r 2=r ，则上式成为 s jkA 12 ej2kr cos(r ,n) dS(2.14 )式 2 s r 2 (2.13 )、 (2.14 )即为目标为刚性体时散射声场的高频近似积分公式，只要设法 将积分求得，就可以求出介质空间中任一点处的散射声场。 式(2.13 )的物理意义可以理解为：目标表面各点在入射声波的激励下，作

11、为次级声源辐射出二次级波， 接收点上的散射声场就是目标表面上所有二次源辐 射的声波在该点的叠加，这些二级次波的相位为k(r 1+r 2), 即取决于声波在介质 中的往返波程。 2费涅尔 (Fresnel) 半波带近似法 式(2.13) 、(2.14 )在知道目标表面方程的情况下是可以求得其散射声场的， 但对于复杂目标表面，计算过程将变得相当繁琐，因此， 工程上常会应用费涅尔 半波带近似法来简化计算工作量。 西北工业大学硕士论文 第二章 分布式目标回波分析方法 基于散射声波的相位取决于入射波在目标表面各点的往返声程这一机理， 我 们可以仿照光学中的方法，应用费涅尔半波带法来求解式 (2.14 )

12、。 图 2-3 菲涅尔半波带示意图 如图 2-3 所式，设入射波从 B点发射， C点为目标表面上距 B点最近的点， 设它距 B点的波程为 r0，以 B点为球心，r 0为半径作球面， 它与目标表面相切于 C 点，接着跳跃式地增加球面半径，每次增加 1/4 波长， 即 r 0+ /4 ，r 0+/2 ， r 0+3/4 ，r 0+， ，这些半径不同的球面把目标表面分割成许多环带 S1，S2， S3， S4， ，SN，它们被称为费涅尔带。可以看出，相邻半波带产生的散射波在 B点的声程差为 /2 ，即相位差为 。因此，若目标表面上共有 N个半波带，且 第i 个波带的贡献用i表示，则 s 1 2 3 (

13、 1)N 1 N (2.15 ) 其中， i jkA 12 ej2krr .cos(r , n)dS i 1,2, ,N (2.16 )若 i 2 Si r2 目标的径向尺寸远远大于波长 ,并且目标表面不太弯曲 , 那么,目标表面可以分割 出许多费涅尔半波带 ,而且,相邻波带的 cos(r , n) / r变化也不大 ,彼此的面积也很 接近, 这样就可以认为第 i 个波带产生的散射波等于相邻两个波带产生散射波的 平均值,即 i ( i 1 i 1)/2 i 2, ,N 1 (2.17) 西北工业大学硕士论文 第二章 分布式目标回波分析方法 将上式代入式 (2.15), 得到总的散射波为 s 1

14、 ( 1 3)/ 2 3 ( 3 5)/2 1 1 ( 1)N1 N (2.18) 2 即总的散射波等于第一个与最后一个费涅尔带所产生的散射波的平均值。 对于大体积目标，最后一个费涅尔带的 cos( r , n ) 0 ，因而它的贡献可以忽 略不计，而第一个费涅尔带的 cos(r , n) 1，于是，总的散射波可以简化表示为 S jkA 12 ej2kr dS(2.19 )这 4 S1 r2 样，所得到的积分表达式就比式 (2.14 )容易计算了。 3边界元法( BEM) 边界积分方程在计算声学中的应用可以追溯到二十世纪六十年代。 最早将边 界积分方程应用于任意形状物体的声辐射的是 Chen和

15、 Schweikert 4 ，此后，很 多学者进行了积极的探索。 在绝大多数的这些早期研究工作中， 对物体表面的网 格化要么是应用常量的网格单元形函数，要么将网格单元近似地以平面来代替。 近年来， 随着边界元方法的迅速发展， 人们将二次单元形函数应用到对声学问题 中边界积分方程的求解， 从而使声学物理量以及边界几何特征能够被更精确地描 述。 1) 自由介质空间中的边界积分公式 图 2-4 介质空间中的物体 如图(2.4) 所示,设有一物体 B 位于无限均匀 介质空间 B中, 物体表面为 S, 介质密度和声速分别 为0 、c. 在稳态时，介质空间 B中的声场应该满 足如下关系： (2.20 西北

16、工业大学硕士论文 第二章 分布式目标回波分析方法 这就是著名的赫姆霍兹 (Helmholtz) 方程，其中， 是速度势， k 是入射波中心 频率波数。对于散射问题，速度势由s及i两部分组成，i为入射波速度势 , s为散射波速度势 , 这两部分在 B中都必须满足赫姆霍兹方程。 物体边界 S 上满足的条件具有一般的形式 (2.21 ) n 其中， n被定义为物体表面上的内法线向量， 、为特定常数。散射波速 度势还应当满足萨默菲尔德 (Sommerfeld) 条件, 即 lim r s jk s0 (2.22) rr s (2.23) 上式反映了在无穷远处速度势所满足的关系。 应用格林第二公式或加权

17、残留量公式， 由式(2.20) 可以得出在边界 S 上的边 界积分方程 : C(P) (P) (Q,P) (Q) (Q,P) (Q) dS(Q) 4 i(P) S n n 上式被称为赫姆霍兹积分方程。其中， (P)为声场中某一定点 P处的速度势， Q 为物体边界上的变点， C(P)为一特定的常量，它取决于 P点的位置， i (P)为入 射波在 P点的速度势, (Q,P)为自由空间中的格林函数 ,它可以表示为 : jkr e (Q,P) (2.24) r 其中， r=|Q-P| ，即 r 为 Q、P两点间的距离。 值得注意的是系数 C(P)的选取，当定点 P在介质空间 B中时， C(P)=4;

18、当 P 点位于物体内部时 ,C(P)=0 ；当 P 点位于物体边界上并且在 P 点有唯一的切平 面( 即 P点位于界面上的光滑部分 ), 则 C(P)=2；而当 P点位于界面上的棱角处 时,则 C(P)由下式确定: C(P) 4 (1 r) dS (2.25) n 西北工业大学硕士论文 第二章 分布式目标回波分析方法 实际上, 上式对于处于任何位置的 P 点都适用。由于没有有效的数值方法判断 P 点是否位于棱角上，因此对于有棱角的物体，只要 P点在边界上， 我们一律用式 (2.25 )计算 C(P) 。 2) 存在界面时的边界积分公式 在实际中,有很多目标位于界面附近 ,比如海面舰艇、 地面建

19、筑物、靠近海面 的潜艇等等，而在这种情况下的边界积分公式与在自由介质空间中的情况有所不 同，所以，需要对其进行专门介绍。 如图 2-5 所示,设 SH是一个无限反射平面 ,物体边界上各点在定点 P的散射波 不仅包括动点 Q的贡献,还包括 Q点关于界面 SH的镜像点的贡献 , 所以, 在求解 P 点的速度势时 , 式(2.23) 的(Q,P) 应当替换为 图 2-5 界面附近物体声散射几何关系 (2.26) r1 e jkr e jkr1 H (Q,P) eRH e r H(Q,P)被称为半无限空间格林函数 , 其中,界面反射系数 RH在界面为刚性时等 于 1, 在界面为软界面时等于 1,r=|

20、Q-P|,r 1=|Q-P1| 。另外，在收发合置的情况 下，入射波应当包含直达波和界面反射波，除非入射方向与界面平行。所以，入 射波在介质中传输的扩展因子为 : e jkr e jkr1 T(r,r1)RH(2.27 ) rr1 西北工业大学硕士论文 第二章 分布式目标回波分析方法 如果物体边界与界面有相接触的部分，如图 2-6 所示 , 物体边界可分为两部 图 2-6 与界面相接触的物体声散射几何关系 分, 一部分是 SC，它与界面 SH相接触, 另一部分是 SO，它与界面不接触，那么， 在 这种情况下的边界积分方程应当为： C(P) (P) H(Q,P) (Q) H (Q,P) (Q)

21、dS(Q) 4 i(P) (2.28 ) SO H n n i 与式(2.23) 不同的是,上式采用的是 H，而不是，并且积分是在 SO上进行的。 当P点在 B或B时,系数 C(P)的选取与自由介质空间中的选取方法相同 ,当P点 在 SO上 ,C(P) 由下式确定 : C(P) 4 (1 r) dS (2.29) SO SC n 当 P点位于 SC上,则 C(P)为: (1/r) C(P) (1 RH ) 2 dS (2.30) SO SCn ( 3) 数值计算公式 前面所介绍的求解散射声场的积分方程法需要知道物体表面边界上的速度 势及其法线方向导数 , 但在散射声场被求解出之前 , 处于物体

22、边界上的这些物理 量却是未知的 , 因此,在求解介质空间中的散射声场时 ,必须先得出物体边界上各 点的速度势及其法线方向导数值 , 再将这些求解出的结果代入边界积分公式 , 进 而求出介质空间中任一点处的散射波速度势。 西北工业大学硕士论文 第二章 分布式目标回波分析方法 问题的关键在于求解物体边界上的速度势及其法线方向导数。 应用数值计算 的方法，可以使这个问题得以解决。 将物体的表面边界 S 以某种方式划分为若干小的曲面片， 这些曲面片被称为 表面边界的单元， 在 S上再选取若干个节点， 我们称这个过程为表面边界的网格 化。设表面边界单元数和节点数分别为 M、N，令 P 点分别置于各个节点

23、上，对 应于每一个节点，应用边界积分公式就会得到一个方程，这样， 当 P 遍取边界上 的 N个所有节点，就得到 N 个线性独立方程组成的方程组，其中的未知量为 2N 个(N 个表面边界速度势及其法线方向导数) ，结合边界条件中给出的 N 个已知 量，最终将会得到一个 N 元恰定的线性方程组，通过求解这个方程组，就可以得 到边界上各节点处的速度势及其法线方向导数。下面将具体介绍这一过程。 在讨论应用单元节点插值方法解边界积分方程之前 , 有必要先简要介绍一下 单元局部坐标的概念。下面以任意三角形和四边形为例进行介绍。 如图 2-7 所示,设 xy 平面上一三角形单元 e 的三个顶点坐标为 P1(

24、x 1,y 1)、 P2(x 2,y 2) 、P3(x 3,y 3) ，该三角形内任一点的坐标为( x，y)，则如下变换： 1 2e x x2 x3 y1 y2 1 y3 1 2e x x3 x1 y y3 y1 2e x x1 x1 y1 y2 1 y2 1 (2.31a 西北工业大学硕士论文 第二章 分布式目标回波分析方法 (2.31b ) x (x1 x3) (x2 x3) x3 y (y1 y3) (y2 y3)y3 将三角形单元 e 映射为平面中的等腰直角三角形 。其中， e为三角形 e 的面积。 、 、一般被称为局部坐标，它们与三角形 e 的形状和大小有关， 由于 、满足 +=1

25、的关系，独立的坐标只有两个，所以该变换并 没有改变维数。而且局部坐标与直角坐标 x、y 的选取方法是无关的，这给处理 问题带来了很大方便。 同理，对于任意四边形单元 ,也有类似的转换关系。 如图2-8 所示,四边形单 元 e 被映射成 平面内的中心在原点的单位正方形 ,a 1、a2、a3、 a4 分别对应 于四边形的顶点 A1、A2、 A3、A4，所用的变换为： 4 xxi Ni ( , ) i 41(2.32 )其 yyi Ni( , ) i1 a4(-1 , 1) a3( 1 , 1) a1(-1, -1) a2( 1, -1) 图 2-8 四边形单元的局部坐标映射关系 中，Ni(, )称

26、为双线性插值基函数 ,它们的表达式为 (2.33) N1( , ) (1 )(1 )/4 N2 ( , ) (1 )(1 )/4 N3( , ) (1 )(1 )/ 4 N4 ( , ) (1 )(1 )/4 这样, 对于一个单元上的任意一点 ,总存在一组局部坐标与之对应 , 反之亦然。即 单元上点的全局坐标与局部坐标具有一一对应的关系。 西北工业大学硕士论文 第二章 分布式目标回波分析方法 实际上，以上变换可以看作以节点 ( 三角形或四边形的顶点 ) 坐标为参数 、 Ni(，)为基函数对单元(三角形或四边形)上任一点坐标的插值。对于空 间的曲面网格单元， 我们可以根据等参数单元的变换的办法，

27、 将这些网格单元映 射为相应的平面几何形状(如平面三角形、四边形等) ，所用的插值基函数相同， 并且对于单元节点处的 及其法线方向导数值， 也可以相同的一组插值基函数为 基函数、节点处的 及其法线方向导数值为参数对单元上任一点处的 及其法线 方向导数进行插值，即 ( , ) Ni( , ) i ( , )Ni( , ) (2.34 ) i n i n i Ni ( ， )与单元的形状密切相关，故它们也被称为单元的形函数(shape function )。 上面对三角形和四边形单元选取的节点个数分别为3、 4，即选取它们的所 有顶点作为节点，对应的形函数为线性形函数。为了达到较高的插值精度， 可

28、以 对单元选取更多的节点， 如对三角形单元选取 6个节点，对四边形单元选取 8 个 A 5 2 A x 图 2-9 等参数单元的映射关系 a4 a8 a3 a5 a7 a1 a6 a2 节点，如图 2-9 所示, 对应的形函数为二次形函数 ,对于三角形单元 : 西北工业大学硕士论文 第二章 分布式目标回波分析方法 N1 (2 1) , N 2 (2 1) , N 3 (2 1) , N 4 4 ,N 5 4 , N 6 4 (2.35) 其中, =1-。对于四边形单元 : N1(1)(1)(1)/ 4，N2(1)(1)(1)/4， N3(1)(1)(1)/4， N4(1)(1)(1)/4 N5

29、(12)(1)/2 ，N6(12)(1)/ 2， N7(12)(1)/2 ，N8(12)(1)/2 (2.36 ) 对于一个空间三维目标， 我们可以将其表面边界划分为若干个曲面三角形和 四边形单元，结合式 (2.35 )、式 (2.36 )提供的形函数，将式 (2.34 )代入边界积 分方程式 (2.23 )，并且把边界积分化为在各表面单元上积分之和，设节点数、 单 元数分别为 向量： K、 M，对应于每一个节点 Pn，在每个单元上将会得到两个单元系数 h nj N dS n 1,2, ,K Sj n (2.37 )其 g nj N dS j 1,2, ,M S 中， N是单元形函数组成的向量

30、，式 (2.23 )被化成如下形式： 式中 MM C(Pn) (Pn)j h nj 4 i (Pn)j g nj (2.38) j 1 j 1 j ( j1, j 2, , j ) , 西北工业大学硕士论文 第二章 分布式目标回波分析方法 j ( j1, j 2, , j ) ( , n j1 n j2 为各个单元对应的节点数 ,如三角形为 6,四边形为 8(二次插值情况 ), 按每 个节点 Pm , 将方程(2.38) 的系数合并 ,得到: H nm h nj C(Pm ) h nj mn (2.39a) Gnm g nj nj (2.39b) HG 4i ( (P1), (P2), , (

31、PK )T n PK 于是就得到 K个方程组成的线性方程组 , 写成矩阵形式 : (2.40) 其中 ii(P1), i(P2), , i(PK) T 对于一个边值问题 , 将式(2.40) 通过矩阵分块、移项等变换将已知量移到方程右 边，未知量移到方程左边，可得到： A X B (2.41) 解此方程即可得到单元节点上未知的速度势 及其法线方向导数值。 进而可由式 (2.23) 求出介质空间中任一点处的速度势 (C(P)=4 )。 4分布源边界点法( DSBP)M 前面所介绍的边界元方法理论上可以求解任意形状几何体声散射的准确解， 所以，由于其应用灵活而被广泛地使用。但是， 当波数为目标内部

32、问题的特征值 (目标的谐振频率)时，将产生非唯一解或奇异解， 尽管人们采用了混合积分方 西北工业大学硕士论文 第二章 分布式目标回波分析方法 程法（ CHIEF）来解决这类问题，但还是避免不了对变量的插值以及对奇异积分 的处理。张胜勇提出一种称为分布源边界点（ DSBP）M的方法 47 ，并将其应用到 对振动体声辐射的研究中。 这种方法是通过在振动体边界节点法线方向上 （背离 分析域）施加特解源，构作出一系列系统方程的特解， 从而对偶地表达出系数矩 阵。 本文主要考虑的是对目标的主动探测， 应当着重研究目标的散射声场。 既然 散射可以看作由入射波在目标表面激发的次级声源的二次辐射， 与辐射过程

33、有着 密切的联系，边界点法在辐射和散射的情况下就有类似的机理。基于这一考虑， 我们将边界点法应用到对目标散射回波的分析中， 下面将详细介绍这一方法的实 现过程。 如图 2-10 所示,S 为目标表面边界，将 S进行网格剖分， p、q 为 S上全部 M 个节点中任两个节点。 在 p 节点的法线方向、 距 p 点一定距离处在 S 内构作一个 不与 S面交割的有限大小的面 （例如半径为 R的圆片、球缺面等等） ，并在 面上施加分布源， 被称为加态面。 在此采用相对于 面的投影面为单位均匀分 布，于是就可以求得其在 q 点的解： 2 Rsin T （q, p） 0 d 0 d T 2 Rsin （2.

34、42） nqT（q,p） 0 d 0nq d 式中, 1 ej r/c 4r 1 j r r c nq r r , r q , 1(r nq) 西北工业大学硕士论文 第二章 分布式目标回波分析方法 图 2-10 分布源边界点法原理图 为源面上的变点， nq为节点 q处的内法线向量。式 (2.42 )可以通过积分上 下限的变换 , 利用标准的二维高斯数值求积来确定。 对于散射问题，若只考虑散射波，则式 (2.23 )所给出的边界积分方程可以 表示为： C(P) s(P) S (P,Q) ns (Q) s(Q) n (P,Q)dS (2.43 ) 其中，s 为散射波速度势。通过表面边界网格化和变量

35、插值，可将上式化为如 下矩阵形式： (2.44 当 q 遍取 S面上的所有节点时，由式 (2.42 )形成的 M维列向量 T(p) 、 T (p) 就是方程 (2.44 )的一个特解，因而满足： n A T(p) B (p) (2.45 ) 同理，当 p 遍取边界 S 上全部节点时，就形成 M个特解构成的矩阵 T 和 T n ， 西北工业大学硕士论文 第二章 分布式目标回波分析方法 同样也满足方程 (2.44 )： A T B T (2.46 ) Tn 由此可以得出对偶关系： 1 A 1 B T T (2.47 )因 Tn 此，在散射体边界节点散射速度势法向导数 s n 已知的情况下，节点上的

36、散 射速度势为： s A 1 B sT Ts(2.48 ) n nn 由式(2.43 )即可求出介质空间中任一点 P处的散射速度势 s(P) 。 分布源边界点法还可以避免边界量的解算， 通过下式直接求得声场中任意一 点 P 处的速度势 s(P) (2.49) 其中 T (P) 是上述 M个分布源在 P点产生的特解构成的 M维行向量。 这使 得求解散射声场的运算量得到很大的降低。 分布源边界点法避免了奇异积分的求解， 对于平面边界， 可以选用圆片或矩 形片等作为加态面 ；对于远离 的节点 q，可近似地用位于 源心的点代替分 布源求特解，从而进一步提高计算速度。从一些算例中发现，加态面 与节点 p

37、 之间的距离、 加态面的曲率半径 R以及夹角 具有较大的选择范围， 这使得该 算法具有较高的稳定性。 5时域积分方程法( TDIEM) 西北工业大学硕士论文 第二章 分布式目标回波分析方法 时域积分方程法 (TDIEM)又叫 Kirchhoff 推迟积分法，通过变化时间步长直 接在时域用数值积分法解积分方程。 它是一种较精确的方法， 计算区域限制在结 构体表面，且可应用于任何几何形状的散射体，其数值积分是依次进行、 而不是 同时进行的。 Friedman 和 Shaw18 用该方法研究了平面波入射下有任意截面的圆 柱障碍物的二维声散射问题。 Mitzner 15 发展了任意形状刚性表面瞬态散射

38、数值 解。本文将这种算法与潜艇目标的表面网格化结构相结合， 建立了一种潜艇的连 续型散射模型， 并通过仿真对潜艇的回波结构进行了分析， 这些内容将在第三章 详细讨论。这里先对其理论基础作简要介绍。 图 2-11 散射问题坐标系 取如图 2-11 所示的坐标系， 则空间任一点的声压可由下式表示： p(r,t)pinc(r,t)1 1p(r,)(nR)p(r,)dS(2.50 ) inc 4 S R n 其中， pinc (r, t)为入射声压，t R/c，R |r r |，R (r r )/ R， 是算 子(1/ R2) (1/ Rc)( / )。 当 r 在物体表面边界 S 上时，有： p(r

39、,t)2pinc(r,t)1 1p(r,)(nR)p(r,) dS( 2.51 ) 2 S R n 对于刚性边界，由于 p(r , )/ n 0 ，所以有： p(r,t) 2pinc (r,t) 1 (n,R) p(r , )dS (2.52) inc 2 S 将表面 S网格化后，为了避免奇异性，先将包含 r 点本身的项从积分中分离出来， 西北工业大学硕士论文 第二章 分布式目标回波分析方法 而余下各项的积分由各单元上的求和来代替， 并假设每个单元上的 R 与声压为常 数，然后对包含 r 点自身区域进行解析积分，最后可得到离散形式： (1 i)p(ri,tj) 2pinc(ri,tj)Sl (

40、n R) p(rl, ijl ) (2.53) l i 2 其中 i i ki1 ki 2 ， iSi ， 4 ki1 、ki2 为主曲率。通过对上面微分方程进行数值求解， 即可求得表面上各单元节 点处的声压，再由离散化的数值计算公式 p(r,tj) pinc(r,tj)Sl(n R) p(rl, jl) (2.54) l4 即可求出空间任一点 r 处的声压。 2.3 计算机仿真 本节将运用所介绍的几种方法对形状规则的简单几何体的散射情况进行讨 论，以便通过对比来研究方法的性能。 (1) 用 BEM求解脉动球的表面声压 设一球体的半径为 R，径向振动速度为 VR，理论分析给出在距球心 r 处的

41、声 压可表示为： p(r) RVRZ0 jkR exp jk(r R) (2.55) r 1 jkR 其中，Z0 0 c为传播媒介的声阻抗 ( 0, c分别为传播媒介的密度及声速 )，k 为 波数。 为了在不同波数情况下对表面声压进行比较， 假设参数 R、VR及 Z0 保持不变。 考虑到球体的对称性， 我们仅对第一卦限的部分进行网格化。 在仿真中使用两种 网格划分方法， 一种将表面划分为两个三角形单元及两个四边形单元， 选取节点 数 共 19 个 ；另 一种将 表 面划 分为 两个 三角 形 单元 及六 个四 边形 单 元， 西北工业大学硕士论文 第二章 分布式目标回波分析方法 图 2-12

42、球面的网格化示意图选 取 节 点数共 39 个，网格化的结果如图 2-12 所示。采用 BEM中的二次形函数对变量插 值，并将所得 BEM结果与理论解进行比较。仿真的结果如表 2-1 所示。 kR 理论解 19 个节点 相对误差 (%) 39 个节点 相对误差 (%) 1 0.5000 +j 0.5000 0.5014+j0.4994 0.28+j0.12 0.5007+j0.4994 0.14+j 0.12 2 0.8000+j0.4000 0.8013+j0.3979 0.16+j0.53 0.8009+j0.4016 0.11+j0.40 3 0.9000+j0.3000 0.8986+

43、j0.3489 0.16+j16.3 0.8986+j0.3091 0.16+j 3.03 3.14 0.9079+j0.2891 0.9179+j0.3794 1.10+j31.25 0.9170+j0.3191 1.00+j10.41 4 0.9412+j0.2353 0.9401+j0.2287 0.12+j2.80 0.9414+j0.2304 0.02+j 2.08 5 0.9615+j0.1923 0.9615+j0.1952 0.00+j1.51 0.9615+j0.1946 0.00+j 1.20 表 2-1 BEM 结果与理论解的对比 由表中所列结果可见，边界元方法在求解表面

44、声压时具有较高的精确度，网 格划分越细致，运算结果的误差就越小。在 kR=3 及 3.14 时误差有所增大，这 是由于边界元算法本身在物体谐振频率上会产生较大误差所致， 这也是该算法的 不足之处。 (2) 用边界点方法求解脉动球的表面声压 对上面所设定的脉动球用分布源边界点和一般边界元两种方法来求解其表 面声压，边界点法采用的整个球面总节点数为 20，边界元法采用图 2-12 中的 39 节点的网格划分方式。在表 2-2 中给出了运算结果。 西北工业大学硕士论文 第二章 分布式目标回波分析方法 kR 理论解 边界元法 相对误差 (%) 边界点法 相对误差 (%) 1 0.5000+j0.5000 0.5007+j0.4994 0.14+j0.12 0.5000+j0.5000 0.00+j0.00 2 0.8000+j0.4000 0.8009+j0.4016 0.11+j0.40 0.7995+j0.3999 0.06+j0.03 3 0.9000+j0.3000 0.8986+j0.3091 0.16+j3.03 0.8995+j0.