水下作业环境一直以来都是各国关注的重点，近年来国内外对无人水下潜航器（unmanned underwater vehicle，UUV）技术进行了广泛研究与应用^[1]. UUV作为典型的水下作业工具，可以完成水下情报监视、侦察、反潜、探测、识别、海洋调查、通信及载荷输送等不同的作业任务^[2]. 目前研制的UUV大多是以电池动力为主，综合考虑成本、放电稳定性、可维护性及安全性，应用锂电池尤为典型^[3-4]. UUV因作业时间受能量制约，适合选用高能量密度电池作为动力源（如常用的三元锂离子电池），无人系统对状态反馈的精度要求高，其中荷电状态（state of charge，SOC）精度易受其他状态的衰减干扰，因此对系统进行状态联合估计必不可少.

UUV精确的电池SOC估计可以提高能量利用率和设备续航时间敏感性^[5]. 工程上常用安时积分法进行SOC粗略估计，不仅依赖传感器和初状态精度，还存在累积误差^[6]. 基于模型的估计方法应用较多，能以较高精度描述特征，所需的数据量较少^[7-8]；鉴于模型内部的非线性关系，应用扩展卡尔曼滤波（extended Kalman filter，EKF）进行SOC估计具有较高精度^[9-11]. 随着充放电的进行，电池老化使得容量和电压关系发生变化. 若不修正，则设备将得到错误反馈. 电池健康状态（state of health，SOH）很好地描述了老化程度，且和SOC之间存在耦合关系，可以利用SOH修正SOC估计精度. SOH反馈主要依赖测试直接获取或模型间接获取，前者须基于实验室环境，因UUV不方便经常拆卸，在实际工程应用中受限；后者可以通过外特性数据来实时估算，具有较高的可实施性，常用方法有差分电压分析^[12]、容量增量分析^[13]、充电数据分析^[14]等，结合高斯过程回归^[15]、神经网络^[16]算法进行估计，但根据设备实际作业工况选取表征十分重要. 王志福等^[17]提出深度置信网络和自适应EKF/H_∞滤波融合算法在卷积神经网络模型下的SOC-SOH联合估计，该方法具有较好的精度和鲁棒性，但算法融合使计算复杂度增加，实际工程应用依赖硬件支撑，且多模块应用对UUV电池舱存在挑战. 赵靖英等^[6]提出基于模型和分数阶理论的SOC-SOH联合估计，在低频计算中能以少参数和高精度进行估计，但十分依赖模型精度，且对UUV内部成组使用的电池系统来说，不一致性可能带来较大的估计误差.

为了解决估计误差及工程实现的问题，提高设备状态反馈的精确性，本文提出适用于UUV电池系统的SOC-SOH联合估计方法. 采用遗忘因子递推最小二乘算法（forgetting factor recursive least squares，FF-RLS）在线辨识模型参数，结合EKF算法进行SOC估计. 通过多维表征结合改进支持向量回归（support vector regression，SVR）进行SOH估计；利用状态间耦合关系联合修正估计，提高状态反馈精度，通过实验验证方法的有效性.

三元电池具有高能量密度和较高的充电效率，非常适用于UUV对尺寸结构的要求，在UUV中应用较多. 其中国外部分UUV应用锂电池的情况如表1所示^[4].

将额定容量为10 A·h的三元镍钴锰酸锂电池分为4组进行实验，获取不同老化状态下的各项参数，记录循环的充放电数据，包括电流、电压、容量等数据. 开展特性测试，包括开路电压（open circuit voltage, OCV）、混合功率脉冲特性(hybrid pulse power characteristic, HPPC)、动态应力测试(dynamic stress test, DST)和美国联邦城市运行工况(federal urban driving schedule, FUDS)；搭建台架，以获取不同老化状态下的电池数据，如图1所示.

将采集步长设置为0.1 s，对4组电池均进行全寿命周期下的老化试验，获取电池全寿命周期下的外特性数据.

