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

西北工业大学硕士论文第二章 分布式目标回波分析方法第二章 分布式目标回波分析方法当声波在传播中遇到目标时，会在目标表面激发起次级声源，它们向周围介质中辐射次级声波，我们把这些次级声波统称为散射波。在收发合置的情况下，返回声源方向的那部分散射波被称为目标回波。目标回波与目标的散射特征密切相关，它是入射波与目标相互作用后才产生的，在此过程中，有关目标本身的某些特征信息会包含在回波中，人们通过对回波的分析处理，可以将目标的特征信息提取出来，再辅以某些先验知识，从而实现对目标的探测。例如，通过接收基阵采集目标对某种入射激励信号的回波，再对阵列的快拍数据进行分析处理，就可以估计出目标的个数、方位、距离、速度等参数。由此可见，研究目标回波的特性，在工程上是很有意义的。目标回波不仅与目标本身的各种特征有关，如形状、大小、材料、结构等，还与入射波的形式有关，入射波时宽、带宽将直接决定回波的结构。宽带入射信号将产生宽带回波，窄带入射信号引起窄带回波，时宽入射信号产生稳态响应，窄脉冲信号产生瞬态响应。本文主要就高频窄带信号在稳态和瞬态两个方面进行介绍。此外，目标回波一般来说包括几何散射波和弹性散射波，几何散射波是由目标本身的外形特征决定的，与目标的材料无关；而弹性散射波是由于目标的弹性声阻抗与介质不同而引起的，与目标本身的材料有关系。本文所关心的主要是目标的几何特征，因此，以后将对几何回波作着重研究。 2.1 离散散射型目标回波所谓离散散射型目标是把目标当作若干个散射很强的亮点组成的的散射体,目标的回波是由所有亮点的子回波的线性叠加。式(1.1）、(1.2）分别给出了单亮点和多亮点目标的传递函数，要得到目标的传递函数，必须先确定式中的三组参数，即幅度因子Am 、 时间延迟m 、相位跳变m（m=1，2，,N），人们通过研究发现了一些几何体的参数。1凸光滑曲面的镜反射 设Pi()为镜反射亮点处入射波的频谱，则由该亮点产生的镜反射波的频谱为： (2.1)其中,R1、R2为镜反射亮点处的目标表面的两个主曲率半径，r为空间散射场中某一点与亮点之间的距离，k为入射波波数，ejkr/r用于修正信道传输的影响，因为无论目标多么复杂，在远场中散射波都以ejkr/r的规律扩展。三个亮点参数分别为： ， =0， =0 (2.2）若满足rR1及rR2，则 (2.3)2 有限长圆柱体L2aa1234图2-1 有限长圆柱体 如图 2-1所示,为声波相对于圆柱体的入射角, L为圆柱体长度，2为直径，取棱角1为参考点，则被照射到的三个亮点的散射参数为： ， ，(2.4a） ， ， (2.4b） ， ，(2.4c）可见，幅度因子和时间延迟是波数k、半径a、圆柱体长度L以及入射角的函数，而相位跳变是一个常数。3无限长圆柱体 根据Sommerfeld-Watson变换分析和实验结果，无限长圆柱回波包括镜反射波和各型表面环绕波组成，在远场高频条件下，回波的主要成分为镜反射回波。镜反射回波可以表示为： (2.5）镜反射亮点反射参数为： , , (2.6)这里参考点取在镜反射点上,M为镜反射系数,当高频(大ka)时可以用平面波正入射到柱体表面的反射系数近似。一旦确定的目标的亮点散射参数，就可以根据式(1.2）模拟出目标的散射回波。 2.2 连续散射型目标回波所谓连续散射型目标是指把分布式目标看作表面各点的散射特性在空间位置上是连续变化的几何体,目标表面各点都对回波有贡献。这种观点更有利于准确地描述目标的散射特性以及精确探测目标的回波结构。由于表面各点的位置、法线向量、声阻抗等参数不同，它们对回波贡献也就不同。对于连续散射型目标回波是通过一些基于面元法的思想来描述和分析的，将目标表面的面积微元对回波的贡献在整个表面进行面积分，就可以获得回波响应。若表面形状很复杂，对表面积分就很繁琐，甚至根本不能得出结果，在这种情况下，可以将表面按某种方式网格化，再运用数值计算的方法就可以得出模拟的回波结果，人们在实际中发现，这些数值方法得出的结果具有足够的精确度，可以满足工程应用的要求。1赫姆霍兹（Helmholtz）积分近似法A r1B dS图2-2 声散射几何关系r2设有一目标置于均匀介质空间中，目标的外表面为封闭曲面S，其外法线向量为n,有一点声源位于点，下面计算介质空间中另一点处的散射声场。如图2-2所示，B点散射声场的赫姆霍兹积分公式为： (2.7）其中，k为波数，s是散射声场的速度势函数，表示在S面上的外法线方向导数, r1、r2分别是A、B点到表面上面元dS的距离。由于被积函数中的s及其法线方向导数是未知量，所以，上式的积分一般不能直接求出。但是，应用边界条件，可将被积函数中的未知量用已知量表示。假设目标表面S是刚性的，则在S上有 (2.8）其中，i是入射波势函数，它可表示为 (2.9）这里A是一个表示振幅的常数（式中略去了时间因子ejt）。由式(2.8）及(2.9），并考虑远场高频条件kr11,我们可以得到 (2.10)其中,表示矢量与法线向量之间的夹角。同样地，在kr21时有 (2.11）作为一种足够精确的近似，可以认为刚性目标表面上的散射声场等于入射声场，于是 (2.12）将式(2.10）、(2.11）、(2.12）代入式(2.7），就可以得到 (2.13）在收发合置的情况下，r1=r2=r，则上式成为 (2.14）式(2.13）、(2.14）即为目标为刚性体时散射声场的高频近似积分公式，只要设法将积分求得，就可以求出介质空间中任一点处的散射声场。式(2.13）的物理意义可以理解为：目标表面各点在入射声波的激励下，作为次级声源辐射出二次级波，接收点上的散射声场就是目标表面上所有二次源辐射的声波在该点的叠加，这些二级次波的相位为k(r1+r2),即取决于声波在介质中的往返波程。2费涅尔(Fresnel)半波带近似法式(2.13)、(2.14）在知道目标表面方程的情况下是可以求得其散射声场的，但对于复杂目标表面，计算过程将变得相当繁琐，因此，工程上常会应用费涅尔半波带近似法来简化计算工作量。基于散射声波的相位取决于入射波在目标表面各点的往返声程这一机理，我们可以仿照光学中的方法，应用费涅尔半波带法来求解式(2.14）。如图2-3所式，设入射波从B点发射，C点为目标表面上距B点最近的点，设它距B点的波程为r0，以B点为球心，r0为半径作球面，它与目标表面相切于C点，接着跳跃式地增加球面半径，每次增加1/4波长， 即r0+/4，r0+/2，r0+/4S3CS4S1S2Br0r0+/2r0+3/4r0+图2-3 菲涅尔半波带示意图r0+3/4，r0+， ，这些半径不同的球面把目标表面分割成许多环带S1，S2，S3，S4， ，SN，它们被称为费涅尔带。可以看出，相邻半波带产生的散射波在B点的声程差为/2，即相位差为。因此，若目标表面上共有N个半波带，且第i个波带的贡献用i表示，则 (2.15）其中， (2.16）若目标的径向尺寸远远大于波长,并且目标表面不太弯曲,那么,目标表面可以分割出许多费涅尔半波带,而且,相邻波带的变化也不大,彼此的面积也很接近,这样就可以认为第i个波带产生的散射波等于相邻两个波带产生散射波的平均值,即 (2.17)将上式代入式(2.15),得到总的散射波为 (2.18)即总的散射波等于第一个与最后一个费涅尔带所产生的散射波的平均值。对于大体积目标，最后一个费涅尔带的，因而它的贡献可以忽略不计，而第一个费涅尔带的，于是，总的散射波可以简化表示为 (2.19）这样，所得到的积分表达式就比式(2.14）容易计算了。3边界元法（BEM）边界积分方程在计算声学中的应用可以追溯到二十世纪六十年代。最早将边界积分方程应用于任意形状物体的声辐射的是Chen和Schweikert4，此后，很多学者进行了积极的探索。在绝大多数的这些早期研究工作中，对物体表面的网格化要么是应用常量的网格单元形函数，要么将网格单元近似地以平面来代替。近年来，随着边界元方法的迅速发展，人们将二次单元形函数应用到对声学问题中边界积分方程的求解，从而使声学物理量以及边界几何特征能够被更精确地描述。SBB图2-4 介质空间中的物体n（1） 自由介质空间中的边界积分公式如图(2.4)所示,设有一物体B位于无限均匀介质空间B中,物体表面为S,介质密度和声速分别为0 、c. 在稳态时，介质空间B中的声场应该满足如下关系： (2.20）这就是著名的赫姆霍兹(Helmholtz)方程，其中，是速度势，k是入射波中心频率波数。对于散射问题，速度势由s及i两部分组成，i为入射波速度势,s为散射波速度势,这两部分在B中都必须满足赫姆霍兹方程。物体边界S上满足的条件具有一般的形式 (2.21）其中，n被定义为物体表面上的内法线向量，、为特定常数。散射波速度势还应当满足萨默菲尔德(Sommerfeld)条件,即 (2.22)上式反映了在无穷远处速度势所满足的关系。应用格林第二公式或加权残留量公式，由式(2.20)可以得出在边界S上的边界积分方程: (2.23)上式被称为赫姆霍兹积分方程。其中，(P)为声场中某一定点P处的速度势，Q为物体边界上的变点，C(P)为一特定的常量，它取决于P点的位置，i(P)为入射波在P点的速度势,(Q,P)为自由空间中的格林函数,它可以表示为: (2.24) 其中，r=|Q-P|，即r为Q、P两点间的距离。值得注意的是系数C(P)的选取，当定点P在介质空间B中时，C(P)=4;当P点位于物体内部时,C(P)=0；当P点位于物体边界上并且在P点有唯一的切平面(即P点位于界面上的光滑部分),则C(P)=2；而当P点位于界面上的棱角处时,则C(P)由下式确定: (2.25) 实际上,上式对于处于任何位置的P点都适用。由于没有有效的数值方法判断P点是否位于棱角上，因此对于有棱角的物体，只要P点在边界上，我们一律用式(2.25）计算C(P)。（2） 存在界面时的边界积分公式在实际中,有很多目标位于界面附近,比如海面舰艇、地面建筑物、靠近海面的潜艇等等，而在这种情况下的边界积分公式与在自由介质空间中的情况有所不同，所以，需要对其进行专门介绍。BSHSPP1rr1BnQ图2-5 界面附近物体声散射几何关系如图2-5所示,设SH是一个无限反射平面,物体边界上各点在定点P的散射波不仅包括动点Q的贡献,还包括Q点关于界面SH的镜像点的贡献,所以,在求解P点的速度势时,式(2.23)的(Q,P)应当替换为 (2.26)H(Q,P)被称为半无限空间格林函数,其中,界面反射系数RH在界面为刚性时等于1,在界面为软界面时等于 1,r=|Q-P|,r1=|Q-P1|。另外，在收发合置的情况下，入射波应当包含直达波和界面反射波，除非入射方向与界面平行。所以，入射波在介质中传输的扩展因子为: (2.27）BBnPP1Qrr1SH图2-6 与界面相接触的物体声散射几何关系SOSC如果物体边界与界面有相接触的部分，如图2-6所示,物体边界可分为两部分,一部分是SC，它与界面SH相接触,另一部分是SO，它与界面不接触，那么，在这种情况下的边界积分方程应当为： (2.28）与式(2.23)不同的是,上式采用的是H，而不是，并且积分是在SO上进行的。当P点在B或B时,系数C(P)的选取与自由介质空间中的选取方法相同,当P点在SO上,C(P)由下式确定: (2.29)当P点位于SC上,则C(P)为: (2.30)（3） 数值计算公式前面所介绍的求解散射声场的积分方程法需要知道物体表面边界上的速度势及其法线方向导数,但在散射声场被求解出之前,处于物体边界上的这些物理量却是未知的,因此,在求解介质空间中的散射声场时,必须先得出物体边界上各点的速度势及其法线方向导数值,再将这些求解出的结果代入边界积分公式,进而求出介质空间中任一点处的散射波速度势。问题的关键在于求解物体边界上的速度势及其法线方向导数。应用数值计算的方法，可以使这个问题得以解决。将物体的表面边界S以某种方式划分为若干小的曲面片，这些曲面片被称为表面边界的单元，在S上再选取若干个节点，我们称这个过程为表面边界的网格化。设表面边界单元数和节点数分别为M、N，令P点分别置于各个节点上，对应于每一个节点，应用边界积分公式就会得到一个方程，这样，当P遍取边界上的N个所有节点，就得到N个线性独立方程组成的方程组，其中的未知量为2N个（N个表面边界速度势及其法线方向导数），结合边界条件中给出的N个已知量，最终将会得到一个N元恰定的线性方程组，通过求解这个方程组，就可以得到边界上各节点处的速度势及其法线方向导数。下面将具体介绍这一过程。在讨论应用单元节点插值方法解边界积分方程之前,有必要先简要介绍一下单元局部坐标的概念。下面以任意三角形和四边形为例进行介绍。如图2-7所示,设xy平面上一三角形单元e的三个顶点坐标为P1(x1,y1)、P2(x2,y2)、P3(x3,y3)，该三角形内任一点的坐标为（x，y），则如下变换：P1xP2P3e2 (0 ,1)1(1 , 0)3 ( 0 ,0 )图2-7 三角形单元的局部坐标映射关系y ， ， (2.31a）或 (2.31b）将三角形单元e 映射为平面中的等腰直角三角形 。其中，e为三角形e的面积。、一般被称为局部坐标，它们与三角形e的形状和大小有关，由于、满足+=1的关系，独立的坐标只有两个，所以该变换并没有改变维数。而且局部坐标与直角坐标x、y的选取方法是无关的，这给处理问题带来了很大方便。同理，对于任意四边形单元,也有类似的转换关系。如图2-8所示,四边形单元e被映射成平面内的中心在原点的单位正方形 ,a1、a2、a3、a4分别对应于四边形的顶点A1、A2、A3、A4，所用的变换为：xyA4(x4,y4)A3(x3,y3)A2(x2,y2)A1(x1,y1)e图2-8 四边形单元的局部坐标映射关系a4(-1 , 1)a3( 1 , 1)a2( 1, -1)a1(-1, -1) (2.32）其中，Ni(,)称为双线性插值基函数,它们的表达式为: (2.33)这样,对于一个单元上的任意一点,总存在一组局部坐标与之对应,反之亦然。即单元上点的全局坐标与局部坐标具有一一对应的关系。实际上，以上变换可以看作以节点(三角形或四边形的顶点)坐标为参数 、Ni（，）为基函数对单元（三角形或四边形）上任一点坐标的插值。对于空间的曲面网格单元，我们可以根据等参数单元的变换的办法，将这些网格单元映射为相应的平面几何形状（如平面三角形、四边形等），所用的插值基函数相同，并且对于单元节点处的及其法线方向导数值，也可以相同的一组插值基函数为基函数、节点处的及其法线方向导数值为参数对单元上任一点处的及其法线方向导数进行插值，即 (2.34）Ni（，）与单元的形状密切相关，故它们也被称为单元的形函数（shape function）。上面对三角形和四边形单元选取的节点个数分别为3、4，即选取它们的所有顶点作为节点，对应的形函数为线性形函数。为了达到较高的插值精度，可以对单元选取更多的节点，如对三角形单元选取6个节点，对四边形单元选取ea5a4a3a2a1 a6yzxA5A1A6A2A3A4zxyA8A7A6A5A4A3A1A2ea8a7a6a5a4a3a2a1 图2-9 等参数单元的映射关系8个节点，如图2-9所示,对应的形函数为二次形函数,对于三角形单元: , , , , , (2.35)其中,=1-。对于四边形单元: ， ， ， ， ， ， (2.36） 对于一个空间三维目标，我们可以将其表面边界划分为若干个曲面三角形和四边形单元，结合式(2.35）、式(2.36）提供的形函数，将式(2.34）代入边界积分方程式(2.23），并且把边界积分化为在各表面单元上积分之和，设节点数、单元数分别为K、M，对应于每一个节点Pn，在每个单元上将会得到两个单元系数向量： (2.37）其中，N是单元形函数组成的向量，式(2.23）被化成如下形式： (2.38)式中 , 为各个单元对应的节点数,如三角形为6,四边形为8(二次插值情况),按每一个节点Pm ,将方程(2.38)的系数合并,得到: (2.39a) (2.39b)于是就得到K个方程组成的线性方程组,写成矩阵形式: (2.40)其中 对于一个边值问题,将式(2.40)通过矩阵分块、移项等变换将已知量移到方程右边，未知量移到方程左边，可得到： (2.41)解此方程即可得到单元节点上未知的速度势及其法线方向导数值。进而可由式(2.23)求出介质空间中任一点处的速度势(C(P)=4)。4分布源边界点法（DSBPM）前面所介绍的边界元方法理论上可以求解任意形状几何体声散射的准确解，所以，由于其应用灵活而被广泛地使用。但是，当波数为目标内部问题的特征值（目标的谐振频率）时，将产生非唯一解或奇异解，尽管人们采用了混合积分方程法（CHIEF）来解决这类问题，但还是避免不了对变量的插值以及对奇异积分的处理。张胜勇提出一种称为分布源边界点（DSBPM）的方法47，并将其应用到对振动体声辐射的研究中。这种方法是通过在振动体边界节点法线方向上（背离分析域）施加特解源，构作出一系列系统方程的特解，从而对偶地表达出系数矩阵。本文主要考虑的是对目标的主动探测，应当着重研究目标的散射声场。既然散射可以看作由入射波在目标表面激发的次级声源的二次辐射，与辐射过程有着密切的联系，边界点法在辐射和散射的情况下就有类似的机理。基于这一考虑，我们将边界点法应用到对目标散射回波的分析中，下面将详细介绍这一方法的实现过程。如图2-10所示,S为目标表面边界，将S进行网格剖分，p、q为S上全部M个节点中任两个节点。在p节点的法线方向、距p点一定距离处在S内构作一个不与S面交割的有限大小的面（例如半径为R的圆片、球缺面等等），并在面上施加分布源，被称为加态面。在此采用相对于面的投影面为单位均匀分布，于是就可以求得其在q点的解： (2.42)式中, , , pz1x1y1Sz2y2x2q图2-10 分布源边界点法原理图R为源面上的变点，为节点q处的内法线向量。式(2.42）可以通过积分上下限的变换,利用标准的二维高斯数值求积来确定。对于散射问题，若只考虑散射波，则式(2.23）所给出的边界积分方程可以表示为： (2.43）其中，s为散射波速度势。通过表面边界网格化和变量插值，可将上式化为如下矩阵形式： (2.44）当q遍取S面上的所有节点时，由式(2.42）形成的M维列向量、就是方程(2.44）的一个特解，因而满足： (2.45）同理，当p遍取边界S上全部节点时，就形成M个特解构成的矩阵和，同样也满足方程(2.44）： (2.46）由此可以得出对偶关系： (2.47）因此，在散射体边界节点散射速度势法向导数已知的情况下，节点上的散射速度势为： (2.48）由式(2.43）即可求出介质空间中任一点P处的散射速度势s(P)。分布源边界点法还可以避免边界量的解算，通过下式直接求得声场中任意一点P处的速度势s(P) (2.49)其中 是上述M个分布源在P点产生的特解构成的M维行向量。这使得求解散射声场的运算量得到很大的降低。分布源边界点法避免了奇异积分的求解，对于平面边界，可以选用圆片或矩形片等作为加态面；对于远离的节点q，可近似地用位于源心的点代替分布源求特解，从而进一步提高计算速度。从一些算例中发现，加态面与节点p之间的距离、加态面的曲率半径R以及夹角具有较大的选择范围，这使得该算法具有较高的稳定性。5时域积分方程法（TDIEM）ORSP(r,t)rrn图2-11 散射问题坐标系时域积分方程法(TDIEM)又叫Kirchhoff推迟积分法，通过变化时间步长直接在时域用数值积分法解积分方程。它是一种较精确的方法，计算区域限制在结构体表面，且可应用于任何几何形状的散射体，其数值积分是依次进行、而不是同时进行的。Friedman和Shaw18用该方法研究了平面波入射下有任意截面的圆柱障碍物的二维声散射问题。Mitzner15发展了任意形状刚性表面瞬态散射数值解。本文将这种算法与潜艇目标的表面网格化结构相结合，建立了一种潜艇的连续型散射模型，并通过仿真对潜艇的回波结构进行了分析，这些内容将在第三章详细讨论。这里先对其理论基础作简要介绍。取如图2-11所示的坐标系，则空间任一点的声压可由下式表示： （2.50）其中， 为入射声压，是算子。当在物体表面边界S上时，有： （2.51）对于刚性边界，由于 ，所以有： (2.52)将表面S网格化后，为了避免奇异性，先将包含r点本身的项从积分中分离出来，而余下各项的积分由各单元上的求和来代替，并假设每个单元上的R与声压为常数，然后对包含r点自身区域进行解析积分，最后可得到离散形式： (2.53)其中 ，ki1、ki2为主曲率。通过对上面微分方程进行数值求解，即可求得表面上各单元节点处的声压，再由离散化的数值计算公式 (2.54)即可求出空间任一点处的声压。 2.3 计算机仿真本节将运用所介绍的几种方法对形状规则的简单几何体的散射情况进行讨论，以便通过对比来研究方法的性能。(1) 用BEM求解脉动球的表面声压设一球体的半径为R，径向振动速度为VR，理论分析给出在距球心r处的声压可表示为： (2.55)其中，为传播媒介的声阻抗(分别为传播媒介的密度及声速)，k为波数。为了在不同波数情况下对表面声压进行比较，假设参数R、VR及Z0保持不变。考虑到球体的对称性，我们仅对第一卦限的部分进行网格化。在仿真中使用两种网格划分方法，一种将表面划分为两个三角形单元及两个四边形单元，选取节点数共19个；另一种将表面划分为两个三角形单元及六个四边形单元，xyzxyz19个节点39个节点图2-12 球面的网格化示意图选取节点数共39个，网格化的结果如图2-12所示。采用BEM中的二次形函数对变量插值，并将所得BEM结果与理论解进行比较。仿真的结果如表2-1所示。kR理论解19个节点相对误差(%)39个节点相对误差(%)10.5000+j0.50000.5014+j0.49940.28+j0.120.5007+j0.49940.14+j0.1220.8000+j0.40000.8013+j0.39790.16+j0.530.8009+j0.40160.11+j0.4030.9000+j0.30000.8986+j0.34890.16+j16.30.8986+j0.30910.16+j3.033.140.9079+j0.28910.9179+j0.37941.10+j31.250.9170+j0.31911.00+j10.4140.9412+j0.23530.9401+j0.22870.12+j2.800.9414+j0.23040.02+j2.0850.9615+j0.19230.9615+j0.19520.00+j1.510.9615+j0.19460.00+j1.20表2-1 BEM结果与理论解的对比由表中所列结果可见，边界元方法在求解表面声压时具有较高的精确度，网格划分越细致，运算结果的误差就越小。在kR=3及3.14 时误差有所增大，这是由于边界元算法本身在物体谐振频率上会产生较大误差所致，这也是该算法的不足之处。(2) 用边界点方法求解脉动球的表面声压对上面所设定的脉动球用分布源边界点和一般边界元两种方法来求解其表面声压，边界点法采用的整个球面总节点数为20，边界元法采用图2-12中的39节点的网格划分方式。在表2-2中给出了运算结果。kR理论解边界元法相对误差(%)边界点法相对误差(%)10.5000+j0.50000.5007+j0.49940.14+j0.120.5000+j0.50000.00+j0.0020.8000+j0.40000.8009+j0.40160.11+j0.400.7995+j0.39990.06+j0.0330.9000+j0.30000.8986+j0.30910.16+j3.030.8995+j0.29990.06+j0.033.140.9079+j0.28910.9170+j0.31911.00+j10.410.9075+j0.28890.04+j0.0740.9412+j0.23530.9414+j0.23040.02+j2.080.9407+j0.23520.05+j0.0450.9615+j0.19230.9615+j0.19460.00+j1.200.9611+j0.19220.04+j0.0560.9730+j0.16220.9745+j0.16710.15+j3.010.9725+j0.16210.05+j0.066.280.9753+j0.15530.9852+j0.17241.02+j11.020.9748+