CN104778754B - 水下目标动态尾迹与流体微元轨道速度的数值仿真方法 - Google Patents

Abstract

本发明属于海洋遥感监测技术领域，具体为风驱动态海洋环境中水下目标动态尾迹及海面流体轨道运动的数值仿真方法。其步骤为：将海面剖分网格单元；在每个时刻、每个网格节点处用线性波理论和海浪谱计算风驱起伏海面的波高和轨道速度矢量。在初始时刻，由水下运动目标的静态尾迹和运动方向，用最小二乘法将该静态尾迹分解为一系列不同频率、不同角度传播的平面前进波的线性叠加；累计各网格节点处所有平面前进波的轨道速度矢量，即可重建水下目标动态尾迹的轨道速度矢量。该方法能快速高效地数值仿真模拟风驱动态海洋环境中水下运动目标的尾迹及海面流体的轨道运动随时间的变化。

Description

水下目标动态尾迹与流体微元轨道速度的数值仿真方法

发明领域

本发明属于海洋遥感监测技术领域，具体涉及风驱动态海洋环境中水下目标动态尾迹及海面流体微元轨道速度的数值仿真方法。

背景技术

水下目标动态尾迹与流体微元轨道速度的数值仿真是海洋流体力学和海洋遥感监测领域的一个难题。合成孔径雷达(SAR)应用于海洋遥感面临复杂的理论建模挑战，由于雷达成像孔径合成需要平台运动一段距离，积累一定时间(孔径合成时间)，期间雷达平台、海面波、水下目标的尾迹始终处于动态变化过程，该动态过程的准确描述对海洋流体力学建模和海洋SAR图像的理解研究甚为重要^[1]。

发明内容

本发明的目的在于提供一种快速、高效、通用的水下目标动态尾迹及流体微元轨道速度的数值仿真方法。

本发明提出的水下目标动态尾迹及流体微元轨道速度的数值仿真方法，是采用最小二乘法结合线性波理论计算水下目标动态尾迹及海面流体微元的轨道速度，其基本思路是：用最小二乘法将静态尾迹分解为一系列不同频率、不同角度传播的平面前进波；累计各网格节点处所有平面前进波的轨道速度矢量，即为水下运动目标尾迹的轨道速度矢量。

该方法也适用于水面舰船动态尾迹的轨道运动仿真计算，无论尾迹是由经典速度势函数得来还是由精确数值方法计算得来，并且计算精度能够通过平面前进波的数目加以控制。

本发明提出的水下目标动态尾迹及流体微元轨道速度的数值仿真方法，具体步骤为：

(1)建立水下目标尾迹的动态几何模型，将模拟海面离散化，具体是剖分为网格单元，每个网格节点的波高随时间变化，每个网格节点的流体微元轨道速度随时间变化；用点源速度势函数模型(例如Rankine卵形体模型)或计算流体力学软件建立水下运动目标的静态尾迹(即某一固定时刻的波高起伏分布)几何模型。一般来说，网格单元为均匀四边形网格即可，网格密度可根据计算精度要求和计算能力的折衷来决定；

(2)将该尾迹分解为一系列不同频率、不同角度传播的平面前进波的线性叠加；其频率步长与角度范围由网格尺寸、水下目标的运动方向确定；用最小二乘法解得的系数分别对应于正弦、余弦级数表示的平面前进波的展开系数；

单一频率、单一传播方向的平面前进波的轨道速度矢量由线性波理论确定；

(3)用最小二乘法计算并累计所有网格节点处所有平面前进波的轨道速度矢量，即为水下运动目标尾迹的轨道速度矢量。

由两组正、余弦级数合成一组包含固定相位的余弦级数，其形式与表示风驱动态海面的余弦级数相同，便于水下目标尾迹和风驱动态海面的统一建模。

