基于射流微分方程組的水炮命中問題求解*

張興江，王 軍

（南京理工大學(xué)先進(jìn)發(fā)射協(xié)同創(chuàng)新中心，南京 210094）

0 引言

隨著海上領(lǐng)土糾紛的升級和維權(quán)活動的增多，消防水炮作為一種非殺傷性武器得到了越來越多的關(guān)注和應(yīng)用。消防水炮可以用來驅(qū)逐外來騷擾和自衛(wèi)反擊，同時又可以避免傷亡和引起戰(zhàn)爭，是海上出巡和海上捕魚作業(yè)時理想的武器裝備。當(dāng)前，我國海洋維權(quán)局勢日益嚴(yán)峻，維權(quán)任務(wù)十分艱巨，中國海監(jiān)海上維權(quán)執(zhí)法力量亟須加強(qiáng)。對此有分析人士認(rèn)為，在中國現(xiàn)有海監(jiān)船無武器配備情況下，高壓水炮是巡航執(zhí)法最有利的裝備，可有效震懾違法船只[1]。

智能消防水炮在滅火系統(tǒng)中有著大量的應(yīng)用與研究，但其較多地集中在射流方程的建立。目前命中問題求解主要有基于數(shù)學(xué)模型和基于計算機(jī)視覺兩種方法。基于計算機(jī)視覺的方法[1]，主要是通過圖像識別對射流落點和目標(biāo)進(jìn)行閉環(huán)控制，但其在調(diào)整速度上仍不及基于數(shù)學(xué)的方法，在此不進(jìn)行分析。水炮中基于數(shù)學(xué)模型的計算方法的研究大多為針對靜止目標(biāo)、固定高度，文獻(xiàn)[2]雖然針對不同高度的目標(biāo)進(jìn)行了研究，但其根本仍然為事先固定好高度，不能跟隨目標(biāo)高度變化實時修正。

對于運(yùn)動目標(biāo)的命中問題求解，在火炮控制中有著大量的研究[3]，但因為求解時間限制，火炮解算中使用更多的是射表法。而對于水炮來說，其解算周期長，精度明顯低于火炮，所以可以使用求解射流微分方程組的方法來求解水炮的命中問題。

1 船用水炮射擊解命中問題描述

解命中問題的實質(zhì)就是確定射流與運(yùn)動目標(biāo)在空間相遇點即提前點的坐標(biāo)。將目標(biāo)航跡的數(shù)學(xué)模型與射流的數(shù)學(xué)模型聯(lián)立求解，可以求得命中點坐標(biāo)，進(jìn)而換算出水炮的射擊諸元。

圖1 解命中問題描述

假設(shè)該時刻對應(yīng)的命中點為Mq。已知經(jīng)過濾波得到射擊瞬時目標(biāo)在大地坐標(biāo)系下的位置坐標(biāo)和速度，目標(biāo)運(yùn)動TF時間后到達(dá)命中點Mq，Mq點在水炮大地坐標(biāo)系下的直角坐標(biāo)為，計算得到Mq點在水炮大地坐標(biāo)系下的極坐標(biāo)為，且對應(yīng)的水平距離和高度分別為dq1和hq1。目標(biāo)在水炮大地坐標(biāo)系下方位角β，高低角α。

1.1 目標(biāo)航跡數(shù)學(xué)模型

考慮到船和目標(biāo)都是慢速目標(biāo)，在一個解算周期內(nèi)往往沒有明顯的機(jī)動，所以解算時將目標(biāo)與船的運(yùn)動看作勻速運(yùn)動。因此，每個解算周期內(nèi)目標(biāo)航跡的數(shù)學(xué)模型便是以目標(biāo)相對于水炮的坐標(biāo)為起點，以目標(biāo)相對于水炮的速度做勻速直線運(yùn)動，即目標(biāo)在水炮大地坐標(biāo)系內(nèi)做勻速直線運(yùn)動。航跡的數(shù)學(xué)模型如下：

其中，Vs、Xs分別為目標(biāo)在水炮大地坐標(biāo)系的速度、位置矢量，d為目標(biāo)距炮口的水平距離，hs為目標(biāo)距炮口的垂直距離。

1.2 射流數(shù)學(xué)模型

目前大多數(shù)水射流運(yùn)動模型都基于牛頓第二定律建立。無論是考慮水射流中的一段射流微團(tuán)或是考慮射流中的水滴，一般都只考慮空氣阻力和重力[2]。針對各種情況下的射流方程也有著許多的研究，但對于水炮不同型號、不同流量目前并沒有一種通用的數(shù)學(xué)模型，常用的方法是利用射流運(yùn)動學(xué)公式進(jìn)行推導(dǎo)，通過射流實測數(shù)據(jù)確定射流方程中的待定系數(shù)，從而獲得適用于此種水炮的射流模型。這種方法雖然通用性差，但針對特定的水炮也可以獲得較高的精度。

