理想声自导鱼雷运动要素的距离方位观测解算方法

李长文,初　磊

(海军潜艇学院,山东 青岛 266199)



理想声自导鱼雷运动要素的距离方位观测解算方法

李长文,初磊

(海军潜艇学院,山东 青岛 266199)

为快速预测鱼雷的声自导弹道,以特定时刻相对本舰的距离、舷角为运动要素,将其弹道参数化为以此要素为初值的微分方程组的解,论证了其解析解与一个代数方程等价,证明了此方程经变换后可用标准牛顿方法计算。考虑声速的作用,建立了发射、接收点分开的主动声纳探测模型,构造了基于距离、方位观测值解算鱼雷运动要素的相对较优的弹道拟合准则。实现了小样本条件下理想声自导弹道估计的快速计算。仿真实验证实,该方法预测的弹道可满足反鱼雷鱼雷射击参数计算的需要。

声自导鱼雷;运动要素解算;弹道拟合准则;反鱼雷鱼雷

水面舰船以反鱼雷鱼雷对抗潜射鱼雷之前,需要计算预先设计的反鱼雷鱼雷射击参数[1-3],其依据是水面舰船确认鱼雷报警后对鱼雷的观测。这些观测一般以主动声纳方式进行,可以给出鱼雷的距离与方位等信息。在一定的条件下,利用这些信息可以对鱼雷的弹道进行估计。基于鱼雷弹道过程,相关条件可以假设为:鱼雷处于直航、声自导、线导等运动状态,对鱼雷弹道的估计除当前时刻之前运动状态的估计外,更主要的是对当前时刻之后运动状态的预测[4-5]。

对应于直航情形,鱼雷运动状态可以用当前时刻位置、航速、航向这几个运动要素表达。在有距离、方位等观测信息的条件下,估计鱼雷弹道状态比较简单,这里所说的简单主要是模型简单、算法简单、计算速度快、精度高等。对应于声自导情形,这类问题将变得复杂,作者未发现相关报导。

声自导运动方式主要有被动、主动、主被动联合等,一般由变周期的分周期控制实现,控制的规则有固定提前角、变提前角等。若固定提前角为零,即直瞄式或尾追式跟踪,这是鱼雷被动声自导运动的主要制导方式,这时,考虑到制导周期很小,可以近似认为鱼雷速度方向总是鱼雷到本舰(噪声源)的方向,以下称之为理想声自导运动,本文只就这一情形进行相关问题的研究。研究主要包括理想声自导鱼雷运动状态的计算、无误差条件下主动声纳对鱼雷的距离方位观测值的计算、考虑误差条件下理想声自导鱼雷运动要素估计等。

1　理想声自导弹道状态及计算

假设水面舰船以速度v等速直航,鱼雷以速度vT对水面舰船进行理想声自导跟踪,速度方向为鱼雷到水面舰船噪声源的方向。以水面舰船舰艉(假设为噪声源)作为其位置参考点,对鱼雷进行观测的某个时刻为零时刻,假设零时刻鱼雷已经处于理想声自导运动状态。由于鱼雷声自导运动的空间范围较小,可以忽略地球椭球形状的影响,将水面舰船的运动看成零时刻所处水平面上的运动。鱼雷对水面舰船的跟踪分为水平面上的跟踪与竖直面上的跟踪,尾追阶段基本上已处于战斗深度,可以忽略其变深运动。假设只进行水平面上的跟踪,则运动规律与其在水面上的投影点的运动规律相同。因此,以零时刻水面舰船位置点为原点,在当地水平面上,以水面舰船速度方向为x轴正方向,其正右方向为y轴正方向,建立一个该水平面上的直角坐标系,以描述水面舰船及鱼雷的运动。为便于公式推导,平面直角坐标系的点(x,y)以及原点到这点的自由向量以复数x+iy表示,i为虚数单位。

如图1,在上述坐标系下,水面舰船总是位于实轴Ox轴上,时刻t到原点的距离为vt,所以表示水面舰船位置点Ov(t)的位置向量OOv(t)复数为实数vt,即