本发明首先将需要计算轨道速度的水下目标尾迹进行几何建模，然后采用最小二乘法结合线性波理论将尾迹分解为一系列不同频率、不同角度传播的平面前进波的线性叠加。在实施方式中按具体实例更具体地介绍了如何进行几何建模、最小二乘法平面前进波级数分解及轨道速度合成。

下面对各步骤的具体细节分别介绍如下：

(1)建立水下目标尾迹的动态几何模型

数值模拟的海面必须首先把海面离散化，即将海面剖分为四边形网格单元，每个网格节点的波高随时间变化，每个网格节点的流体微元轨道速度随时间变化。可根据计算精度要求和计算能力的折衷来决定剖分的网格密度。

在计算海面起伏高度时，假设模拟的海面面积为L_x×L_y，数值离散为N_x×N_y个网格节点，某网格节点(x,y)处在t时刻的波高z(x,y,t)由水下目标尾迹的波高z_w(x,y,t)和风驱海面起伏z_s(x,y,t)两部分构成：

z(x,y,t)＝z_w(x,y,t)+z_s(x,y,t) (1)

风驱海面的起伏通过叠加不同频率、不同方向的平面前进波得到：

其中

是平面前进波的波幅，k_n、ω_n、φ_m、ψ_mn分别为平面前进波的波数、角频率、传播方向和初始相位。S(ω_n,φ_m)是由实验观测数据统计得到的海浪方向谱。由于海浪谱能 量主要集中在较窄的波谱范围，本发明采用等量谱方法^[3]获取最佳的离散格式。波数k_n与角频率ω_n之间满足重力波的弥散关系：

其中g是重力加速度。

水下运动目标的静态尾迹，即某一确定时刻t₀的尾迹z_w(x,y,t₀)一般由点源速度势函数模型(例如Rankine卵形体模型)或计算流体力学软件得到。

假设该尾迹可分解为N_k个不同频率、M_φ个不同角度传播的平面前进波的线性叠加：

对于确定时刻t₀，和式通项中最后两项可以合并为一项相位项，故可简写为：

其中ψ'_mn＝ω_n't₀+ψ_mn，上式可进一步分解为不含相位项的形式，即一组正弦和一组余弦级数的形式，便于后续最小二乘法的计算：

式中

α₀＝A₀'

α_mn＝A_mn'cosψ'_mn

β_mn＝A_mn'sinψ'_mn (8)

求得系数α₀、α_mn、β_mn，即可将水下运动目标某一时刻的尾迹分解为余弦级数的 形式，或者说以余弦级数的形式重建水下运动目标在时刻t₀的尾迹。

(2)最小二乘法求解

水下目标尾迹余弦级数重建的关键是展开系数α₀、α_mn、β_mn的计算。本发明采用最小二乘法求解，即使得原尾迹z_w(x,y,t₀)与重建尾迹η_w(x,y,t₀)在最小二乘的意义下误差最小。除原理简单明了之外，该方法一个重要优点是海面网格节点可以非均匀分布，对于需要使用非均匀网格模拟复杂海面现象的模型同样适用。

本发明中最小二乘法的误差函数e定义如下：

为使误差函数最小，展开系数α₀、α_mn、β_mn必须满足如下方程组

代入(6)式可得如下最小二乘法线性方程组

(P^TP)X＝P^TY (11)

其中上标T代表矩阵或向量的转置运算，如果令N_xy＝N_x×N_y，N_kφ＝N_k×N_φ， 则P为N_xy×(2N_kφ+1)阶矩阵，可写为如下形式：

(12)

由于排版空间的限制，上式用C(x)代表cos(x)，用S(x)代表sin(x)，例如

求解最小二乘法线性方程组(13)，即得所需正、余弦级数的展开系数

(3)尾迹流体微元轨道速度的计算

水下运动目标的尾迹在给定时刻、给定位置的流体微元轨道速度也同样由不同频率、不同传播方向的平面前进波的轨道速度线性叠加确定。