等效电路模型（equivalent-circuit model，ECM）可以使用简单元件模拟电池输出特性，适用于对外特性的研究. 一般ECM中串联RC网络越多，模型精度越高，但建模所需的参数增加，计算复杂度越高. 为了较好地满足电池电化学极化、浓差极化不同时间常数的极化过程，建立二阶ECM，如图2所示.

图中，U_oc为开路电压，U_t为端电压，U₁ 、U₂为2个极化电压，R_s为欧姆内阻，I_B为干路电流，R₁ 、R₂为2个极化内阻，C₁ 、C₂为2个极化电容. 由基尔霍夫电压定律，易得

(1)$ \frac{{{\mathrm{d}}{U_i}\left( t \right)}}{{{\mathrm{d}}t}} = \frac{{{I_{\mathrm{B}}}\left( t \right)}}{{{C_i}}} - \frac{{{U_i}\left( t \right)}}{{{R_i}{C_i}}} ; \;i = 1,2.$

(2)$ {U_{\mathrm{t}}}\left( t \right) = {U_{{\mathrm{oc}}}}\left( t \right) - \sum\limits_{i = 1}^2 {{U_i}\left( t \right)} - {I_{\mathrm{B}}}\left( t \right) {R_{\mathrm{s}}} . $

通过OCV测试获取当前状态下的SOC-OCV关系，采用7阶多项式拟合，结果如图3所示.

建立相关状态空间方程的离散形式如下：

(3)$ \begin{gathered} \left[ \begin{gathered} {U_1}\left( k \right) \\ {U_2}\left( k \right) \\ {\mathrm{SOC}}(k) \\ \end{gathered} \right] = \left[ {\begin{array}{*{20}{c}} {\exp \;( - \Delta t/({R_1}{C_1}))}& 0& 0 \\ 0& {\exp \;( - \Delta t/({R_2}{C_2}))}& 0 \\ 0& 0& 1 \end{array}} \right] \times \\ \left[ \begin{gathered} {U_1}\left( {k - 1} \right) \\ {U_2}\left( {k - 1} \right) \\ {\mathrm{SOC}}(k - 1) \\ \end{gathered} \right]+\left[ {\begin{array}{*{20}{c}} {{R_1}(1 - \exp \;( - \Delta t/({R_1}{C_1})))} \\ {{R_2}(1 - \exp \;( - \Delta t/({R_2}{C_2})))} \\ {\Delta t/{C_{\mathrm{n}}}} \end{array}} \right]{I_{\mathrm{B}}}(k) , \\\end{gathered} $

(4)$ {U_{\mathrm{t}}}\left( k \right) = {U_{{\mathrm{oc}}}}\left( {{\mathrm{SOC}}(k)} \right) - {U_1}\left( k \right) - {U_2}\left( k \right) - {I_{\mathrm{B}}}(k){R_{\mathrm{s}}} . $

式中：k为离散点，Δt为采样间隔，C_n为电池容量.

结合采集的电池数据，采用RLS进行参数辨识. 将ECM转化为ARMA模型，将式（2）经过拉普拉斯变换为

(5)$ {U_{\mathrm{t}}}\left( s \right) = {U_{{\mathrm{oc}}}}\left( s \right) - {I_{\mathrm{B}}}\left( s \right)\left( {{R_{\mathrm{s}}}+\sum\limits_{i = 1}^2 {\frac{{{R_i}}}{{1+s{R_i}{C_i}}}} } \right) . $

模型的传递函数为

(6)$ G\left( s \right) = \frac{{{U_{{\mathrm{oc}}}}\left( s \right) - {U_{\mathrm{t}}}\left( s \right)}}{{{I_{\mathrm{B}}}\left( s \right)}} = {R_{\mathrm{s}}}+\sum\limits_{i = 1}^2 {\frac{{{R_i}}}{{1+s{R_i}{C_i}}}} . $

