張成舉， 王聰， 曹偉， 王金強(qiáng)

(哈爾濱工業(yè)大學(xué) 航天學(xué)院, 黑龍江 哈爾濱 150001)

0 引言

當(dāng)水下航行體運(yùn)動(dòng)速度到達(dá)50 m/s時(shí)，航行體頭部或整體會(huì)被空泡包裹，此時(shí)航行體稱為超空泡航行體。超空泡航行體高速運(yùn)動(dòng)時(shí)會(huì)出現(xiàn)兩種情況：一是航行體未被空泡全部包裹，航行體滑行力不為0；二是航行體被空泡全部包裹，航行體滑行力為0.海洋環(huán)境復(fù)雜多變，超空泡航行體受環(huán)境噪聲和測(cè)量噪聲的影響較大，因此超空泡航行體的運(yùn)動(dòng)控制具有極大的不確定性，所以超空泡航行體的運(yùn)動(dòng)控制是當(dāng)今的研究熱點(diǎn)。

在超空泡航行體運(yùn)動(dòng)控制方面，已有許多學(xué)者開展了研究，例如：Dzieiski等[1]建立了縱平面內(nèi)4狀態(tài)2自由度模型，對(duì)其基準(zhǔn)問題開展了研究；Vanek等[2]考慮超空泡航行體運(yùn)動(dòng)過程中的空泡存在記憶效應(yīng)，基于反饋線性化方法提出了雙閉環(huán)控制系統(tǒng)，證明了該控制系統(tǒng)的跟蹤控制有效性；Qiang等[3]改進(jìn)了反步法，提出狀態(tài)反饋控制，對(duì)超空泡航行體的運(yùn)動(dòng)穩(wěn)定性開展了研究；Nguyen等[4]考慮空泡形態(tài)變化會(huì)對(duì)超空泡航行體運(yùn)動(dòng)有重要影響，提出反饋控制，對(duì)超空泡航行體運(yùn)動(dòng)穩(wěn)定性開展了研究；Lü等[5]基于保代價(jià)理論提出自適應(yīng)滑?？刂破?，針對(duì)超空泡航行體存在非匹配不確定性問題開展了控制研究；Ruzzence等[6]運(yùn)用最優(yōu)方程獲取航行體的最優(yōu)彈道，對(duì)航行體轉(zhuǎn)彎運(yùn)動(dòng)和下潛運(yùn)動(dòng)開展了研究；Mao等[7]對(duì)超空泡航行體提出基于反步法的自適應(yīng)控制器，設(shè)計(jì)尾舵效率計(jì)算方法，利用李雅普諾夫方法證明了該方法的有效性；Zhao等[8]基于超空泡航行體空泡突變特性，對(duì)航行體變深控制進(jìn)行了研究；白濤等[9]提出基于混合卡爾曼濾波的變深運(yùn)動(dòng)控制算法，針對(duì)超空泡航行體的測(cè)量噪聲干擾開展了研究，采用高斯白噪聲測(cè)量了噪聲，仿真結(jié)果表明該控制系統(tǒng)提高了系統(tǒng)的運(yùn)動(dòng)穩(wěn)定性；龐愛平等[10]提出采用魯棒控制算法來保證超空泡航行體各個(gè)狀態(tài)變量在指定范圍內(nèi)變化，通過仿真證明該系統(tǒng)具有較好的穩(wěn)定性；呂瑞等[11]針對(duì)超空泡航行體提出基于增益自適應(yīng)變結(jié)構(gòu)控制算法，經(jīng)過仿真驗(yàn)證了該控制算法的有效性；王京華等[12]提出基于空泡記憶效應(yīng)的超空泡航行體控制算法，通過仿真驗(yàn)證了該控制算法的有效性；李洋等[13]提出針對(duì)超空泡航行體的反演控制算法，通過仿真驗(yàn)證了該控制系統(tǒng)的穩(wěn)定性。

根據(jù)以上文獻(xiàn)可知，采用濾波器同時(shí)降低環(huán)境干擾與測(cè)量干擾的研究較少。本文針對(duì)超空泡航行體受到環(huán)境噪聲和測(cè)量噪聲的影響，設(shè)計(jì)基于無跡卡爾曼濾波器(UKF)的最優(yōu)控制器，通過仿真分析驗(yàn)證了該控制器的有效性。

1 超空泡航行體動(dòng)力學(xué)模型