在无限水深的假设条件下，对于单一频率、单位振幅的平面前进波，其速度势函数由线性波理论可表示为^[2]：

其对应的波高函数可表示为：

η_nm(x,y,t)＝cos[k_n(xcosφ_m+ysinφ_m)-ω_nt+ψ_nm] (14)

其轨道速度可表示为：

其中

(14)式中的η_nm与(6)式中的余弦函数完全相同，故将(14)式代入(6)式即可得到用最小二乘法重建尾迹的波高分布。同理，将(16)式代入(6)式即可得到用最小二乘法重建的海面流体微元的轨道速度矢量。

本发明中，具体模拟仿真时，用一个漂浮的球体标识流体微元的运动轨迹；用多个漂浮的球体标识流体微元的速度聚束现象；

漂浮的球体的球心与海面某一位置的流体微元重合，漂浮的球体的北极与该位置海 面的法线方向重合，球体的漂浮轨迹与流体微元的运动轨迹重合；

每一漂浮球体确定各自的漂浮轨迹，反映了海面相关区域流体的轨道速度变化规律。

本发明中，当海面风速大于3m/s时，需要模拟风场对海面波高起伏和流体微元轨道速度的影响；风驱海面的波高和海面流体微元轨道运动用线性波理论计算。

本发明方法能快速高效地数值仿真模拟风驱动态海洋环境中水下运动目标的尾迹及海面流体的轨道运动随时间的变化。除原理简单明了之外，该方法一个重要优点是海面网格节点可以非均匀分布，对于需要使用非均匀网格的复杂海面流体模型同样适用。

附图说明

图1是二维平面前进波的流体轨道运动示意图。

图2是三维平面前进波及其上多个漂浮小球的动态示意图。

图3是风驱动态海面上流体微元轨道速度矢量示意图。

图4是水下运动目标静态尾迹的最小二乘法重建示意图。其中，(a)潜艇静态尾迹，

(b)最小二乘法重建尾迹(N_k＝20，M_φ＝20)，(c)最小二乘法重建尾迹(N_k＝30，M_φ＝50)。

图5是风驱动态海面流体微元轨道速度重建示意图。其中，(a)风驱海面及其轨道运动矢量(线性波理论)，(b)风驱海面及其轨道运动矢量(最小二乘法重建1，N_k＝15，M_φ＝30)，(c)风驱海面及其轨道运动矢量(最小二乘法重建2，N_k＝25，M_φ＝40)。

图6是水下运动目标尾迹流体微元轨道速度矢量示意图。

图7是风驱动态海洋环境中水下运动目标动态尾迹示意图。

具体实施方式

先以二维单频平面前进波为例。图1为理想平面前进波在某一时刻的流体微元轨道运动示意图，平面前进波以相速度V_p向右传播，可见在波面不同位置蓝色流体微元轨道速度矢量V_orb(粉红色箭头)指向不同。如果用黄色漂浮小球代表某一特定流体微元，则其运动轨迹接近圆形(尽管平面前进波的波面在不断向右运动)。

图2是单一频率三维平面前进波的动态运动示意图，分别是t＝3s至t＝6s时刻的图像。平面前进波的周期T＝8s，相速V_p＝12.5m/s。水面有5个漂浮小球，标志所在点的流体微元轨道运动。当水面处于平静状态时小球沿x轴等间距15m分布；当水面存在平面 前进波时，可见小球随波浪动态起伏并且以各自的轨道速度运行，轨迹为图1所示圆形。小球之间的间距不断变化，反映了海面流体微元轨道运动所导致的速度聚束现象。