令$ {\tau _i} = {R_i}{C_i} $（i = 1,2），$ {E_{\mathrm{t}}}\left( s \right) = {U_{{\mathrm{oc}}}}\left( s \right) - {U_{\mathrm{t}}}\left( s \right) $，由双线性变换原理，可得z平面离散传递函数为

(7)$ {Z_{\mathrm{T}}}(z) = \frac{{{\beta _0}+{\beta _1}{z^{ - 1}}+{\beta _2}{z^{ - 2}}}}{{1+{\alpha _1}{z^{ - 1}}+{\alpha _2}{z^{ - 2}}}} . $

式中：$ {\alpha _1} $、${\alpha _2} $、$\beta _0 $、$\beta _1 $、$\beta _2 $为待定系数. 离散方程为

(8)$ \begin{split} {E_{\mathrm{t}}}(k) =& {\alpha _1}{E_{\mathrm{t}}}(k - 1)+{\alpha _2}{E_{\mathrm{t}}}(k - 2) - \\& {\beta _0}{I_{\mathrm{B}}}(k) - {\beta _1}{I_{\mathrm{B}}}(k - 1) - {\beta _2}{I_{\mathrm{B}}}(k - 2) \end{split} . $

令

(9)$ \left. {\begin{array}{*{20}{c}} {\bf\textit{φ}} \left( {{k}} \right){\text{ = }}\left[ {{E_{\mathrm{t}}}(k - 1),{E_{\mathrm{t}}}(k - 2),{I_{\mathrm{B}}}(k),{I_{\mathrm{B}}}(k - 1),{I_{\mathrm{B}}}(k - 2)} \right] ,\\ {{\boldsymbol{\theta}} = [{\alpha _1},{\alpha _2}, - {\beta _0}, - {\beta _1}, - {\beta _2}]} ,\end{array}} \right\} $

则式（8）可以转化为FF-RLS可辨别的形式：

(10)$ {E_{\mathrm{t}}}(k) = {\bf\textit{φ}} ({{k}}){{\boldsymbol{\theta}} ^{\mathrm{T}}}+e(k) . $

模型参数方程组为

(11)$ \left. {\begin{array}{*{20}{l}} {{R_{\mathrm{s}}} = \dfrac{{{\beta _0} - {\beta _1}+{\beta _2}}}{{1 - {\alpha _1}+{\alpha _2}}}} ,\\ {{R_{\mathrm{s}}}+{R_1}+{R_2} = \dfrac{{{\beta _0}+{\beta _1}+{\beta _2}}}{{1+{\alpha _1}+{\alpha _2}}}}, \\ {\tau _1}{\tau _2} = \dfrac{{\Delta {t^2}(1 - {\alpha _1}+{\alpha _2})}}{{4(1+{\alpha _1}+{\alpha _2})}}, \\ {\tau _1}+{\tau _2} = \dfrac{{1 - {\alpha _2}}}{{\Delta t(1+{\alpha _1}+{\alpha _2})}}, \\ {{R_{\mathrm{s}}}\left( {{\tau _1}+{\tau _2}} \right)+{R_1}{\tau _2}+{R_2}{\tau _1} = \dfrac{{\Delta t({\beta _0} - {\beta _2})}}{{1+{\alpha _1}+{\alpha _2}}}} .\end{array}} \right\} $

解耦方程组，可得R_s、R₁、R₂、C₁、C₂.

SOH反映了电池老化程度，是缓变量. 为了探究老化过程中与SOH相关度高的外在表征，开展深度老化实验. 鉴于电池实际工作过程中放电工况无法预知，但充电工况相对平稳固定，且单个循环内可以认为SOH基本不变，因此可以通过实验结果及理论分析从充电工况数据中提取与老化相关的表征.

随着老化的加剧，电池内阻增大，恒电流在相同时间内充入容量相等，但电压变化增大. 分析不同老化程度下的充电数据，结果如图4所示. 图中，U为电压.