以超空泡航行體縱平面內(nèi)穩(wěn)定問題為研究對(duì)象，通過對(duì)航行體建模受力進(jìn)行分析，可知超空泡航行體主要受力如下：

1) 航行體所受重力G，重力方向指向地心。

2) 超空泡航行體發(fā)動(dòng)機(jī)的推力T，推力作用方向與航行體縱軸重合。

3) 空化器轉(zhuǎn)動(dòng)過程中為航行體提供的升力[1]：

(1)

式中：ρ為水的密度；Rn為空化器半徑；v為航行體航行速度；c0為空化器初始阻力系數(shù)；σ為航行體空化數(shù)；vy為航行體縱向速度；Lc為空化器頂端至航行體質(zhì)心之間的距離；ωz為航行體俯仰角速度；δc為空化器轉(zhuǎn)角。

4) 超空泡航行體航行過程中，其空泡區(qū)域分為全沾濕區(qū)、半沾濕區(qū)和空泡包裹區(qū)。航行體空泡區(qū)域劃分如圖1所示，其中δf為航行體尾翼轉(zhuǎn)角。

圖1 超空泡航行體空泡區(qū)域示意圖Fig.1 Configuration of supercavitating vehicle’s cavitation

空化器轉(zhuǎn)動(dòng)會(huì)使空泡軸線偏移，重力作用也會(huì)導(dǎo)致空泡上漂。對(duì)空泡形態(tài)預(yù)測(cè)利用較多的是Logvinovich[14]提出的空泡形態(tài)半理論和半經(jīng)驗(yàn)預(yù)測(cè)公式。本文采用此公式用于空泡形態(tài)預(yù)測(cè)，在文獻(xiàn)[11]提出的空泡狀態(tài)方程中，空泡截面為軸對(duì)稱橢圓，假設(shè)空泡在任意位置均為圓形，t時(shí)刻空泡半徑為

(2)

(3)

空泡受重力影響會(huì)導(dǎo)致空泡軸線發(fā)生偏移，在空化數(shù)σ<0.1且弗勞德數(shù)Fr較大的條件下，重力導(dǎo)致的空泡軸線偏移量為

(4)

式中：g為重力加速度；LRmax為空泡長度,LRmax=2Rn/(1.92/σ-3)；ζ=x/LRmax，x為空化器與空泡截面的距離。

空化器轉(zhuǎn)動(dòng)對(duì)空泡形態(tài)影響較大，并且航行體俯仰運(yùn)動(dòng)對(duì)空泡產(chǎn)生極大的非線性干擾，基于獨(dú)立膨脹原理，空化器轉(zhuǎn)動(dòng)導(dǎo)致空泡軸線偏移量為

hc=0.82(1+σ)(θ+δc)Rn(0.46-σ+ζ)，

(5)

式中：θ為航行體俯仰角。

當(dāng)超空泡航行體部分處于沾濕狀態(tài)時(shí)，航行體尾部會(huì)與空泡之間發(fā)生接觸和碰撞，航行體產(chǎn)生滑行力。根據(jù)細(xì)長體理論，將滑行力簡化為工程中計(jì)算細(xì)長體浸入自由液面，采用如下考慮重力和空化器轉(zhuǎn)動(dòng)的滑行力計(jì)算公式[2]：

(6)

式中：

(7)

(8)

R為航行體半徑，h′為沾濕深度，αp為航行體沾濕角，yc為空泡軸線與航行體軸線之間的縱向偏差，zc為空泡軸線與航行體軸線之間的側(cè)向偏差，ω為航行體旋轉(zhuǎn)角速度。

5) 尾舵包括俯仰舵和偏航舵，當(dāng)尾翼處于沾濕狀態(tài)時(shí)，尾舵才能發(fā)揮作用；航行體運(yùn)動(dòng)過程中，尾翼與空泡之間的關(guān)系處于不確定狀態(tài)；假定航行體在縱平面內(nèi)做小角度機(jī)動(dòng)，此時(shí)可認(rèn)為尾翼是一種楔形空化器，與空化器估算公式相同，則尾翼升力估算公式為

(9)

式中：η為尾翼升力與空化器升力的相似準(zhǔn)數(shù)，這里取0.5；Lf為尾翼與航行體質(zhì)心之間的距離。

建立動(dòng)量定理和動(dòng)量矩定理，可得超空泡航行體縱向運(yùn)動(dòng)模型矩陣：

(10)

(11)

I(t,τ)為尾舵沾濕率，

(12)