图3是风驱动态海面t＝3s至t＝6s时刻的图像。海面风速为10m/s，风向为侧风(与x轴夹角为-90°，即吹向观察者)，所用的海浪谱为经典PM(Pierson-Moskowitz)谱。海面每点的波高和轨道速度矢量是公式(2)所描述的不同频率、不同方向平面前进波的叠加，水面漂浮小球的分布同图2。可见风驱海面随时间动态起伏，海面流体微元的轨道运动对小球之间位置、间距等的影响非常复杂，同时可见流体微元轨道运动所导致的速度聚束现象(速度矢量在某些局部区域汇聚、指向同一方向，在某些局部区域发散、指向不同方向)。

图4是水下运动目标尾迹的最小二乘法重建示意图。原始的潜艇尾迹如图4a所示，由流体计算软件生成。所用模型是澳大利亚Collins级潜艇，艇长78.245m，舷宽7.8m，航速10节(5.144m/s)，航向负x轴，潜深6m。图4b和图4c是用最小二乘法重建的尾迹，可见随着余弦级数所用频率数目N_k和传播角数目M_φ的增加，重建尾迹与原尾迹的吻合度越来越高。

图5验证了最小二乘法重建风驱海面轨道运动矢量图的正确性。图中风驱动态海面处于t＝0.4s时刻，风速10m/s，风向与x轴夹角为180°。图5(a)的轨道运动矢量图是经典线性波理论的结果(同图3)，即由(2)式和(16)式得到。图5(b)和图5(c)用最小二乘法重建了海面波高和轨道运动矢量，即首先将图5(a)的海面用最小二乘法展开为(6)式所描述的平面前进波的叠加，再由展开系数、公式(16)、(2)和(6)重建海面波高和海面轨道速度矢量。图5(b)用的级数较少，分别为N_k＝15，M_φ＝30，可见无论海面波高、小球轨道运动所处的位置、还是轨道速度矢量均存在很大误差。图5(c)将级数的阶数分别提高到N_k＝25，M_φ＝40，可见与5(a)的吻合度已非常高。

图6是最小二乘法重建的潜艇尾迹的轨道运动矢量图，潜艇模型参数同图4。可见Bernoulli水丘、Kelvin臂等不同位置由流体微元轨道运动所导致的速度聚束现象。

图7是风驱海洋环境中水下运动目标尾迹t＝3s至t＝6s时刻的图像，潜艇模型参数同图4，海面风速10m/s，与x轴夹角为-90°。可见在风力和潜艇兴波的共同作用下，海面的动态起伏与海面流体微元的轨道运动展现了非常复杂的模式。

参考文献：

[1].Y-Q.Jin and F.Xu,Polarimetric scattering and SAR informationretrieval.Singapore:John Wiley&Sons,Mar.2013.

[2].J.N.Newman,Marine Hydrodynamics.Cambridge,Mass:M.I.T.Press,1977.

[3].F.L.W.Tang and C.F.Lin,“Practical method for evaluationdirectional spectra after shoaling and refraction,”in Proc.20th Conf.CoastalEng.,Taipei,pp.780-793,1986。

Claims (3)

1.一种水下目标动态尾迹与流体微元轨道速度的数值仿真方法，其特征在于具体步骤为：

(1)建立水下目标尾迹的动态几何模型，将模拟海面离散化，剖分为网格单元，每个网格节点的波高随时间变化，每个网格节点的流体微元轨道速度随时间变化；用点源速度势函数模型或计算流体力学软件建立水下运动目标的静态尾迹；

步骤(1)建立水下目标尾迹的动态几何模型的具体过程为：

假设模拟海面面积为L_x×L_y，数值离散为N_x×N_y个网格节点，某网格节点(x,y)处在t时刻的波高z(x,y,t)由水下目标尾迹的波高z_w(x,y,t)和风驱海面起伏z_s(x,y,t)两部分构成：

z(x,y,t)＝z_w(x,y,t)+z_s(x,y,t) (1)

风驱海面的起伏通过叠加不同频率、不同方向的平面前进波得到：

其中

