張興江，王 軍

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

0 引言

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

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

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

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

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

圖1 解命中問題描述

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

1.1 目標航跡數(shù)學模型

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

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

1.2 射流數(shù)學模型

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

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

圖2 射流在射流坐標系的運動

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

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

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

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

2 命中方程求解

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

2.1 初始條件的求取

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

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

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

軌跡長S0=0；高低角αt初值選取為目標坐標對應的高低角，即

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

圖3 合成速度示意圖

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

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

解算周期開始獲得的水炮方位角β0，對方位角作線性預測，假設解算周期結束對應架位β，可得水炮的方位速度前饋為

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

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

所以射流初速為

2.2 求解流程

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

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

當求解終止時，此時i對應的序列值則為射流在射流坐標系的落點坐標和速度，飛行時間：

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

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

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

利用新的目標未來點和修正量確定初值。

3 實驗仿真分析

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

仿真時相關參數(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。

因為解算過程縱風和船速最終會分解到對應軸，本質相當于船靜止。所以設定固定橫風vwx=2m/s，船速vb=0，縱風vyz=0。

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

圖4 射流軌跡

3.1 靜止目標、勻速運動和勻加速運動目標的對比

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

靜止目標航跡：位置（10，50，-3），勻速運動目標航跡：起始點（10，50，-3），速度（3，2，0.1），勻加速運動目標航跡：起始點（10，50，-3），速度（3，2，0.1），加速度（1，-1，0）。

3.1.1 單點的解算誤差

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

圖5 靜止目標解算點誤差和射擊諸元

圖6 勻速運動目標解算點誤差和射擊諸元

圖7 勻加速目標解算點誤差和射擊諸元

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

3.1.2 整個解算周期內勻速轉動的解算誤差

對于解算周期內勻速轉動水炮來說，以目標做勻速運動為例，解算誤差如圖8所示。

圖8 勻速運動目標解算周期誤差和射擊諸元

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

3.2 不同解算周期仿真對比

以目標勻速運動為例，不同的解算周期對于解算點的誤差影響不大，都保持在0.1 m以內，而對于每個解算周期內勻速轉動的情況有明顯影響，仿真結果如圖9所示。

圖9 勻速運動目標不同解算周期誤差

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

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

3.3 不同解算精度的對比

求解過程中迭代誤差設定直接影響解算精度，當設定0.1 m時，解算點精度在0.1 m內，當設定0.5 m時，解算點解算精度在0.5 m內，而在整個解算周期內，解算誤差仿真結果如圖10所示。

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

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

3.4 解算時間

按目標靜止，設定距離，解算時間如表1所示。

表1 不同距離解算時間

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

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

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

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

5 結論

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

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

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