0 引　言

随着深海资源的开采探测越来越频繁，船舶动力定位（DP）系统越来越受到重视。DP主要通过推进装置抵御外界风、浪、流的干扰，实现海上勘探、钻井、布管等作业。海浪干扰可分为一阶波浪引起的高频分量和二阶波浪引起的低频分量。其中高频分量使船舶推进器产生不必要的机械磨损和燃料消耗，需要进行滤除，因此动力定位系统滤波器的研究一直是一个热点^[1]。

20世纪60年代第1代动力定位系统采用了陷波滤波器，其结构简单，易于实现，但是会产生相位滞后，使系统稳定性变差。80年代后，第2代动力定位系统采用卡尔曼滤波器 ^[2]，其设计基于船舶模型，实时更新滤波参数，不会产生相位滞后，滤波效果比陷波滤波器更好。但卡尔曼滤波算法是一种线性算法，针对复杂海况下非线性的船舶系统滤波效果并不理想。为解决线性化问题，采用扩展卡尔曼滤波器，按首摇角度均分为36个工作点，进行分段线性化，但是需要调节的参数很多，过程噪声、测量噪声的协方差矩阵只能根据设计者的经验来确定，系统的稳定性无法得到保证。近些年，非线性滤波器算法的研究成为了重点^[3-7]。

无源性理论最早由Popov提出^[8]，有效解决了非线性问题，1974年Desoer从理论上给出了无源系统的稳定性证明^[9]。将无源性理论引入动力定位系统滤波器设计中，可以直接解决非线性问题，不需要进行模型线性化，滤波更精确；其滤波器参数可以根据海浪参数进行确定，调节简单，易于实现；其稳定性可以由李亚普诺夫稳定性判据进行证明，能实现滤波算法收敛稳定^[10]。

1 无源性理论

m输入l输出非线性系统如下式所示：