2 超空泡航行體最優(yōu)控制器設(shè)計(jì)

采用精確線性化方法，利用非線性反饋?zhàn)儞Q實(shí)現(xiàn)輸入輸出的精確線性化[15]，首先將方程(10)式改寫為

(13)

式中：

(14)

(15)

(16)

(17)

(18)

y1、y2為輸出變量y的分量，n1、n2為充分光滑的標(biāo)量函數(shù)。

對(duì)于此非線性多輸入多輸出(MIMO)系統(tǒng)，具有相對(duì)階向量(λ1,λ2)=(2,2)，總相對(duì)階[15]為λ1+λ2=4，這時(shí)不考慮內(nèi)動(dòng)態(tài)的穩(wěn)定性，輸出變量與輸入的關(guān)系為

(19)

式中：Lχ、LMi(i=1,2)分別為向量χ(x)和M(x)的李導(dǎo)數(shù)。

通過反饋可得

(20)

將(20)式代入(19)式，可得

(21)

將(21)式精確線性化后，可得超空泡航行體狀態(tài)方程：

(22)

(23)

式中：ε(t)為線性二次型最優(yōu)控制器；k1、k2、k3、k4為控制參數(shù)，其選擇應(yīng)使兩個(gè)解耦子系統(tǒng)對(duì)應(yīng)特征多項(xiàng)式的所有根均位于左半復(fù)平面；h1、h2、h3、h4為向量h(t)的分量，則可得誤差方程為

(24)

式中：

(25)

選取如下線性二次型性能指標(biāo)離散泛函數(shù)：

(26)

式中：Q為對(duì)稱正定狀態(tài)加權(quán)矩陣；R為對(duì)稱正定控制加權(quán)矩陣。

當(dāng)方程(26)式取得最小值時(shí)，可得線性二次型最優(yōu)控制反饋控制器為

ε(t)=-R-1ΘTPh(t)，

(27)

式中：P為滿足Riccati方程(28)式的正定對(duì)稱矩陣，

HTP+PH-PΘR-1ΘTP+Q=0.

(28)

綜上所述，聯(lián)立方程(20)式、方程(23)式和方程(27)式，可得基于精確線性化后的最優(yōu)控制器表達(dá)式為

(29)

3 帶有環(huán)境噪聲和測(cè)量噪聲的最優(yōu)控制算法仿真分析

令：超空泡航行體初始縱向速度vy=-4 m/s，俯仰角θ=3°，航行深度y=-4 m,俯仰角速度ωz=7°/s；正弦跟蹤方程為y(t)=-sint-5；航行體航行速度v=75 m/s，質(zhì)量m=22 kg，重力加速度g=9.81 m/s2,空化器半徑Rn=0.019 1 m，航行體柱段半徑R=0.050 8 m，航行體總長L=1.8 m，空化數(shù)σ=0.03，水的密度ρ=998 kg/m3,Lc=17L/28,Lf=11L/28.

采用經(jīng)輸入輸出精確線性化后的最優(yōu)控制器對(duì)超空泡航行體進(jìn)行位置跟蹤控制；環(huán)境噪聲和測(cè)量噪聲幅值均為0.02 m/s的高斯白噪聲[9]。考慮到較小的控制誤差，經(jīng)過多次仿真實(shí)驗(yàn)，為保證控制精度，選取控制參數(shù)k1=k2=k3=k4=20，

(30)

將Q和R值代入方程(29)式，滿足P為正定矩陣。采用以上參數(shù)進(jìn)行仿真，驗(yàn)證算法的有效性，仿真時(shí)間為20 s，仿真結(jié)果如圖2～圖8所示。

圖2 航行體深度變化(無濾波器)Fig.2 Change of depth (without filter)

圖3 航行體縱向速度變化(無濾波器)Fig.3 Change of vertical speed (without filter)

由圖2～圖8可知，超空泡航行體在環(huán)境噪聲和測(cè)量噪聲干擾下出現(xiàn)持續(xù)無規(guī)則振蕩。由于環(huán)境噪聲和測(cè)量噪聲的干擾，航行體滑行力呈現(xiàn)無規(guī)則振蕩，給航行體運(yùn)動(dòng)控制帶來了困難。

圖4 航行體俯仰角速度變化(無濾波器)Fig.4 Change of pitch angular velocity (without filter)

圖5 航行體俯仰角變化(無濾波器)Fig.5 Change of pitch angle (without filter)

圖6 航行體滑行力變化(無濾波器)Fig.6 Change of sliding force (without filter)