顯然射流的運(yùn)動可以分解為垂直面和水平面兩種運(yùn)動，二者互不影響，求解射流方程時，為便于分析，將水炮大地坐標(biāo)系繞Zp軸旋轉(zhuǎn)得到射流坐標(biāo)系Xh-Yh-Zh，使射流沿Yh-Zh平面射出，則射流垂直面內(nèi)的運(yùn)動將完全在Yh-Zh平面內(nèi)，水平面內(nèi)的運(yùn)動將完全在Xh-Yh平面內(nèi)，即Yh軸反映射流的水平距離，Zh軸反映射流的垂直距離，Xh軸反映射流的橫向偏移量，如圖2所示。

圖2 射流在射流坐標(biāo)系的運(yùn)動

射流運(yùn)動模型的不同關(guān)鍵在于空氣阻力的選取不同，不同的空氣阻力公式可以構(gòu)造不同的射流方程，其能達(dá)到的精度也各有不同。HATTON等[4-5]提出了一種空氣阻力模型，證明其可以模擬大流量消防炮的射流軌跡，而船用水炮正是一種大流量消防炮，故本文用此種模型建立射流方程，其空氣阻力模型為：

式中，S為射流軌跡長，滿足

其中，F(xiàn)yz為Yh-Zh平面內(nèi)單位質(zhì)量的空氣阻力，F(xiàn)x為Xh-Yh平面內(nèi)單位質(zhì)量的空氣阻力，待定系數(shù)k、b、A、B可以參考文獻(xiàn)[4-6]由實驗數(shù)據(jù)獲得。

按運(yùn)動學(xué)規(guī)律，射流的數(shù)學(xué)模型最終建立如下：

2 命中方程求解

命中方程的求解即聯(lián)立射流微分方程組式（5）和目標(biāo)航跡數(shù)學(xué)模型式（1）進(jìn)行數(shù)值求解，可以同時獲得射流與目標(biāo)相遇時間tf和水炮射擊諸元。其中求解微分方程組的方法很多，如歐拉矩形法、梯形法、龍塔庫塔法等，在此，選用常用的四階龍塔庫塔法進(jìn)行數(shù)值求解[7-8]。

2.1 初始條件的求取

要求解射流微分方程組，初始條件至關(guān)重要。在水炮的射流命中問題中，已知條件有：

這是典型的兩點邊值問題，如何從邊值條件（6）出發(fā)，求解出水炮射擊諸元和射流飛行時間tf，正是水炮命中問題求解的任務(wù)。在此，以射流落點為依據(jù)，使用變步長的四階龍塔庫塔法，將邊值問題變?yōu)槌踔祮栴}迭代求解。

2.1.1 初始高低角αt和軌跡長S0

軌跡長S0=0；高低角αt初值選取為目標(biāo)坐標(biāo)對應(yīng)的高低角，即

船用水炮在射擊時往往采用固定的流量和壓力，因此，射流相對炮口初速vp0便也確定[9]，后面也不考慮vp0變化的情況。射流初速主要影響因素有：vp0、風(fēng)速、船的運(yùn)動、水炮的轉(zhuǎn)動。對于船的姿態(tài)角在解算周期內(nèi)變化很小，造成的影響可以忽略。

圖3 合成速度示意圖

其中，vw為風(fēng)速大小，βw為風(fēng)向，vb為船速大小，βw為船航向，βq為目標(biāo)點相對于水炮的方位角。由于風(fēng)的垂直分量比較小，因此，風(fēng)只在水平面上分解為橫風(fēng)和縱風(fēng)，同時船基本沒有垂直方向的運(yùn)動。

假設(shè)水炮在每個解算周期內(nèi)的運(yùn)動為勻速轉(zhuǎn)動，忽略加減速過程，那么由于水炮轉(zhuǎn)動產(chǎn)生的射流初速度便不能忽略?？紤]水炮在持續(xù)射擊的過程中，每個運(yùn)動周期水炮高低轉(zhuǎn)動的角度很小，所以水炮轉(zhuǎn)動對vy、vz產(chǎn)生的影響可以忽略。

解算周期開始獲得的水炮方位角β0，對方位角作線性預(yù)測，假設(shè)解算周期結(jié)束對應(yīng)架位β，可得水炮的方位速度前饋為