由实验结果可知，随着充放电循环的进行，相同充电时间内的电压变化幅度增大，意味着充电使电池变化同电压差的时间逐渐减少，充电时间和电池容量存在正相关性. 考虑到实际工程应用中初始充电电压的不确定性，应选择较大的起始电压到恒流截止电压的充电时间为第1个表征因子(H1).

随着电池老化的加剧，电池恒压充电时间逐渐增加^[18]. 分析不同老化程度下的电池恒压阶段充电数据，结果如图5所示. 图中，I为充电电流.

由实验结果可知，由于充电截止电流相同，充电时间增加意味着老化后期的电池在恒压阶段的充电电流下降速率减缓，相同时间内电流均值和电池容量存在负相关性. 考虑到实际工程应用中可能出现的充电周期不完整性，选取恒压充电开始后一小段时间内的电流平均值作为第2个表征因子(H2).

随着电池老化的加剧，容量和电压关系会发生变化，分析电池容量增量（incremental capacity, IC）是可行的方法^[19]. 分析电池不同老化程度下的IC与电压的对应关系，结果如图6所示. 图中，Q为充电过程中的电池容量.

由实验结果可知，IC曲线具有独特的形状、高度和位置，它反映了电池充放电过程中的电化学反应. 其中存在2个明显的高峰，表示当前电压下电池容量随电压变化的斜率，有研究表明峰值的变化可能与电池的活性材料损失有关. 随着老化的加剧，两峰高度逐渐降低，对应的电压位置逐渐增大，因此可以从以上变化趋势中提取相关表征^[20]. 由于两峰有明显的跟随性质，为了避免多重共线性影响，考虑提取第2峰的高度、位置作为第3、4个表征因子（H3、H4）.

为了验证表征因子的合理性，对因子进行相关性分析. 结合标准组1#电池的测试结果，对比如表2所示.

由实验电池表征因子与电池SOH的相关性可知，本文所提的因子可以对电池SOH进行表征. 考虑到所提的因子和相关性验证都是基于本文实验电池的数据，无法说明方法具备普适性. 为了验证表征因子的泛用型，联立美国国家航空航天局卓越预测中心（National Aeronautics and Space Administration，NASA）公开数据集中的2 A·h的18650型电池进行验证，结果如图7所示.

以上分析可知，所提表征与电池容量间具有良好的相关性，基于笔者^[21]先前的研究，联立前文所提的4种因子，建立基于四维输入的SVR模型. 采用粒子群优化算法（particle swarm optimization, PSO）自适应确定参数，选取2个实验标准组数据进行SOH估计，具体流程如下. 1）结合采集的电池数据，选取当前循环状态下的四维因子（H1、H2、H3、H4）. 2）将四维因子作为输入，实时通过PSO优化SVR模型参数. 3）输出当前循环下的SOH估计结果. 估计结果如图8所示.

对以上2组电池结果进行估计精度计算，均方根误差（RMSE）分别为0.003 9、0.006 8，平均绝对百分比误差（MAPE）分别为0. 351%、0. 6108%，可以证明所建立的SOH估计模型具有较好的效果.

结合式（3）、（4），令

(12)$ \left. {\begin{array}{*{20}{c}} {{\boldsymbol{A}} = \left[ {\begin{array}{*{20}{c}} {\exp \;( - \Delta t/({R_1}{C_1}))}&0&0 \\ 0&{\exp \;( - \Delta t/({R_2}{C_2}))}&0 \\ 0&0&1 \end{array}} \right]} ,\\ {\begin{array}{*{20}{c}} {{\boldsymbol{x}}\left( {{k}} \right) = \left[ \begin{gathered} {U_1}\left( k \right) \\ {U_2}\left( k \right) \\ {\mathrm{SOC}}(k) \\ \end{gathered} \right]},\;{{\boldsymbol{B}} = \left[ {\begin{array}{*{20}{c}} {{R_1}(1 - \exp \;( - \Delta t/({R_1}{C_1})))} \\ {{R_2}(1 - \exp \;( - \Delta t/({R_2}{C_2})))} \\ {\Delta t/{C_{\mathrm{n}}}} \end{array}} \right]} .\end{array}} \end{array}} \right\} . $