根據(jù)上述對(duì)超空泡航行體在環(huán)境噪聲和測(cè)量噪聲干擾下的運(yùn)動(dòng)狀態(tài)分析可知，超空泡航行體的運(yùn)動(dòng)是極其不穩(wěn)定的。因此，在超空泡運(yùn)動(dòng)控制中有必要加入卡爾曼濾波器，使航行體減少干擾噪聲的影響。

圖7 航行體空化器轉(zhuǎn)角變化(無濾波器)Fig.7 Change of cavitor angle (without filter)

圖8 航行體尾翼轉(zhuǎn)角變化(無濾波器)Fig.8 Change of tail angle (without filter)

4 基于UKF的最優(yōu)控制算法

本文針對(duì)以上超空泡航行體的運(yùn)動(dòng)不穩(wěn)定性，提出基于UKF的最優(yōu)控制算法，解決超空泡航行體的噪聲干擾問題。UKF采用卡爾曼線性濾波框架，對(duì)于一步預(yù)測(cè)方程，使用如下無跡變換處理均值和協(xié)方差的非線性傳遞問題：

(31)

式中：x(k+1)為k+1時(shí)刻的系統(tǒng)狀態(tài)；f(x(k),W(k))為非線性狀態(tài)方程函數(shù)，W(k)為環(huán)境噪聲；y(k)為觀測(cè)變量；h(x(k),V(k))為非線性觀測(cè)方程函數(shù)，V(k)為測(cè)量噪聲；W(t)為協(xié)方差S1的連續(xù)時(shí)間高斯白噪聲；V(t)為協(xié)方差S2的連續(xù)時(shí)間高斯白噪聲。

具體處理步驟如下：

1) 獲取超空泡航行體的一組采樣點(diǎn)，設(shè)置如下初值：

x(i)(k|k)=
[(k|k)(k|k)+T(k|k)-T]，

(32)

式中：x(i)(k|k)為k時(shí)刻超空泡航行體狀態(tài)值，i為第i個(gè)采樣點(diǎn)，i=1,2,…,2n，n為狀態(tài)維數(shù)；(k|k)為k時(shí)刻狀態(tài)估計(jì)值；為半正定矩陣，D(k|k)為狀態(tài)x(k|k)的方差，λ為縮放比例參數(shù)，λ=α2(n+a)-n，用于減小總的預(yù)測(cè)誤差；α為控制采樣點(diǎn)分布狀態(tài)的參數(shù)，用于調(diào)整Sigma點(diǎn)與均值點(diǎn)的距離。

2) 計(jì)算2n個(gè)Sigma點(diǎn)集的一步預(yù)測(cè)：

x(i)(k+1|k)=f[k,x(i)(k|k)]，

(33)

式中：x(i)(k+1|k)為超空泡航行體k+1時(shí)刻的預(yù)測(cè)狀態(tài)；f[k,x(i)(k|k)]為k時(shí)刻的非線性狀態(tài)方程。

3) 預(yù)測(cè)超空泡航行體的系統(tǒng)狀態(tài)量與方差：

(34)

式中：wm為采樣點(diǎn)均值的權(quán)值；wc為采樣點(diǎn)協(xié)方差的權(quán)值。

采樣點(diǎn)均值的權(quán)值和協(xié)方差的權(quán)值計(jì)算公式如下：

(35)

4) 根據(jù)一步預(yù)測(cè)值，再次運(yùn)用無跡變換產(chǎn)生新的Sigma點(diǎn)集：

(36)

5) 量測(cè)更新方程，得到Sigma點(diǎn)集的觀測(cè)預(yù)測(cè)值，通過加權(quán)求和得到超空泡航行體系統(tǒng)預(yù)測(cè)的均值、狀態(tài)量和觀測(cè)量之間的協(xié)方差、新息序列的協(xié)方差分別為

(37)

6) 計(jì)算卡爾曼增益矩陣：

(38)

7) 計(jì)算系統(tǒng)的狀態(tài)更新和協(xié)方差更新：

k(k+1|k+1)=k(k+1|k)+Kk(yk-k)，
D(k+1|k+1)=

(39)

5 基于UKF的最優(yōu)控制仿真分析

