专利名称：快速预报声呐目标强度的方法
技术领域：
本发明涉及一种声学目标散射特性预报方法，是一种快速预报水声目标强度(target strength，TS)的方法。
背景技术：
本发明以波的传播理论为基础，适用范围为水声工程中的水中声目标强度特性预报和工程计算。
对水下目标回波特性进行预报是水声工程的重要课题，在声呐设计、潜艇和水下兵器的隐身和反隐身设计、精确声制导、水下目标识别、海洋开发等应用中具有重要意义。回波特性研究的核心问题是目标对声波的散射，该问题的解一般通过解析方法或数值方法得到。严格的解析方法具有理论指导意义，但仅能处理一些简单形状的目标，对于形状和材料复杂的工程实际问题通常无能为力。为此，在长期的研究中提出了各种渐近方法如GTD，以及数值方法如有限元/边界元方法(FEM/BEM)、时域有限差分法(FDTD)等。随着计算机技术的发展，数值方法受到高度重视并得到了广泛的应用。
例如，应用FDTD对收发合置(monostatic)条件下水中目标在不同方位角的目标强度计算过程如下1.在已知目标周围确定计算域，对计算域进行网格划分，网格尺度不应大于波长的十分之一，为了得到较高的计算精度，网格通常更小。
2.根据目标声学特性得到适用于FDTD计算的边界条件表达形式。应用FDTD迭代公式和计算域边缘上应满足的吸收边界条件(ABC)，计算目标周围的近场分布。
3.通过近场-远场变换求出远场的后向散射声强度。
4.对目标每一方位重复上述计算，得到目标强度随目标所取方位的变化规律。
由此可见，要计算目标强度必须进行大量的迭代运算以得到计算域的场分布，计算量巨大。在计算中为了保存大量的声学量和辅助变量，对存储器的容量要求也很高。这一缺点在高频、大目标应用中尤为突出，即使采用大规模的并行计算，也难以满足工程应用的需要。
突破现有框架，寻求新的高效数值计算方法是水声工程中的一项迫切课题。
发明内容
本发明的目的在于解决现有技术中存在的迫切课题，提出一种快速预报声呐目标强度的方法。根据已知的目标几何形状、声学性质、入射声波频率和入射方向，利用Fourier衍射定理快速计算水中声目标的目标强度；为达到上述目的，我们作如下研究正问题和反问题已知目标几何形状和声学性质并给定入射波频率和入射方向，计算声压分布、目标强度、指向性等声场特性，是水声学中典型的正问题。反之，由测量所得声场分布重建产生这一声场的散射体几何形状和声学参数分布是一个逆问题。本发明要解决的目标强度预报是声学正问题，它建立在对下列反问题解法即计算机层析成像技术(computerized tomography，CT)的逆向思维基础上。
衍射CT和Fourier衍射定理医学中X射线CT的数学基础是Radon变换。在算法实现中，X射线投影数据与目标断面的二维Fourier变换由Fourier切片定理(Fourier slice theorem)相联系。在层析成像中用波长有限的声波取代X射线照射目标，则必须考虑衍射效应，Fourier切片定理不再成立，而要用Fourier衍射定理(Fourier diffraction theorem)。利用声场测量数据中携带的目标信息反演物体截面的声学参数分布(声成像)称为衍射CT。
考虑单频声波入射情况，Fourier衍射定理可通过下列Helmholtz方程的解导出(▿2+k02)ps(r)=-k02o(r)[pi(r)+ps(r)]---(1)]]>pi(r)是入射场声压，ps(r)是散射场声压，二者之和为总声压。k0为波数(或空间角频率，即单位长度内相位改变的弧度数)。散射体函数o(r)包含了目标的几何形状和声学参数分布。在弱散射假设下方程(1)的解可表示为ps(r)＝∫g(r-r′)o(r′)pi(r′)dr′，ps＜＜pi(2)其中g(r-r′)是Green函数。对上式两端作Fourier变换，经推导可得到以下定理。
Fourier衍射定理图1中物体o(r)受方向为φ的单频率平面波照射，沿直线T1T1′测得的散射场分布称为前向投影。前向投影的一维Fourier变换给出物体二维Fourier变换O(k)在频域中粗线所示低频半圆弧上的取值。圆弧半径等于入射波空间角频率k0。类似地，在入射波同侧沿直线T2T2′测得的散射场称为后向投影，其一维Fourier变换给出O(k)在频域中虚线所示高频半圆弧上的取值。
改变入射波方向可得到许多投影，若样本充分覆盖频域就能通过二维Fourier反变换得到o(r)，于是重建了目标图像。图2给出间隔45°的8个投影对应的频域样本分布实例，为了使图示清楚，每个投影仅取样10点。在实际使用时应提高样本密度。
图2(a)和(b)分别对应于透射型和反射型衍射CT。
解决正问题的新方法Fourier衍射定理的逆向应用先考虑二维问题将无限长柱形目标置于无界海水中，用平面波照射，入射方向与柱轴垂直。取目标周围的矩形区域为计算域。将区域内目标和海水的ρc分布按正方形网格离散化，由此得到代表计算域内ρc分布的广义数字图像，如图3(a)所示。ρ和c分别为介质密度和声速。计算该图像的二维离散Fourier变换可得到图3(b)中的频谱。这里为了图示方便取了较大的象素尺寸(参看以下表1)，因此在高频区呈现不规则性，实际使用时提高空域网格密度可克服。
由于目标图像的空域和频域表示完全等价，因此在二维DFT(频谱)中包含了目标的全部几何和声学参数信息。
根据Fourier衍射定理，在给定入射波频率的条件下，空域内沿入射波方向在目标后(前)侧一条直线上的均匀取样位置对应于频域内相应方向上某一高频(低频)半圆弧上的非均匀取样位置。空域中该直线上的声场分布是相应半圆弧上非均匀分布的二维DFT值的一维离散Fourier反变换。
由此得到目标强度预报的新方法根据已知入射平面波频率、入射方向、目标几何形状、声学参数分布，对目标二维DFT在高频半圆弧上的非均匀分布的频域值作一维IDFT得到空域投影，由后向投影的中点可求得目标强度。目标二维DFT域中低频半圆弧对应的前向投影中点则给出目标前向散射声压。
推广到三维问题计算三维目标广义图像的三维DFT，空域中三维目标受平面波照射，其后(前)侧平面上的投影对应于三维DFT域中相应方向上高频(低频)半球面上的非均匀取样。目标强度可通过二维IDFT计算得到。
根据上述的研究成果，本发明采用下述技术方案一种快速预报声呐目标强度的方法，其特征在于将目标和周围介质的声学参数看成一个平面或立体分布，即广义图像，计算其二维或三维DFT，得到目标的几何形状和声学的频域表示；根据Fourier衍射定理，给定的入射平面波频率和入射方向可确定频域中的一个圆或球面；分别取位于低频侧和高频侧两个半圆弧或半球面上的频域样本，计算它们的一维或二维IDFT，求得前向和后向散射强度。由后向散射特性得到声呐目标强度。
上述方法的具体步骤如下a.在二维情况下，将目标置于无界声场中，用平面波照射，取目标周围的正方形区域为计算域；b.将区域内目标和周围介质的ρc分布按正方形网格离散化，得到代表计算域内ρc分布的图像；c.计算上述图像的二维离散Fourier变换，得到目标的频域表示；d.根据入射平面波的频率和入射方向确定频域内非均匀分布的样本，这些样本位于以入射波波数为半径的高频或低频半圆弧上；e.对频域内的样本作一维IDFT，得到目标空间的一个投影；f.由投影中点及其近邻的平均声压计算目标强度或前向散射声压；g.在三维情况下，步骤a改为取正立方体为计算域；步骤b改为按正立方体将目标离散化；步骤c改为计算三维DFT；步骤d所指非均匀分布样本位于以入射波波数为半径的高频或低频半球面上；步骤e改为对频域内的样本作二维IDFT，得到目标空间的投影；由投影中点及其近邻的平均声压可计算目标强度或前向散射波声压。
本发明与现有技术相比较，具有如下的特点和显著优点本发明基于衍射CT成像的逆向思维提出了一种新的目标散射特性预报方法，可用于水声工程计算目标强度。其计算复杂度大大低于现有的其他数值方法，主要原因如下1.目标广义图像的二维(三维)DFT包含目标的全部几何和声学性质信息；频域取样位置则包含入射波的全部信息。该方法的核心是根据Fourier衍射定理所揭示的这一事实对目标DFT进行简单的取样操作，以此取代解波动方程的复杂过程，避免了一般数值方法中耗费大量时间和内存的迭代运算。
2.仅需计算所要求的声场量即后向(或前向)散射波参数，而不必首先计算目标周围计算域内的整个场分布。
3.充分利用快速Fourier变换实现了高效数值计算。
图1是Fourier衍射投影定理示意图。
图2是空域8个投影对应的频域样本分布图。
图3是图像的空域和频域示图。
图4是对应于前向和后向投影的频域取样点分布图。
图5是钢质椭圆目标后向声散射特性曲线图。
具体实施例方式本发明的一个优选实施例是对水下目标快速预报其水声目标强度。以二维情况为例说明实施本发明提出的声纳目标强度预报方法。步骤如下1.将无限长柱形目标置于无界海水中，用平面波照射，入射方向与柱轴垂直，构成一个二维问题。取目标周围的正方形区域为计算域。
2.将区域内目标和目标周围介质(如海水)的ρc分布按正方形网格离散化，得到代表计算域内ρc分布的图像。
3.计算上述图像的二维离散Fourier变换，得到目标的频域表示(即频谱)。
4.根据入射平面波的频率和方向确定频域内非均匀分布的样本(参看图1)。
5.对频域内的样本作一维IDFT，得到目标空间的投影。
6.根据投影中点及其近邻的平均声压计算目标强度或前向散射声压。取中点近邻平均是为了进行空间平滑，平滑窗口大小可根据实际情况确定。
发明人已用本技术计算了海水中钢质椭圆目标的前向和后向声散射特性。表1给出入射声波和海水声学参数、目标材料和几何参数、数值计算中的采样参数。入射波沿正x轴方向。图4给出频域取样点的位置，圆弧半径为u0＝k0/2π＝1/λ＝1m-1，(a)是前向散射情况，(b)是后向散射。计算中取L＝8，即目标图像的象素大小为0.125×0.125m2。
图5为钢质椭圆目标的后向散射特性预报结果，在计算中将空域采样密度提高到L＝15，旋转目标使方位角θ由0°～360°变化，每隔1°计算后向散射声强度。其余参数见表1。(a)为目标处于正横方向即θ＝0°的后向投影；(b)为目标声强度随方位角θ的变化，其中实线是采用本发明的方法所得结果。作为比较，图中用点划线给出FDTD计算结果(网格尺寸同样是L＝15)。
本发明在计算复杂度方面有巨大优势。以上计算在Pentium IV，2.79GHz，内存512MB的计算机上进行，采用MATLAB R13。FDTD计算耗时达10,538s，程序中使用变量占内存35MB以上；而用本发明的方法仅需7s，程序使用变量占内存2.8MB。
表1 目标和入射波参数名称符号数值单位海水密度ρ 1026kg/m3海水中声速 c 1500m/s钢的密度ρ27800kg/m3钢中声速c25653m/s钢的杨氏模量E 19.5×1010Newton/m2钢的泊松比 σ 0.28-椭圆长轴a 2.3 m椭圆短轴b 0.6 m短轴与x轴夹角 θ 50 °入射波频率 f 1500Hz波长λ＝c/f 1 m每一波长象素数 L 8 -空域采样间隔δ＝λ/L1/8 m波数(空间角频率)k0＝ω/c 6.28radians/m空间频率u0＝f/c1 m-1计算域大小 Ax＝Ay15 m每一投影的样本数N＝[Ax/δ]+1 121 -频域采样间隔f＝1/Ax1/15m-权利要求
1.一种快速预报声呐目标强度的方法，其特征在于将目标和周围介质的声学参数看成一个平面或立体分布，即广义图像，计算其二维或三维DFT，得到目标的几何形状以及声学的频域表示；根据Fourier衍射定理，给定的入射平面波频率和入射方向确定频域中的一个圆或球面；分别取位于低频侧和高频侧两个半圆弧或半球面上的频域样本，计算它们的一维或二维IDFT，求得前向和后向散射强度，由后向散射特性得到声呐目标强度。
2.根据权利要求1所述的快速预报声呐目标强度的方法，其特征在于具体步骤如下a.在二维情况下，将目标置于无界声场中，用平面波照射，取目标周围的正方形区域为计算域；b.将区域内目标和周围介质的ρc分布按正方形网格离散化，得到代表计算域内ρc分布的图像；c.计算上述图像的二维离散Fourier变换，得到目标的频域表示；d.根据入射平面波的频率和入射方向确定频域内非均匀分布的样本，这些样本位于以入射波波数为半径的高频或低频半圆弧上；e.对频域内的样本作一维IDFT，得到目标空间的一个投影；f.由投影中点及其近邻的平均声压计算目标强度或前向散射声压；g.在三维情况下，步骤a改为取正立方体为计算域；步骤b改为按正立方体将目标离散化；步骤c改为计算三维DFT；步骤d所指非均匀分布样本位于以入射波波数为半径的高频或低频半球面上；步骤e改为对频域内的样本作二维IDFT，得到目标空间的投影；由投影中点及其近邻的平均声压可计算目标强度或前向散射声压。
全文摘要
本发明涉及一种目标散射特性的快速预报方法。它是将目标和周围介质的声学参数看成一个平面或立体分布，即广义图像，计算其二维或三维DFT，得到目标几何形状和场特性的频域表示。根据Fourier衍射定理，由给定的入射平面波频率和入射方向确定频域中的一个圆或球面，分别取位于低频侧和高频侧两个半圆弧(半球面)上的频域样本，计算它们的一维或二维IDFT，求得前向和后向散射强度。该方法的核心是对代表目标和介质的广义图像频谱进行简单的取样操作，以此取代解波动方程的复杂过程，从而避免一般数值方法中耗费大量时间和内存的迭代运算。此外这一新方法也不必先求出计算域内的场分布，而且充分利用快速Fourier变换完成主要的数值计算。该方法大大降低了目标散射特性预报的计算复杂度，特别适合于工程实际应用。本发明可用于水声工程，进行声纳目标强度快速预报。
文档编号G01S15/88GK1605885SQ200410084278
公开日2005年4月13日 申请日期2004年11月18日 优先权日2004年11月18日
发明者王朔中, 冯玉田, 陆健峰 申请人:上海大学