SOC估计可按以下流程开展.

1） 参数初始化，设置初始协方差矩阵P₀、初始状态x(0)、初始模型参数.

2） 引入遗忘因子进行参数迭代更新.

3） 结合更新的模型参数，进行状态初步预测：

(13)$ \hat {\boldsymbol{x}} \left( k \right) = {\boldsymbol{A}} {\boldsymbol{x}}\left( {k - 1} \right)+{\boldsymbol{B}} {I_{\mathrm{B}}}\left( {k - 1} \right) . $

4） 预测协方差矩阵P(k)：

(14)$ \hat {\boldsymbol{P}} \left( k \right) = {\boldsymbol{A}} {\boldsymbol{P}}\left( {k - 1} \right) {{\boldsymbol{A}}^{\mathrm{T}}}+{\boldsymbol{Q}}\left( {k - 1} \right) . $

式中：Q为噪声矩阵.

5） 更新卡尔曼增益K(k)，其中C(k)为输出矩阵：

(15)$ \left. {\begin{array}{*{20}{l}} {{\boldsymbol{C}}\left( k \right) = \left[ {\begin{array}{*{20}{c}} 1&1&{\dfrac{{\delta {U_{{\mathrm{oc}}}}\left( {{\mathrm{SOC}}\left( k \right)} \right)}}{{\delta {\mathrm{SOC}}\left( k \right)}}} \end{array}} \right]}, \\ {{\boldsymbol{K}}\left( k \right) = \hat {\boldsymbol{P}} \left( k \right){\boldsymbol{C}}^{\mathrm{T}}{\left( k \right)}{{\left[ {{\boldsymbol{C}}\left( k \right)\hat {\boldsymbol{P}} \left( k \right){\boldsymbol{C}}^{\mathrm{T}}{\left( k \right)}+{{r}}\left( k \right)} \right]}^{ - 1}}} .\end{array}} \right\} $

式中：r(k)为测量噪声.

6） 状态更新：

(16)$ {\boldsymbol{x}}\left( k \right) = \hat {\boldsymbol{x}} \left( k \right)+{\boldsymbol{K}}\left( k \right) e\left( k \right) . $

式中：$ e\left( k \right) $为电压估计值与测量值的误差.

7） 更新协方差矩阵：

(17)$ {\boldsymbol{P}}\left( k \right) = \left( {{\boldsymbol{E}} - {\boldsymbol{K}}\left( k \right) {\boldsymbol{C}}\left( k \right)} \right) \hat {\boldsymbol{P}} \left( k \right) . $

式中：E为单位矩阵.

8） 令$ k = k+1 $，进入下一次循环.

随着电池老化的加剧，SOH不断降低，电池最大可用容量不断减小，这会使SOC估计的误差逐渐增大. 通过充放电老化试验数据，探究SOC-OCV^[22]在不同SOH下对应关系的改变，选取电池初始循环、第210次循环、第450次循环、第630次循环、第810次循环以及第1 000次循环数据进行分析，结果如图9所示.

从结果可以看出，随着循环数的增加，SOC-OCV曲线特征明显改变. 为了提高电池的估计精度，对不同SOH下的SOC估计进行修正.

随着充放电循环的不断进行，电池容量将会出现衰减，这会导致在不同老化状态下的SOC估计出现误差. 为了提高状态反馈精度，提出SOC-SOH联合估计，流程如图10所示.