令：超空泡航行體初始縱向速度vy=-4 m/s，俯仰角θ=3°，航行深度y=-4 m,俯仰角速度ωz=7°/s；α通常選取一個(gè)較小的正實(shí)數(shù)(一般10-4≤α≤1)，α=0.01；a為待選參數(shù)，需要保證矩陣(n+λ)D(k|k)為半正定矩陣；κ=0；β用于合并方程中高階項(xiàng)的動(dòng)差，調(diào)節(jié)β可改善方差的精度，對(duì)于高斯分布，β=2時(shí)最優(yōu)。選取的控制器參數(shù)k1=k2=k3=k4=20，Q和R與第3節(jié)相同，對(duì)選取的深度信號(hào)y(t)=-sint-5進(jìn)行深度跟蹤，仿真時(shí)間20 s時(shí)得到仿真結(jié)果如圖9～圖15所示。

圖9 航行體深度變化(含濾波器)Fig.9 Change of depth (filter)

由圖9～圖14可知，通過運(yùn)用UKF對(duì)噪聲干擾處理，超空泡航行體處于穩(wěn)定運(yùn)動(dòng)狀態(tài)。

由圖9可知，經(jīng)過UKF處理后的跟蹤信號(hào)與理想信號(hào)相差較小，誤差能夠控制在0.1 m內(nèi)，而且未出現(xiàn)跟蹤延遲現(xiàn)象。

由圖10可知，超空泡航行體俯仰角速度急劇反向轉(zhuǎn)變，在4 s時(shí)趨于穩(wěn)定狀態(tài)，經(jīng)過濾波器處理后俯仰角速度跟蹤誤差控制較小。

圖10 航行體俯仰角速度變化(含濾波器)Fig.10 Change of pitch angular velocity (filter)

圖11 航行體縱向速度變化(含濾波器)Fig.11 Change of vertical speed (filter)

由圖11可知，超空泡航行體縱向速度變化穩(wěn)定，跟蹤誤差在0.1 m/s范圍內(nèi)。

由圖12可知，超空泡航行體俯仰角逐漸較小，變化率逐漸減小，在4 s時(shí)航行體俯仰角趨于穩(wěn)定，經(jīng)過UKF處理后的俯仰角跟蹤誤差控制較小。

圖12 航行體俯仰角變化(含濾波器)Fig.12 Change of pitch angle (filter)

圖13 航行體滑行力變化(含濾波器)Fig.13 Change of sliding force (filter)

圖14 航行體空化器轉(zhuǎn)角變化(含濾波器)Fig.14 Change of cavitor angle (filter)

由圖13可知，超空泡航行體起始狀態(tài)時(shí)航行體包裹面積較大，滑行力較小。開始運(yùn)動(dòng)后，滑行力呈現(xiàn)先增大、后減小現(xiàn)象，之后航行體滑行力減小為0，此時(shí)航行體除空化器外，均為空泡包裹，經(jīng)過UKF處理后的滑行力信號(hào)與理想狀態(tài)相差無幾。在航行體滑行力穩(wěn)定之前，誤差控制在0.1 kN范圍內(nèi)，但是從起始狀態(tài)到穩(wěn)定狀態(tài)，經(jīng)過UKF處理后的滑行力信號(hào)仍然比理想信號(hào)延遲。

由圖14和圖15可知，航行體空化器轉(zhuǎn)角與尾翼轉(zhuǎn)角偏轉(zhuǎn)范圍較小，與理想信號(hào)較為接近。

圖15 航行體尾翼轉(zhuǎn)角變化(含濾波器)Fig.15 Change of tail angle (filter)

綜合以上分析可知，通過運(yùn)用UKF對(duì)噪聲進(jìn)行有效處理，可使超空泡航行體達(dá)到穩(wěn)定運(yùn)動(dòng)狀態(tài)。

6 結(jié)論

本文以最優(yōu)控制為基礎(chǔ)，提出超空泡航行體UKF最優(yōu)控制策略，通過仿真實(shí)驗(yàn)得到如下結(jié)論：

1)在環(huán)境噪聲和測(cè)量噪聲干擾下，超空泡航行體運(yùn)動(dòng)極其不穩(wěn)定。

2)在UKF作用下，航行體運(yùn)動(dòng)狀態(tài)改變明顯，達(dá)到穩(wěn)定狀態(tài)所需時(shí)間減小。

3)將UKF引入狀態(tài)觀測(cè)器的設(shè)計(jì)，可以獲得環(huán)境噪聲和測(cè)量噪聲干擾下的最優(yōu)狀態(tài)估計(jì)量，從而實(shí)現(xiàn)航行體的高精度控制，減少輸入輸出線性化后帶來的誤差。