所以射流在炮口處的線速度：

其中，T為解算周期，Vβmax為水炮方位轉(zhuǎn)動最大角速度，L為水炮身管長。

所以射流初速為

2.2 求解流程

2）使用四階龍塔庫塔法求解射流微分方程組式（5），設(shè)定初始步長為0.01 s。每次求解過程中有關(guān)射流諸元序列值：

i為以步長h的計算序列號。求解終止條件為：

當(dāng)求解終止時，此時i對應(yīng)的序列值則為射流在射流坐標(biāo)系的落點坐標(biāo)和速度，飛行時間：

3）利用求得的tf和目標(biāo)航跡模型式（1）獲得目標(biāo)未來點坐標(biāo)，誤差分別為：

4）對高低角修正，新的射流高低角為

αt為上一次龍塔庫塔法求解的高低角初值。

利用新的目標(biāo)未來點和修正量確定初值。

3 實驗仿真分析

針對上述數(shù)學(xué)模型和命中求解流程，使用matlab2014進(jìn)行實驗仿真，驗證求解方法的可行性，并進(jìn)行相關(guān)的分析。

仿真時相關(guān)參數(shù)如下：

各空氣阻力系數(shù)：k=0.000 8，b=0.023，A=20，B=0.04，g=9.8（在此省略單位），水炮部分參數(shù)：L=0.5 m，vp0=50 m/s，解算周期 T=0.2 s，解算精度0.1 m。

因為解算過程縱風(fēng)和船速最終會分解到對應(yīng)軸，本質(zhì)相當(dāng)于船靜止。所以設(shè)定固定橫風(fēng)vwx=2m/s，船速vb=0，縱風(fēng)vyz=0。

水炮高低角為30°時，射流軌跡如圖4所示。

圖4 射流軌跡

3.1 靜止目標(biāo)、勻速運(yùn)動和勻加速運(yùn)動目標(biāo)的對比

各種目標(biāo)仿真參數(shù)如下：

靜止目標(biāo)航跡：位置（10，50，-3），勻速運(yùn)動目標(biāo)航跡：起始點（10，50，-3），速度（3，2，0.1），勻加速運(yùn)動目標(biāo)航跡：起始點（10，50，-3），速度（3，2，0.1），加速度（1，-1，0）。

3.1.1 單點的解算誤差

對于解算點來說，解算誤差和射擊諸元仿真如圖5～圖7所示。

圖5 靜止目標(biāo)解算點誤差和射擊諸元

圖6 勻速運(yùn)動目標(biāo)解算點誤差和射擊諸元

圖7 勻加速目標(biāo)解算點誤差和射擊諸元

從以上仿真結(jié)果可以看出，因設(shè)定迭代誤差為0.1 m，假定目標(biāo)求解過程做勻速運(yùn)動，所以當(dāng)目標(biāo)實際靜止或做勻速運(yùn)動時，解算誤差小于0.1 m；而當(dāng)目標(biāo)做勻加速運(yùn)動時，可以看到解算誤差大于0.1 m，且具有隨時間誤差增大的趨勢，這是由于目標(biāo)航跡模型假定不一致造成的，增大的誤差主要是加速度產(chǎn)生的位移，符合公式，增大趨勢是因為隨著目標(biāo)越來越遠(yuǎn)，射流飛行時間越來越長，加速度產(chǎn)生的影響也越來越大。

3.1.2 整個解算周期內(nèi)勻速轉(zhuǎn)動的解算誤差

對于解算周期內(nèi)勻速轉(zhuǎn)動水炮來說，以目標(biāo)做勻速運(yùn)動為例，解算誤差如圖8所示。

圖8 勻速運(yùn)動目標(biāo)解算周期誤差和射擊諸元

與解算點誤差進(jìn)行對比，解算周期內(nèi)實際射流落點和目標(biāo)誤差變化趨勢與解算點誤差變化趨勢相同，只是誤差變大。原因在于水炮在每個解算周期內(nèi)勻速轉(zhuǎn)動，因水炮高低角不變，所以射流的水平距離不變，隨著目標(biāo)的運(yùn)動，誤差便會越來越大，到達(dá)下一個周期，誤差便會重新變小，然后逐漸變大，所以呈現(xiàn)的結(jié)果就是鋸齒狀誤差。而對于靜止目標(biāo)來說，因沒有轉(zhuǎn)動，所以誤差不變。

3.2 不同解算周期仿真對比

以目標(biāo)勻速運(yùn)動為例，不同的解算周期對于解算點的誤差影響不大，都保持在0.1 m以內(nèi)，而對于每個解算周期內(nèi)勻速轉(zhuǎn)動的情況有明顯影響，仿真結(jié)果如圖9所示。