是平面前进波的波幅，k_n、ω_n、φ_m、ψ_mn分别为平面前进波的波数、角频率、传播方向和初始相位；S(ω_n,φ_m)是由实验观测数据统计得到的海浪方向谱；波数k_n与角频率ω_n之间满足重力波的弥散关系：

其中g是重力加速度；

水下运动目标的静态尾迹，即某一确定时刻t₀的尾迹z_w(x,y,t₀)由点源速度势函数模型或计算流体力学软件得到；

(2)将静态尾迹分解为一系列不同频率、不同角度传播的平面前进波的线性叠加；其频率步长由网格尺寸决定，角度范围由水下目标的运动方向确定；最小二乘法解得的系数分别对应于正弦、余弦级数表示的平面前进波的展开系数；

单一频率、单一传播方向的平面前进波的轨道速度矢量由线性波理论确定；

步骤(2)的具体过程为：

假设该尾迹可分解为N_k个不同频率、M_φ个不同角度传播的平面前进波的线性叠加：

对于确定时刻t₀，和式通项中最后两项可以合并为一项相位项，故可简写为：

其中ψ_mn"＝ω_n't₀+ψ_mn'，上式可进一步分解为不含相位项的形式，即一组正弦和一组余弦级数的形式，便于后续最小二乘法的计算：

式中

求得系数α₀、α_mn、β_mn，即可将水下运动目标某一时刻的尾迹分解为余弦级数的形式，或者说以余弦级数的形式重建水下运动目标在时刻t₀的尾迹；

(3)用最小二乘法计算并累计所有网格节点处所有平面前进波的轨道速度矢量，即为水下运动目标尾迹的轨道速度矢量；

步骤(3)的具体过程为：

最小二乘法的误差函数e定义如下：

为使误差函数最小，展开系数α₀、α_mn、β_mn必须满足如下方程组：

代入(6)式得如下最小二乘法线性方程组：

(P^TP)X＝P^TY (11)

其中，上标T代表矩阵或向量的转置运算， 令N_xy＝N_x×N_y，N_kφ＝N_k×N_φ，则P为N_xy×(2N_kφ+1)阶矩阵，可写为如下形式：

上式中C(x)代表cos(x)，S(x)代表sin(x)，

求解最小二乘法线性方程组(11)，即得所需正、余弦级数的展开系数

尾迹流体微元轨道速度的计算：

水下运动目标的尾迹在给定时刻、给定位置的流体微元轨道速度同样由不同频率、不同传播方向的平面前进波的轨道速度线性叠加确定；

在无限水深的假设条件下，对于单一频率、单位振幅的平面前进波，其速度势函数由线性波理论表示为：

其对应的波高函数表示为：

η_nm(x,y,t)＝cos[k_n(xcosφ_m+ysinφ_m)-ω_nt+ψ_nm] (14)

其轨道速度表示为：

其中：

(14)式中的η_nm与(6)式中的余弦函数完全相同，将(14)式代入(6)式即得到用最小二乘法重建尾迹的波高分布；同理，将(16)式代入(6)式即得到用最小二乘法重建的海面流体微元的轨道速度矢量。

2.根据权利要求1所述的水下目标动态尾迹与流体微元轨道速度的数值仿真方法，其特征在于：

用一个漂浮的球体标识流体微元的运动轨迹；

用多个漂浮的球体标识流体微元的速度聚束现象；

漂浮的球体的球心与海面某一位置的流体微元重合，漂浮的球体的北极与该位置海面的法线方向重合，球体的漂浮轨迹与流体微元的运动轨迹重合；

每一漂浮球体确定各自的漂浮轨迹，反映了海面相关区域流体的轨道速度变化规律。

3.根据权利要求2所述的水下目标动态尾迹与流体微元轨道速度的数 值仿真方法，其特征在于：

当海面风速大于3m/s时，需要模拟风场对海面波高起伏和流体微元轨道速度的影响；风驱海面的波高和海面流体微元轨道运动用线性波理论计算。