根据ECM和BMS采集到的数据，通过HPPC确定初始参数，采用FF-RLS在线辨识ECM参数. 选取H1、H2、H3、H4表征因子. 采用SVR建立四维输入的SOH估计模型，利用PSO优化选取模型参数. 建立基于EKF的SOC估计模型，由辨识的ECM参数实时更新状态空间方程矩阵，通过SOH修正SOC估计模型. 由2.3节所推导的估计流程可知，不同SOH状态下对应的SOC-OCV曲线特征会直接影响参数更新精度，结合3.1节探究可得SOH与SOC-OCV曲线特征相关，可由SOH实时更新当前状态下的SOC-OCV关系，修正估计结果. 具体如下：在不同的电池老化状态下，由SOH值修正曲线的7阶多项式参数，获取精确的SOC-OCV对应关系.

为了验证所提SOC-SOH联合估计方法的泛用性，对电池采用不同测试工况（DST、FUDS）的数据进行算法验证. 电池真实SOC初始为80%，根据不同的初始估计值对电池SOC进行估计，其中第630次循环下标准组电池的估计结果如图11所示. 图中，t_e为估计时间，SOC₀为SOC的初始估计值.

以SOC₀=70%、DST工况为例，修正后估计误差由2.47%降为0.92%，不同工况下修正后的收敛时间和精度都得到了优化，证明了算法的通用性.

为了验证方法的有效性，采用不同老化状态（第450次循环、第930次循环）数据进行验证，结果如图12所示.

从第450~930次循环的结果可以看出，不同老化状态的联合修正都提升了模型的精度，证明了所提方法在电池寿命周期内的有效性.

为了量化联合估计算法的修正精度，以系数RMSE和MAPE对第630次循环下的估计效果进行评估，结果如表3所示. 表中，RMSE₁、RMSE₂分别为第1、2组标准组电池的RMSE，MAPE₁、MAPE₂分别为第1、2组标准组电池的MAPE.

从表3可以看出，通过当前的SOH状态对SOC估计进行修正后，算法的收敛速度和精度都得到了提升. 在FUDS/DST测试工况下，2组标准组电池修正前、后的RMSE和MAPE都减小了，且无论初始值设置如何，都能够完成SOC估计.

以上是以初始额定容量定义的SOC，该定义不能表示电池在当前状态还能放出多少电量，须进行真实值转化. 以充满状态SOC为100%，可以进行如下转换：

(18)$ {{\mathrm{SOC}}_{{\text{real}}}} = 100{\text{%}} - \frac{{\left( {100{\text{%}}- {{\mathrm{SOC}}_{{\mathrm{est}}}}} \right) {C_{{\mathrm{n}}0}}}}{{{C_{{\mathrm{n}}0}} {\mathrm{SOH}}}} . $

式中：SOC_real为基于当前最大容量定义的SOC，C_n0为初始额定容量，SOC_est为估计到的基于初始容量定义的SOC. 通过转化可以得到电池当前阶段下的真实状态，提高设备的状态灵敏性.

（1）所提的4个状态表征因子与电池容量间有较强的相关性，2组标准电池SOH估计结果RMSE分别为0.003 9、0.006 8，MAPE分别为0.351%、0.610 8%，证明SOH估计模型具有较好的估计精度.

（2）所提的联合估计方法在不同测试工况下都有较好的效果，2组电池在630次循环FUDS工况下估计的RMSE由0.027 8、0.016 3变化为0.010 3、0.008 5，MAPE由3.58%、2.03%变化为1.07%、0.8%；在DST工况下的估计RMSE由0.016 5、0.018 4变化为0.014 1、0.006 4，MAPE由1.73%、2.34%变化为1%、0.73%，证明了方法的通用性.

（3）所提的联合估计方法在不同老化状态下的估计精度都得到了提升，证明了所提方法在电池寿命范围内有较高的可靠性. 与其他方法相比，所提的SOC-SOH联合估计模型可以通过多维的外在表征实时反馈出精确的UUV内部电池状态，且模型计算复杂度低，对硬件要求不高，可以作为有效的UUV嵌入式BMS算法进行应用.