圖9 勻速運(yùn)動目標(biāo)不同解算周期誤差

其中，按從左往右、從上往下順序，解算周期分別為：0.1 s、0.2 s、0.5 s、1 s。

從圖中可以看出，隨著解算周期的增大，解算誤差越來越大，原因在于解算周期越大，周期內(nèi)水炮勻速轉(zhuǎn)動的時間越長，積累的誤差便會越大。同理，可以知道對于固定解算周期，目標(biāo)速度越快，誤差越大。

3.3 不同解算精度的對比

求解過程中迭代誤差設(shè)定直接影響解算精度，當(dāng)設(shè)定0.1 m時，解算點精度在0.1 m內(nèi)，當(dāng)設(shè)定0.5 m時，解算點解算精度在0.5 m內(nèi)，而在整個解算周期內(nèi)，解算誤差仿真結(jié)果如圖10所示。

圖10 解算精度0.5時解算誤差

與圖8對比，迭代誤差0.1 m和0.5 m時，解算誤差變化不大。所以可以根據(jù)解算誤差適當(dāng)?shù)倪x擇迭代誤差，以減少迭代次數(shù)和時間。

3.4 解算時間

按目標(biāo)靜止，設(shè)定距離，解算時間如表1所示。

表1 不同距離解算時間

從表1可以看出，目標(biāo)距離越遠(yuǎn)，解算時間越長。原因在于步長和高低角的調(diào)整不理想。但總體來說，選擇適當(dāng)?shù)慕馑阒芷诤途嚯x，解算時間是可以滿足要求的。同時降低解算精度可以明顯減少解算時間。

4 船姿態(tài)角變化的影響

船在海上行駛時，勢必會因為海面的運(yùn)動造成船姿態(tài)角的變化，相應(yīng)的水炮姿態(tài)角也會發(fā)生變化。由前面所述，通過命中方程求解可以求得水炮在水炮大地坐標(biāo)系下的射擊諸元，但因為水炮的姿態(tài)運(yùn)動，水炮最終輸出的射擊諸元不能不發(fā)生變化。如果水炮驅(qū)動控制系統(tǒng)自帶有處理水炮姿態(tài)變化的程序，那么只需輸出水炮在大地坐標(biāo)系下的射擊諸元即可；但如果其沒有，在命中問題求解后輸出的射擊諸元必須進(jìn)行相應(yīng)的處理，以抵消水炮姿態(tài)角變化帶來的影響。

通過船體測量系統(tǒng)和相應(yīng)的姿態(tài)濾波等處理程序，可以獲得水炮射擊時的水炮姿態(tài)角，然后只需將水炮的射擊諸元由水炮大地坐標(biāo)系進(jìn)行姿態(tài)坐標(biāo)變換，轉(zhuǎn)換到水炮的載體坐標(biāo)系即可[3]，此時水炮在水炮載體坐標(biāo)系下的射擊諸元即為輸出到水炮驅(qū)動控制系統(tǒng)的射擊諸元。

5 結(jié)論

提出了基于射流微分方程組水炮的實時命中問題求解方法，介紹了目標(biāo)航跡數(shù)學(xué)模型和射流數(shù)學(xué)模型的建立，對于影響射流的因素：風(fēng)、船的運(yùn)動、水炮轉(zhuǎn)動等統(tǒng)一處理，使用四階龍塔庫塔法求解微分方程組，并通過修改步長和高低角不斷迭代，直至滿足精度要求，獲得水炮射擊諸元。通過matlab仿真驗證，其解算誤差在0.1 m內(nèi)，證明了方法的可行性，為水炮自動控制系統(tǒng)的開發(fā)提供了一種可行思路。

通過實驗對不同目標(biāo)航跡、不同解算周期和解算精度等情況的仿真，得出了不同情況下的解算誤差，為水炮的實際應(yīng)用，如數(shù)學(xué)模型的建立、控制策略的選擇、解算精度和解算時間的合理優(yōu)化等，提供了一定的參考價值。

但在命中求解過程中仍有許多可以改進(jìn)的地方[10-14]，比如龍塔庫塔法步長的選擇、高低角的修正、初值的確定等，并且射流微分方程可以進(jìn)行簡化，上述的微分方程組具有7個變量，可以將其減少或分離，以減少解算的時間。同時在每個解算周期內(nèi)水炮勻速轉(zhuǎn)動的控制策略會造成一定誤差，可以嘗試其他控制策略以優(yōu)化誤差。