$S:\left\{ {\begin{aligned}& {\dot x = f(x,u)}{\text{，}}\\& {y = h(x,u)}{\text{。}}\end{aligned}} \right.$

(1)

其中： $x \in {R^n}$ 为状态变量； $u \in {R^m}$ 为系统输入变量； $y \in {R^l}$ 为系统输出变量；f为关于 $(x,u)$ 局部的Lipschitz函数。

无源性定义：如果存在半正定函数 $S(x)$ 和正定函数 $W(x)$ 满足对于任意初始状态 $x(0)$ 和时间 $\tau $ 公式（2）成立，则称系统是严格无源的。

$\int_0^\tau {{u^{\rm T}}} y{\rm d}t \geqslant S(x(\tau )) - S(x(0)) + \int_0^\tau {W(x){\rm d}t}{\text{，}}$

(2)

无源性证明：

KYP引理：设 ${{G}}(s) = {{C}}{(sI - {{A}})^{\rm T}}{{B}} +{{D}}$ 是传递函数矩阵，满足可控可观性，当且仅当存在矩阵 ${{P}} = {{{P}}^{\rm T}} > 0$ ，L，M和 $\varepsilon $ 满足

${{PA}} + {{A}^{\rm T}}{{P}}= - {{{L}}^{\rm T}}L - \varepsilon {{P}}{\text{，}}$

${{PB}}= {{{C}}^{\rm T}} - {{{L}}^{\rm T}}{{W}}{\text{，}}$

${{{W}}^{\rm T}}{{W}} ={{D}} + {{{D}}^{\rm T}}$

时， ${{G}}(s)$ 是严格正实的，进而无源的。

2 数学模型

无约束的船舶模型是六自由度的，如图1所示。

水面船舶可以只考虑纵荡、横荡、首摇等3个自由度运动。为描述船舶运动，建立地球固定坐标系 ${X_e}O{Y_e}$ 和随船坐标系 ${X_b}O{Y_b}$ ，如图2所示。

定义向量 ${\eta} = [x,y,\varphi ] \in {R^3}$ ，分别表示地球固定坐标系下的纵荡位移、横荡位移、首摇角度；向量 $\upsilon = [u,\nu ,r] \in {R^3}$ ，分别表示随船坐标系下的纵荡速度、横荡速度、首摇角速度。满足

$\dot{\eta} = {{J}}(\varphi ){\upsilon}{\text{，}} $

(3)

其中 ${{J}}(\varphi )$ 为转换矩阵，满足

${{J}}(\varphi ) = \left[ {\begin{array}{*{20}{c}}{\cos \varphi }&{ - \sin \varphi }&0\\{\sin \varphi }&{\cos \varphi }&0\\0&0&1\end{array}} \right]{\text{。}}$

(4)

船舶动力学模型如式（5）–式（8）所示：

$\dot{\xi} = {\varOmega}{\xi}+ {\varSigma} { w}{\text{，}}$

(5)

$\dot b = - {{{T}}^{ - 1}}b + {{ w}_b}{\text{，}}$

(6)

${{M}}\dot{\upsilon} + {{D}}{\upsilon} = {\tau}+ {{J}}(\varphi )b + {{ w}_l}{\text{，}}$

(7)

$y ={\eta}+ {\varGamma}{\xi}+ {{ w}_y}{\text{。}}$

(8)

其中： ${ M}\in {R^{3 \times 3}}$ 为系统惯性矩阵； ${ D}\in {R^{3 \times 3}}$ 表示水动力阻尼系数矩阵； $\tau \in {R^3}$ 表示控制力矩矩阵； $b \in {R^3}$ 表示缓慢变化的误差，包括风、流、二阶波浪，以及其他未建模动态误差； ${{ w}_l} \in {R^3}$ 为过程噪声向量。 ${ T} \in {R^{3 \times 3}}$ 为三维对角矩阵，表示时间常数； ${{ w}_b},{ w} \in {R^3}$ 为环境高频白噪声向量。 ${\xi} \in {R^{6 \times 3}}$ 为高频位置速度向量； ${\varOmega}\in {R^{6 \times 6}}$ ， ${\varSigma}\in {R^{6 \times 1}}$ 为常数矩阵，与海浪特性相关。 ${{ w}_y} \in {R^3}$ 为测量噪声向量， ${\varGamma} \in {R^{3 \times 3}}$ 是常数矩阵。

3 无源观测器设计

$\mathop {\hat {\xi}}\limits^ \cdot = {\varOmega}\hat {\xi}+ {{ K}_1}\tilde y{\text{，}}$

(9)

$\mathop {\hat{\eta}}\limits^. ={ J}(\varphi )\dot {\upsilon}+ {{ K}_2}\tilde y{\text{，}}$

(10)

$\mathop {\hat b}\limits^. = - {{ T}^{ - 1}}\hat b + {{ K}_3}\tilde y{\text{，}}$

(11)

${ M}\mathop {\hat {\upsilon}}\limits^. = -{ D}\hat {\upsilon} + {{ J}^{\rm T}}(\varphi )\hat b + \tau + {{ K}_4}\tilde y{\text{，}}$

(12)

$\tilde y = \hat{\eta}+ {\hat{\eta}_w}{\text{。}}$

(13)

其中：K₁，K₂，K₃，K₄为调节参数； $\hat{\xi}$ ， $\hat b$ ， $\hat b$ ， $\hat \eta $ 分别为 ${\xi}$ ，v，b， ${\eta}$ 的估计值； $\tilde y = y - \hat y$ 为输出估计误差。

4 无源性和稳定性证明 4.1 滤波器误差动态特性

定义变量 ${\hat \eta _o} = {\left[ {\begin{array}{*{20}{c}}{{{\hat {\xi}}^{\rm T}}}&{{{\hat{\eta}}^{\rm T}}}\end{array}} \right]^{\rm T}}$ ，得到式（14）–式（15）：

$\mathop {{{\hat \eta }_o}}\limits^ \cdot = {A_o}{\hat \eta _o} + {B_o}{ J}(\varphi )\hat{\upsilon}+ {K_o}\tilde y{\text{，}}$

(14)

$\hat y = {C_o}{\hat \eta _o}{\text{。}}$

(15)

其中：

${A_o} = \left[ {\begin{array}{*{20}{c}}\varOmega &0\\0&0\end{array}} \right]$ ， ${B_o} = \left[ {\begin{array}{*{20}{c}}0\\I\end{array}} \right]$ ， ${C_o} = \left[ {\begin{array}{*{20}{c}}\varGamma &I\end{array}} \right]$ ， ${K_o} = {\left[ {\begin{array}{*{20}{c}}{{K_1}}&{{K_2}}\end{array}} \right]^{\rm{T}}}$ 。

定义滤波器估计误差 $\tilde{\upsilon}={\upsilon} - \hat {\upsilon} $ ， $\tilde b = b - \hat b$ ， ${\tilde \eta _o} = {\eta _o} - {\hat \eta _o}$ ，得到滤波器动态误差如式（16）–式（17）所示：

$\mathop {\tilde x}\limits^. ={ A}\tilde x +{{BJ}}(\varphi )\tilde{\upsilon}{\text{，}}$

(16)

$\tilde z = { C}\tilde x{\text{。}}$

(17)

其中 $\tilde x = {\left[ {\begin{array}{*{20}{c}}{{{\tilde \eta }_o}}&{\tilde b}\end{array}} \right]^{\rm T}}$ ， $\tilde z = {K_4}\tilde y - \tilde b$ ，

${ A} = \left[\!\!\!\!{\begin{array}{*{20}{c}}{{A_o} - {K_o}{C_o}}&0\\{ - {K_3}{C_o}}&{ - {T^{ - 1}}}\end{array}}\!\!\!\!\right]$ ， ${ B} = \left[\!\!{\begin{array}{*{20}{c}}{{B_o}}\\0\end{array}}\!\!\right]$ ， ${ C} = {K_4}{C_o} - I$ 。定义新的误差动态变量： ${\varepsilon _z} = - {{ J}^{\rm T}}(\varphi )\tilde z$ ， ${\varepsilon _\upsilon } ={ J}(\varphi )\tilde \upsilon $ 。

滤波器误差动态特性如图3所示。

4.2 子系统H₁无源性分析

选取正定函数V满足

$V = \frac{1}{2}{\tilde {\upsilon}^{\rm T}}{ M}\tilde{\upsilon}{\text{，}}$

(18)

对时间求导得：

$\dot V = - \frac{1}{2}{\tilde {\upsilon}^{\rm T}}({ D} + {{ D}^{\rm T}})\tilde{\upsilon} - {\tilde z^{\rm T}}{ J}(\varphi )\tilde{\upsilon} {\text{，}}$

(19)

代入 ${\varepsilon _z} = - {{ J}^{\rm T}}(\varphi ){\tilde z_o}$ ，得到：

$\dot V = - \frac{1}{2}{\tilde {\upsilon}^{\rm T}}({ D} + {{ D}^{\rm T}})\tilde {\upsilon}+ {\tilde{\upsilon} ^{\rm T}}{\varepsilon _z}{\text{，}}$

(20)

进行积分得到：

$\int_{{t_0}}^t {\varepsilon _z^{\rm T}} \tilde{\upsilon} (\tau ){\rm d}\tau \geqslant \alpha {\tilde {\upsilon}^{\rm T}}\tilde{\upsilon}+ \beta {\text{。}}$

(21)

其中：

$\alpha = \displaystyle\frac{1}{2}{\lambda _{\min }}({ M}) > 0$ ， $\beta = \displaystyle\frac{1}{2}\int_{{t_0}}^t {{{\tilde{\upsilon}}^{\rm T}}({ D} + {{ D}^{\rm T}})} \tilde{\upsilon}{\rm d}\tau \geqslant 0$

由无源性定义式（2）可知子系统H₁是严格无源的。

4.3 子系统H₂无源性分析

选取合适的增益矩阵 ${K_1},{K_2},{K_3},{K_4}$ 使H₂满足KYP引理，实现子系统无源性。

设 ${K_1},{K_2},{K_3},{K_4}$ 满足：

${K_1} = \left[ {\begin{array}{*{20}{c}}{diag({k_{1i}})}\\{diag({k_{2i}})}\end{array}} \right]$ ， ${K_2} = diag({k_{3i}})$ ， ${K_3} = diag({k_{4i}})$ ， ${K_4} = diag({k_{5i}})$ ， $\quad i = 1,2,3$ 。

由于 ${K_3},{K_4}$ 为对角矩阵，所以H₂可以表示为解耦传递函数：

$\tilde z(s) = {H_o}(s){H_B}(s){\varepsilon _\upsilon }(s){\text{，}}$

(22)

其中 ${H_o}(s) = {C_o}{(sI + {A_o} - {K_o}{C_o})^{ - 1}}{B_o}$ ，

${H_B}(s) = {K_4} + {(sI + {T^{ - 1}})^{ - 1}}{K_3}$ 。

${H_o}(s)$ 对角元素 $h_o^i(s)$ 满足：

$h_o^i(s) \!\!=\!\!\! \displaystyle\frac{{{s^2} \!+\! 2{\xi _i}{\omega _{oi}}s \!+\! \omega _{oi}^2}}{{{s^3} \!\!+\!\! (\!{k_{2i}} \!\!+\!\! {k_{3i}} \!\!+\!\! 2{\xi _i}{\omega _{oi}}\!){s^2} \!\!+\!\! (\!\omega _{oi}^2 \!\!+\!\! 2{\xi _i}{\omega _{oi}}{k_{3i}} \!\!-\!\! {k_{1i}}\omega _{oi}^2\!)s \!\!+\!\! \omega _{oi}^2{k_{3i}}}}$ 。采用陷波滤波器结构来进行滤波：

$h_{od}^i = \frac{{{s^2} + 2{\xi _i}{\omega _{oi}}s + \omega _{oi}^2}}{{({s^2} + 2{\xi _{ni}}{\omega _{oi}}s + \omega _{oi}^2)(s + {\omega _{ci}})}}{\text{，}}$

令 $h_o^i \equiv h_{od}^i$ 得：

${k_{1i}} = - 2{\omega _{ci}}({\xi _{ni}} - {\xi _i})/{\omega _{oi}}{\text{，}}$

${k_{2i}} = 2{\omega _{oi}}({\xi _{ni}} - {\xi _i}),{k_{3i}} = {\omega _{ci}}{\text{，}}$

$1/T \ll {k_{5i}}/{k_{4i}} < {\omega _{oi}} < {\omega _{ci}}{\text{。}}$

其中： ${\omega _{ci}}$ 为滤波器截止频率； ${\omega _{oi}}$ 为高频波浪峰值频率； ${\zeta _{ni}}$ 为相对阻尼系数； ${\xi _i}$ 为阻尼系数；T为环境干扰力时间常数。

由复合系统无源性理论可知，子系统H₁H₂具有无源性，总系统无源。

4.4 稳定性证明

选取Lyaponov函数V满足：

$V = {\tilde{\upsilon}^{\rm T}}{ M}\tilde{\upsilon}+ {\tilde{ x}^{\rm T}}{ P}\tilde { x}{\text{，}}$

(23)

两边求导并代入KYP引理公式得：

$\dot V = - {\tilde {\upsilon}^{\rm T}}({ D} + {{ D}^{\rm T}})\tilde{\upsilon} - {\tilde { x}^{\rm T}}{ Q}\tilde{ x}{\text{。}}$

(24)

根据Lyaponov稳定性判据可知，系统稳定。

5 实验仿真

本文以动力定位船舶Northern Clipper为模型进行仿真验证^[12]，船长L=76.2 m，重量 $m = 4.591 \times {10^6}\;{\rm kg}$ ，船舶参数：

${ M} = \left[ {\begin{array}{*{20}{c}}{5.3122 \times {{10}^6}}&0&0\\0&{8.2831 \times {{10}^6}}&0\\0&0&{3.7454 \times {{10}^9}}\end{array}} \right]{\text{，}}$

${ D}= \left[ {\begin{array}{*{20}{c}}{5.0242 \times {{10}^4}}&0&0\\0&{2.7229 \times {{10}^5}}&{ - 4.3933 \times {{10}^6}}\\0&{ - 4.3933 \times {{10}^6}}&{4.1894 \times {{10}^8}}\end{array}} \right]{\text{。}}$

时间常数 $T = diag\left[ {\begin{array}{*{20}{c}}{1000}&{1000}&{1000}\end{array}} \right]$ 。

5.1 四种海况高频波浪仿真

根据国际海事组织(IMO)定义4种典型海况参数如表1所示。4种海况能量频谱函数如图4所示。由仿真结果可知，当海况由平静变为恶劣时，高频波浪峰值频率逐渐减小，幅值逐渐增大，符合实际海况情况。

5.2 卡尔曼滤波器仿真

卡尔曼仿真原理可以参见文献[11]。仿真时间500 s。在4种海况下仿真结果如图5所示。黑色曲线为实测值，浅色曲线为滤波值。分析可得，当海况相对平稳时（海况1和海况2）卡尔曼滤波器能够有效滤除高频波浪噪声，滤波曲线平缓；但当海况复杂时（海况3和海况4），卡尔曼滤波器滤波效果变差，仍然有高频分量未被滤除。

5.3 非线性无源滤波器仿真

仿真时间500 s，4种海况下仿真结果如图6所示。黑色曲线为实测值，浅色曲线为滤波值。4种海况下滤波器的陷波作用如图7所示，结果显示4种海况下滤波器的陷波频率和高频波浪噪声的峰值频率均相近，能够有效的滤除波浪高频成分。分析可得，无源非线性滤波器在4种海况下均能够有效滤除高频分量，对比图5可得，当海况平稳时，2种滤波器滤波效果近似，但当高海况和恶劣海况时，无源非线性滤波器效果明显更好。

6 结　语

针对动力定位系统中高频环境干扰的问题，设计了无源非线性滤波器，并与卡尔曼滤波器分别在不同海况下进行了滤波性能对比分析，结果显示，无源非线性滤波器滤波效果更强。