Ov(t)=vt

(1)

假设鱼雷运动状态以t时刻相对水面舰船的距离D(t)、舷角X(t)表示,则其绝对位置点G(t)的复数为

G(t)=vt+D(t)eiX(t)

(2)

式中:i为虚数单位。

图1　理想声自导运动

1.1鱼雷运动状态的微分方程组

鱼雷运动的速度向量复数为

G′(t)=v+[D′(t)+D(t)iX′(t)]eiX(t)

(3)

式中:G′(t),D′(t),X′(t)表示这些时间变量函数的导数。

由前面关于理想声自导的假设,鱼雷速度方向相对Ox轴的角为X(t)+π,所以其速度向量复数的另一表示为vTei(X(t)+π),将这一表示与G′(t)相等,则

-vT=ve-iX(t)+D′(t)+D(t)iX′(t)

(4)

改成实部与虚部表示的分量形式:

(5)

此即鱼雷相对水面舰船运动状态所满足的微分方程组。给定任一时刻t0的距离D(t0),舷角X(t0),总有唯一解,(D(t),X(t)),t≥0,因此(D(t0),X(t0))可以作为理想声自导鱼雷的运动要素。运动状态的计算就是解这个微分方程组,在解算鱼雷运动要素的搜索算法中要不断调用解这个方程组的计算程序。若用数值方法计算,一方面计算误差随着时间变长而增加,更主要的是计算速度远不能满足需要。因此,必须寻找解微分方程组(5)的精确快速算法。

考虑到对称性,只研究0≤X(t)≤π的情形。若X(t0)=0(或π),即自水面舰船正前方(或正后)方跟踪目标,鱼雷将作直线运动,且舷角保持X(t)=X(t0),D(t)=D(t0)-(v+vT)(t-t0)(正前方),或D(t)=D(t0)-(v-vT)(t-t0)(正后方),自t0时刻至鱼雷追上水面舰船时刻的跟踪时间tg=D0/(v+vT)(正前方),或tg=D0/(v-vT)(正后方)。下面假设t0>0,D(t0)>0,0

(1/D)dD/dX=-kcscX-cotX

(6)

式中:k=vT/v,此方程有解析解:

D=D(t0)sinX(t0)tank[X(t0)/2]cotk(X/2)cscX

(7)

只要能计算X(t),就可以用此式计算D(t)。因为式(6)正比于cosk-1(X/2)/sink+1(X/2),所以,D→0⟺k>1,X→π,即当且仅当鱼雷速度大于水面舰船速度,理想声自导运动的鱼雷才能追上目标,且总是自正后方追上目标,以下假设k>1。

将式(7)代入式(5)可得X(t)的微分方程:

csc2X(t)cotk[X(t)/2]X′(t)=r

(8)

式中:

r=[v/D(t0)]cscX(t0)cotk[X(t0)/2]

(9)

式(8)两边自t0到t积分,得:

(10)

此为X(t)所满足的方程,计算鱼雷弹道状态的关键就是解这个方程。式(10)中的时间满足-t0≤t-t0≤tg,tg为跟踪时间,只需在式(10)令t→t0+tg,X(t)→π,即可得:

(11)

式(11)即为计算理想声自导跟踪时间的公式,对X(t0)=0(或π)的情形同样成立。

1.2鱼雷弹道状态的计算

前文已将微分方程组(5)的求解变为一个关于未知数X(t)的代数方程(10)的求解,虽然这个方程无解析解,但可以用比较标准的方法精确快速地计算。

利用r,tg的表达式可将式(10)改为

(12)

令z=cot[X(t)/2]∈(0,∞),则

X(t)=2arccotz

方程(12)变为

(13)

因为g′(z)=zk-2+zk>0,所以g(z)为(0,∞)上单调上升函数。因为g(0)<0,g(∞)=∞,所以g(z)=0在(0,∞)上有唯一解。

用Newton方法计算一元方程的数值解十分快速和精确,对于单调上升(下降)的函数,只要其二阶导数大于(小于)0,起点的选择可不作要求。对于函数g(z),考虑到g″(z)=(k-2)zk-3+kzk-1,若k≥2,则g″(z)>0,无论方程的解小于z0(t0

定理1设b>a>0,函数y=g(z)在[a,b]上有二阶连续导数,且g′(z)>0,z ∈[a,b],则存在α>0,使f(u)=g(uα),满足f′(u)>0,f″(u)>0,u∈[a1/α,b1/α]的充要条件为p=minz∈[a,b]zg″(z)/g′(z)>-1,且可取α=1/(1+q),q∈(-1,p]。

证明∀α>0,u∈[a1/α,b1/α]

f′(u)=g′(uα)αuα-1>0

f″(u)=g″(uα)(αuα-1)2+g′(uα)α(α-1)uα-2⟹

f′(u)>0,f″(u)≥0,∀u∈[a1/α,b1/α]⟺

f″(u)≥0,∀u∈[a1/α,b1/α]⟺

uαg″(uα)/g′(uα)+1-1/α≥0,u∈[a1/α,b1/α]⟺

0<1/α≤1+zg″(z)/g′(z),z∈[a,b]⟺

p=minz∈[a,b]zg″(z)/g′(z)>-1

且

α≥1/(1+q),q ∈(-1,p]

证毕。

式(13)的g(z)满足:

g′(z)>0,0

取α=1/(1+k-2)=1/(k-1),则

满足:

f′(u)>0,f″(u)≥0,u∈[0,u0],u0=z0k-1

可以用u0为起点的Newton方法计算方程f(u)=0的数值解u*,迭代格式为

ul+1=ul-f(ul)/ f′(ul),l=0,1,…

(14)

停止条件为l>lmax或|f(ul)|<ε,实验程序内取lmax=5,ε=10-15,实际计算时取uα=eαlnu,还能提高计算速度。图2是k<2的一次实验对应的g(z),f(u)的图形,图中,u,f(u),z,g(z)均为无量纲量。

g(z)=0的数值解为z*=(u*)α,以此可以计算t≤t0+tg时刻的X(t),D(t),之后以式(2)计算鱼雷绝对位置点复数,以公式:

vG(t)=-vTeiX(t)

(15)

计算鱼雷速度向量复数。

图2　k=1.2113时g(z),f(u)曲线

2　主动声纳对理想声自导鱼雷的观测

设水面舰船以主动声纳探测鱼雷,信号发射点P(t)相对舰艉的有向线段数量为Ls,信号接收点Q(t)相对舰艉的有向线段的数量为Lr,则其对应的复数为

P(t)=vt+Ls, Q(t)=vt+Lr

(16)

Ls或Lr的符号为正表示相应的点位于舰艉前方,符号为负表示相应的点位于舰艉后方,为简单起见,假设这些点与鱼雷位于同一水平面上。记主动信号发射时刻为ts,j,j=1,2,…,ts,j时刻由P(ts,j)点发射的信号到达鱼雷的时刻记为tr,j,假设到达信号在G(tr,j)点被鱼雷立即反射,反射信号到达接收点的时刻记为ta,j。假设声波在水中的传播可近似为直线,平均声速为常数,记为vs,则

|G(tr,j)-P(ts,j)|=vs(tr,j-ts,j)

(17)

|Q(ta,j)-G(tr,j)|=vs(ta,j-tr,j)

(18)

给定j,z=tr,j-ts,j,由式(17)有:

(19)

其解z*可利用Newton方法计算,格式同式(14),起点可取1.5|G(ts,j)-P(ts,j)|/vs,方程(19)左端函数的导数计算公式为

ta,j为方程(18)的解,可以用上述方法计算。但由于接收点Q(t)等速直航,与方程(17)不同,故方程(18)为一个关于δ=ta,j-tr,j的二次方程,因此有简单的解析计算公式。

定理2以A(t),B(t)表示2个水下运动物体t时刻位置点的空间坐标,记w=-1(A发B收),1(B发A收),设B以速度vB=(vBxvByvBz)匀速运动,‖vB‖

式中:U=A(tA)-B(tA),q=U(vB/ vs)

证明由已知:

‖A(tA)-B(tB)‖=vs(tA-tB)w

因为B(t)=B(0)+vBt,所以

B(tB)=B(tA)-vB(tA-tB)

‖U‖2+2(UvB)(tA-tB)+‖vB‖2(tA-tB)2

用求根公式计算tA-tB,因为‖vB/vs‖<1,根式前面的正负号与w相同。证毕。

利用上述定理可得:

(20)

式中:Uj=G(tr,j)-Q(tr,j),q=ksReUj,ks=v/vs。

接收点于ta,j时刻收到信号的理论距离及到达方位(此处以舷角表示)为

d0j=|G(tr,j)-Q(ta,j)|

(21)

β0j=arg[G(tr,j)-Q(ta,j)]

(22)

假设无系统误差,只考虑随机观测误差,这些理论值的观测值为

dj=d0j+Δdj

(23)

βj=β0j+Δβj

(24)

为简单起见,假设这些观测值都在可观测范围内,Δdj,Δβj,j=1,2,…,相互独立,且

3　理想声自导运动要素的距离方位观测估计

所谓距离方位观测解算,就是用当前时刻t0之前(ta,j

即可以只用观测值dj,j=1,…,n,求出(D(t0),X(t0))的极大似然估计,且这个估计与以Jd为目标函数的最小二乘法估计相同。

同样也可以只用βj,j=1,…,n,给出(D(t0),X(t0))的极大似然估计,且与以

为目标函数的最小二乘法估计相同。

为便于反鱼雷鱼雷发现鱼雷以及发现鱼雷时刻的态势有利于后续的跟踪,只要有可能,就应利用所有观测信息对鱼雷运动状态给出精确的预测[6]。因此必须综合考虑样本(dj,βj),j=1,…,n,其似然函数为

对于σd,σβ已知的情形,这对应于以

为目标函数的加权最小二乘法估计,对于σd,σβ未知的情形,Jd,β不能作为目标函数,只能以L为目标函数求(D(t0),X(t0))的极大似然估计,这其中包含对σd,σβ的估计。

数值实验发现,搜索以上这几种目标函数的最优点需要较长的时间,且对于样本容量n 较小的情形,估计的精度不足。为此,考虑一种弹道拟合的方法,记

zj=djexp(iβj), j=1,…,n

(25)

为ta,j时刻鱼雷相对于接收点的观测位置复数,则复随机变量zj的数学期望E(zj)为

E(zj)=E[(d0j+Δdj)exp(i(β0j+Δβj))]=

E(d0j+Δdj)E[exp(i(β0j+Δβj))]=

d0jexp(iβ0j)E[exp(iΔβj)]

由正态分布特征函数在点1处的值可知:

(26)

所以

(27)

因为

|zj-E(zj)|2=|(d0j+Δdj)exp[i(β0j+Δβj)]-

所以zj的方差Var(zj)为

(28)

式(27)、式(28)说明距离方位观测条件下,观测位置点并非理论位置点的无偏估计,且是异方差的,为平衡各次观测,可令

(29)

为目标函数,以

minJ(D(t0),X(t0))

s.t. (D(t0),X(t0))∈(0,dmax)×(0,π)

的解D*(t0),X*(t0)作为D(t0),X(t0)的估计,此处称之为弹道拟合准则。数值实验发现,这个最优问题的求解速度较其他方法快,且以最近的15～25个观测样本,能给出满足射击参数计算等要求的弹道估计。图3为一次计算的图形(输入参数,略)。

图3　一次仿真实验预测的鱼雷弹道

4　结束语

本文构造了对理想声自导运动的鱼雷运动要素进行基于距离方位观测的估计方法,用仿真实验进行了验证。该方法对于解算等速直航运动或其他参数化弹道的鱼雷运动要素同样适用。但这类方法的

实质为非线性最优化问题,其算法的速度一般不能满足实时计算的需要,对于给定的算法速度,必须设置发射反鱼雷鱼雷的时刻较当前时刻t0有一点延迟,给计算预留一定的时间。正是这一原因,需要尝试一些提高算法速度的办法。

[1]李本昌,李长文.预定相遇态势的潜射鱼雷射击及其技术方法[J].指挥控制与仿真,2013,35(2):5-8.

LI Ben-chang LI Chang-wen.Firing and technical method based on given meeting situation of sub-launched torpedo[J].Command Control & Simulation,2013 35(2):5-8.(in Chinese)

[2]李长文,李本昌,初磊.预定相遇态势的鱼雷二次转角射击通用模型与算法[J].鱼雷技术,2013,21(4):306-312.

LI Chang-wen,LI Ben-chang,CHU Lei.General model and algorithm of two-time turn angle shooting of torpedo for predetermined encountering situation[J].Torpedo Technology,2013,21(4):306-312.(in Chinese)

[3]由大德,徐德民.反鱼雷鱼雷拦截概率影响因素仿真分析[J].鱼雷技术,2010,18(4):312-315.

YOU Da-de,XU De-min.Analysis on parameters influencing intercept probability of anti-torpedo torpedo[J].Torpedo Technology[J],2010,18(4):312-315.(in Chinese)

[4]陈颜辉,朱伟良,杜毅.潜射鱼雷弹道预测模型与仿真[J].海军工程大学学报,2013,25(2):57-61.

CHEN Yan-hui,ZHU Wei-liang,DU Yi.Trajectory prediction model and simulation of sub-launched torpedo[J].Journal of Naval University of Engineering,2013,25(2):57-61.(in Chinese)

[5]姜凯峰,周明,林宗祥.舰载反鱼雷鱼雷作战效率仿真研究[J].计算机仿真,2011,28(6):9-12.

JIANG Kai-feng,ZHOU Ming,LIN Zong-xiang.Operational efficiency based on simulation of shipboard anti-torpedo torpedo[J].Computer Simulation,2011,28(6):9-12.(in Chinese)

[6]李宗吉,张西勇,练永庆.反鱼雷鱼雷自适应滑模导引律[J].弹道学报,2015,27(1):24-26.

LI Zong-ji,ZHANG Xi-yong,LIAN Yong-qing.Antit-orpedo-torpedo adaptive sliding mode guidance[J].Journal of Ballistics,2015,27(1):24-26.(in Chinese)

Method of Solving Motion Factors of Ideal Acoustic Homing Torpedo by Distance and Bearing Observations

LI Chang-wen,CHU Lei

(Navy Submarine Academy,Qingdao 266199,China)

To predict the trajectory of acoustic homing torpedo fastly,the distance and bearing relative to this ship for specified time were taken as the torpedo motion factors,and the trajectory was parameterized to the solution of differential equations taking the factors as initial value.The result demonstrates that the analytic solution is equivalent to an algebraic equation,and this transformed algebraic-equation can be calculated using standard Newton method.A detection model of active sonar with the separation of sending and receiving point was established by considering sound speed influence.Based on the observation value of distance and azimuth,a relatively superior criterion named trajectory fitting for solving the motion factor was constructed.The fast calculation was realized for estimating the ideal acoustic homing trajectory under the conditions of small samples.Simulation results show that the torpedo trajectory predicted by this method can meet the parameters calculation for firing anti-torpedo torpedo.

acoustic homing torpedo;solution to motion factor;trajectory fitting criterion;anti-torpedo torpedo

2016-01-25

李长文(1962- ),男,副教授,研究方向为随机数学及军事运筹学。E-mail:lcwqtxy2996@sina.com。

O315

A

1004-499X(2016)03-0043-05