基于最大似然估计和混合梯度优化的射手模型辨识

吴骏雄， 林德福， 王辉， 袁亦方

(1.北京理工大学 无人机自主控制技术北京市重点实验室, 北京 100081；2.北京特种机电研究所, 北京 100012)

0 引言

光纤图像制导武器使用光纤图像与指令双向传输体制，可在复杂的交战条件下攻击坦克、舰船、地面或低速空中目标。弹载图像导引头能够将图像探测器敏感的战场实时图像传输到控制地面站，射手可根据显示器中的图像选择攻击目标以及攻击部位，并通过图像跟踪器自动跟踪，也可以通过控制手柄进行手动跟踪直至命中目标，即所谓的“人在回路”制导方式[1]。因此，了解射手对导引头跟踪误差的响应并对射手本身进行建模，不仅可以更深入地了解射手的感知过程，而且能够对制导控制系统设计[2]以及培训模拟设备评价起重要作用[3]。

对人进行建模的研究始于20世纪50年代[4]，其中最具代表性的是McRuer等提出的精确模型，该模型已在诸多领域应用[5]。随着现代控制理论的发展，Kleinman等[6]建立了人的最优控制模型，Hess[7]提出了能充分体现人的信号处理过程的Hess结构模型。近来，将智能技术应用到人的建模问题中逐渐成为研究热点，其中模糊控制和神经网络技术已经成功应用到飞行员模型辨识中[8-9]。

人响应的具体特征是由模型参数决定的[10]，可以利用参数辨识技术从实验数据中得到。目前人的模型参数辨识方法可分为两类：一类是针对人对时不变控制对象的响应；另一类是针对人对环境或者控制对象动力学突然变化时的响应[11]。傅里叶系数[12]和最大似然估计[13]是辨识时不变控制对象人的行为典型方法，前者需先估计非参数化的系统频率响应，再通过拟合非参数化模型得到模型参数；后者直接通过拟合时域数据来辨识模型参数。与傅里叶系数法相比，最大似然估计方法仅使用1个步骤来辨识模型，减小了参数估计的偏差和方差。针对时变控制对象人的响应辨识问题主要有最大似然、小波变换[11]和Kalman滤波[14]等方法。“人在回路”制导体制下射手的控制对象是导引头内回路，在飞行过程中其动力学特性几乎不变。因此，本文采用最大似然估计方法来辨识射手模型。

基于上述分析，本文拟针对射手控制行为进行建模辨识。目前国内外关于射手模型的文献较少，未有详细的关于射手建模的相关文献。其中赵军民等采用辅助变量法对模型参数进行辨识[2]，但其认为模型的神经延迟环节固定，实验条件也仅考虑在跟踪横向运动的坦克目标，未能充分激励射手的响应。因此，本文将设计具有较宽频率范围的激励信号，而且为了提高辨识精度，采用输出误差结构的最大似然法。此外，针对参数估计过程梯度优化算法对初值敏感且易发散的特性，采用混合梯度算法进行优化，即利用遗传算法逼近最优值，利用梯度算法的快速性加快算法运行，并使用单纯形法保证算法的稳定性。

1 参数估计问题

1.1 控制任务

光纤图像制导系统采用大容量光纤传输系统，能够将电视导引头的图像实时传输到地面站，且能在武器地面站实现传统的导引头图像跟踪能力，使得跟踪器可以重复使用。相比于弹载图像跟踪器，地面站图像跟踪器的功能可以设计得更加强大和全面，同时降低武器使用成本。人工控制模式下导引头的控制框图如图1所示。

图1中：qt为目标的惯性空间视线角；qs为导引头光轴相对于惯性空间的角度；Δq为视场角误差信号；Hp(s)为射手模型的传递函数，s为拉普拉斯算子；n为射手模型中的残差信号；u为射手的输出；K为信号放大器；G(s)和H(s)分别为稳定回路校正网络和速率陀螺。当目标偏离光轴时，射手通过显示屏观测到视场角误差Δq，并操控手柄使得光轴指向目标，实现对目标的跟踪。射手在导引头跟踪回路中相当于1个反馈控制器，输入信号为屏幕上显示的视场角误差。因此射手的目标跟踪任务可认为是单回路补偿任务，控制对象即为图1中的虚线框内部分。不失一般性，可假设G(s)≈1，H(s)≈1，且导引头稳定回路的带宽约为20～30 Hz，稳定回路的传递函数可近似为100/(s+100)[15]，则控制对象的传递函数Hc(s)可表示如下：

(1)

(1)式中的增益K决定了跟踪回路的带宽和稳定裕度。导引头跟踪回路带宽远小于稳定回路带宽，为保证导引头系统回路具有足够的稳定裕度，在这里选取增益K=2.

1.2 射手模型

本文采用精确模型[16]来描述射手行为。如图1所示，射手模型的响应包含传递函数Hp(s)的响应和残差信号n两部分，其中残差信号代表射手模型的非线性部分。射手对显示运动提示作出响应的线性部分形式如下：

(2)

式中：Kp为视觉感知的增益；τL为视觉感知超前时间常数；τl为视觉感知滞后时间常数；ωnm为系统自然频率；ξnm为系统阻尼；τd为视觉感知时间延时。 (2)式表明神经肌肉系统的动力学限制了射手响应，神经肌肉系统认为是1个2阶质量弹簧阻尼系统。

McRuer的交叉定理[16]表明，人会调整自身的配平动力学来适应控制对象的动力学，使得系统开环传递函数在交叉频率附近具有类似单积分器系统的频率特性。对于(1)式所示控制对象的动力学，射手需要在更高的频率处进行配平[17]。因此图1中射手传递函数表示如下：

(3)

综上所述，共有5个参数需要被估计，参数向量θ=[Kp,τL,τd,ξnm,ωnm]T。

2 射手模型辨识

下面基于最大似然估计方法来得到射手模型的参数，为简化整个参数估计过程，采用输出误差的结构来构造似然函数。为了解决非线性优化问题的全局最优问题，并避免出现局部最优现象，即找到观测值出现的最大概率，在使用经典的梯度优化算法之前，使用遗传算法来决定梯度算法的初值。

2.1 最大似然估计方法

最大似然估计是一种时域参数估计方法，在众多领域都有着应用[18]。最大似然估计是一类概率性贝叶斯估计方法，通过引入观测量的条件概率密度，构造1个以观测数据和未知参数为自变量的似然函数，以观测值出现的概率最大作为估计准则，获得系统模型的参数估计值。

在应用最大似然估计方法时，需要将射手传递函数模型表示为状态空间形式。在转化过程中，射手模型中的纯延时环节为非线性环节，可以将延时环节近似为1个高阶的传递函数，其中1种近似方法为Padé近似[19]，其具体形式如下：

(4)

式中：T为延迟时间；ka=0，…，λ，λ为Padé近似的阶次。为保证延时模型在高频处的精确描述，且不使近似射手模型阶次过高，本文采用5阶延时环节Padé近似[13]：

(5)

将(5)式代入(3)式中，可得近似的射手传递函数为

(6)

式中：b0，…，b6和a0，…，a6分别为对应的传递函数系数。

由于引入了Padé近似，射手传递函数的系数为未知辨识参数的非线性函数，例如(6)式中分母的第1项和最后1项系数为

(7)

将射手传递函数(6)式转换为状态空间的标准型，可得

(8)

式中：x为状态向量；状态矩阵A(θ)、控制矩阵B(θ)和输出矩阵C(θ)的表达式分别为

(9)

(8)式给出的射手模型仅在输出方程中包含1个额外的噪声项，即假设没有过程噪声的影响，从而可以采用输出误差结构的最大似然估计，大大简化了辨识步骤。在文献[20]关于辨识人的模型研究中，一般假设人的响应残差部分是1个零均值高斯噪声，并通过实验验证该假设。本文同样假设射手残差具有相同的特性，即认为(8)式的模型中残差是高斯量测噪声，其均值和方差定义如下：

(10)

式中：n(k)为残差信号n的离散采样时间序列，k为序列号，k=1,…,m，m为采样点个数。

最大似然估计方法中，需要寻找模型参数向量θ的估计值，该估计值使得似然函数L(θ)具有最大值。L(θ)定义为预测误差的条件概率密度函数：

L(θ)=f(v(1)，v(2)，…，v(k)，…，v(m)|θ),

(11)

式中：v(k)为预测误差，定义为量测的离散射手输出信号u(k)和离散射手模型输出信号(k)之差，

v(k)=u(k)-(k).

(12)

由(10)式定义的残差性质，可得似然函数为

(13)

图2所示为输出误差辨识策略的基本框架，该策略需要最小化1个二次型罚函数，该函数与输出误差相关，其中输出误差为射手实际输出和在相同输入情况下辨识模型的仿真输出之差。因此，罚函数定义为实际数据和仿真数据的匹配度，即参数估计值和实际估计值之间的差别。输出误差辨识策略采用迭代搜索算法在参数空间中寻找最优参数值，使得罚函数值最小，则最优参数值即为最终的参数估计值。

在应用最大似然估计方法时，若估计的参数使得似然函数最大，则为模型参数向量θ的估计值；而在输出误差策略下使用最大似然估计方法，需要用求取似然函数对数形式的最小值代替求取似然函数的最大值，从而简化为一个更加直接的优化问题。当对数似然取得全局最小时，得到的参数向量即为最大似然估计值，表示为ML. 因此，输出误差策略下的最大似然估计表示为

(14)

(14)式概括了本文射手模型的参数估计问题，它定义了1个强非线性优化问题。下面采用混合梯度优化算法来得到该优化问题的全局最优解。

2.2 混合梯度优化算法

传统求解最大似然极小值的方法有Levenberg-Marquardt、高斯- 牛顿和Newton-Raphson等方法。本文提出一种混合梯度优化算法，采用遗传算法和高斯- 牛顿梯度优化算法相结合的方式来估计射手传递函数的参数，具体策略如图3所示。

在没有先验信息情况下，该混合优化算法首先应用遗传算法进行初始迭代，然后用高斯- 牛顿和单纯形法算法来优化(14)式。由于较大的初始误差可能会导致参数估计发散，在计算初始时刻使用遗传算法以减小算法对初始值的敏感。此外，高斯- 牛顿算法具有较快的收敛速度，但当离真实解较远时会产生不准确的梯度信息、导致发散。因此，为增强算法的鲁棒性，当高斯- 牛顿算法发散时，切换为Nelder-Mead单纯形法进行迭代，直至高斯- 牛顿算法能够收敛。

2.2.1 遗传算法

遗传算法固有的随机性以及在初始参数向量集内的最优解开始搜索，导致该算法有很高的概率找到优化问题的全局最优值[21]。遗传算法首先需要创建1个初始种群，即初始参数向量集，该种群在参数的上下界随机选取，在此用矩阵表示为

(15)

通过目标函数(14)式评估种群中每个个体的适应度。图4所示为对当前种群每个个体运用3个基因操作从而产生新种群的过程，具体步骤如下：

1) 复制(elite)：从当前种群中选择适应度强的个体，并复制给下一代；

2) 交叉(crossover)：从当前种群中任选2个个体，随机选择基因中交叉点的位置，模拟进化过程中的繁殖现象，得到全新基因的个体；

3) 变异(mutation)：当前个体基因小概率随机选择变异点，并发生变化从而产生新的个体。

不断重复上述步骤，直至终止条件满足。

由于遗传算法使用概率搜索技术，导致1个搜索点到另1个搜索点的转移方式和关系具有不确定性，也导致每次迭代出的结果不一致。因此，从测试数据中来估计射手模型的参数，并不适合仅仅使用遗传算法。

2.2.2 无约束高斯- 牛顿优化算法

遗传算法迭代得到的模型参数估计值，在很大程度上接近于优化问题的全局最优解。为进一步得到精确的参数估计值，将遗传算法的解作为高斯- 牛顿优化算法的初始参数估计值，这种基于梯度的优化算法是解决最大似然优化问题的一种经典求解方法[22]。高斯- 牛顿优化算法的参数迭代更新方程由下式给出：

(i+1)=(i)+Δ,Δ

(16)

(17)

Mθ=θ(i)为Fisher信息矩阵，

(18)

(16)式中需要求解Fisher信息矩阵的逆矩阵，则要求Fisher信息矩阵可逆。通过求解Fisher信息矩阵可以得到Cramér-Rao的下界(CRLB)[23]，即参数估计可实现的最小方差，表示为Jc. 当未知模型参数过多时或数据中信息量不足时，信息矩阵M容易变成病态矩阵甚至奇异矩阵，导致迭代步长在某些方向上过大或者信息矩阵不可逆。因此，将信息矩阵M进行奇异值分解并求逆，有

(19)

式中：tks、uks分别为矩阵T和U的第ks列；δ1，…，δr为信息矩阵的r个奇异值。将奇异值从大到小排列，有

δ1≥δ2≥…≥δks≥δr.

(20)

一般而言，信息矩阵中最小的几个奇异值所包含的矩阵信息较少，且当奇异值过小时容易导致信息矩阵病态甚至不可逆，故可将最小的几个奇异值舍去。假设有η个奇异值满足舍去的标准：

(21)

式中：ε为计算机的计算精度；δmax为最大奇异值。经过上述处理后，信息矩阵的逆矩阵阶次从r阶降低到r-η阶，故该方法又称为降阶矩阵求逆法[24]。

此外，为防止高斯- 牛顿优化算法发散而终止迭代，将该算法切换成Nelder-Mead单纯形法。Nelder-Mead单纯形法是一种多维直接搜索的局部优化方法，在寻优过程中不必计算目标函数的梯度，只是针对一定图形的顶点，按照一定规则进行搜索，该方法操作简单且一直是收敛的，但是计算速度较慢，具体过程详见文献[25]。

(22)

3 仿真实验及结果分析

3.1 激励信号

下面介绍射手目标跟踪补偿实验，目标视场角信号驱动目标在视场内运动，射手根据视场误差角控制手柄，使光轴对准目标。为了防止射手预测信号的趋势，激励信号应表现一定的随机性，McRuer等[16]在实验中，将激励信号设计为由若干正弦信号叠加而成，其形式为

(23)

式中：Akf为正弦信号幅值；ωkf为正弦信号频率；φkf为正弦信号相位；kf=1,…,N，N为正弦信号个数。

为了得到较高精度的射手动力学模型，激励信号必须具有充分宽的频带和高的信噪比，激励信号的每个部分即Akf、ωkf和φkf都需要精确设计。因此将激励信号的频谱设计成1阶低通滤波器形状[12]，每个频率点的幅值为

(24)

表1 激励信号参数

3.2 实验装置和实验过程

实验在实验室开发的仿真系统上进行，该系统主要包括显示子系统、动力学仿真子系统和控制手柄，仿真系统原理图如图5所示。其中动力学仿真子系统在嵌入式系统RTX基础上开发，用以采集手柄的输出并实时仿真导引头系统动力学，计算频率为1 000 Hz；显示子系统采用vista2D+微软基础类库(MFC)方式开发，用以显示视场角误差，其显示的信息通过用户数据报协议(UDP)和动力学仿真子系统通信得到，然后射手根据显示信息操纵手柄，将导引头光轴对准目标以减小跟踪误差。

实验中共有5个参试人员，在采集数据前，所有参试人员需进行相当次数的训练，直到他们的跟踪表现稳定在一定水平。参试人员的跟踪水平可以通过误差信号均方根来评价，每次训练后参试人员都会被告知他们的跟踪评分，以便在初始熟悉期间激励改善他们的跟踪表现，且在达到一定熟练程度后用以保持稳定的表现。为了提高辨识结果精度，每个参试人员重复5次实验并采集相应的测量数据。

3.3 实验结果

实验中图1所示的所有信号都被记录下来，为了减少信号中噪声部分的影响并改善射手模型参数估计的精度，将每个参试人员的5次重复实验信号求平均值。图6所示为平均视场角误差信号Δq、射手输出信号u、目标视场角(激励信号)qt、光轴相对于惯性系夹角qs. 由图6可知参试人员成功将光轴对准了目标。

3.3.1 算法性能分析

表2 遗传算法参数域的上界和下界

遗传算法是一种基于概率启发式的搜索算法，在一定迭代步长限制下，每次迭代结果不完全相同。对单个参试人员的实验结果重复若干次遗传算法辨识运算，将迭代结果最小的10个解代入高斯- 牛顿优化算法进一步优化。图7(a)给出了这10个解的全部迭代过程，由图7(a)可知，重复10次遗传算法最终给出的值都不同，但通过高斯- 牛顿算法优化后最终都得到了相同的优化结果，即得到全局最小对数似然函数值，其对应的模型参数即为辨识结果。因为两个优化算法中残差方差的初始估计值都设置为1，所以得到的目标函数值较大，经过前述交替迭代过程最终得到各自的估计值，具体结果如图7(b)所示。

为了对比说明混合优化算法的优势，选择高斯- 牛顿优化算法结合单纯形法的复合优化算法[24]来估计模型参数。在表2参数域内随机选择100个初始参数向量，并分别使用高斯- 牛顿优化算法结合单纯形法的复合优化算法进行寻优，图8以升序方式给出了相应的结果。从图8可知，收敛至全局最小的参数向量小于一半，剩下的参数向量将导致收敛发散或者收敛至局部最小值。由此可见，采用混合优化算法能大大地减小初值选取工作，且能够精确地收敛至全局最小值。

3.3.2 结果分析

表3 模型参数辨识结果(参试人员1)

表4给出了本次实验所有参试人员的辨识结果，图10为相应的频率响应曲线。从图10可知，所有参试人员有着一致的辨识结果。其中不同参试人员的传递函数略有不同，这是因为每个参试人员有着不同的控制策略且个体之间的差异导致的，这种现象在该类型实验中普遍存在。

表4 所有参试人员的辨识结果

相关性系数ρu一般用于评价辨识模型的准确性，能够用来表示辨识模型的输出再现量测输出的能力。如果相关系数接近于1，则表示辨识模型的输出能很好地拟合量测输出，辨识模型可以较真实地反映实际对象的动态特性；反之，如果相关系数接近于0，则表示辨识结果较差，未能捕捉到实际对象的动态特性。相关系数ρu的定义[26]如下：

(25)

表4中参试人员1的相关系数为86%，可知模型输出能够较好地拟合量测输出，辨识模型能够反映射手动态特性。图11给出的参试人员1量测输出和模型输出的对比图也证明了该结论。表5给出了所有参试人员的相关系数。由表5可见，本实验得到的每个参试人员辨识模型相关性系数均大于70%，可以认为辨识模型能够较好地拟合射手的动力学特性。

相关系数参试人员12345ρu^u/%86.1872.6980.9585.7492.52

3.3.3 残差分析

射手输出除了线性输出部分，还包含模型为高斯噪声的残差，残差可以由辨识线性射手模型输出和实际量测射手的输出之差得到。图12(a)给出了参试人员1的残差时域图，其相应的功率谱密度(PSD)在图12(b)中给出。由图12可以看到，残差信号并不为白噪声，其特性近似为1个3阶低通滤波器：

(26)

式中：Kn为残差强度，Kn=0.26；ωn为残差滤波器的转折频率，ωn=0.8 rad/s；ξn为残差滤波器的阻尼系数，ξn=0.06.

射手模型残差的概率密度函数如图13所示。由图13可见：有残差信号的分布接近于零均值的高斯分布，同时也验证了之前的假设；此外残差信号谱特性并没有影响到信号的正态性。

4 结论

本文针对光纤制导武器中射手模型估计问题，采用基于输出误差框架下最大似然估计方法对射手模型参数和残差特性进行辨识。为提高优化算法的性能和鲁棒性，引入混合梯度优化算法，并成功应用于实验得到的数据中，经过多轮次迭代计算都能够找到全局最优解。辨识结果表明：算法能够从时域数据中精确地辨识得到射手模型参数，射手辨识模型能够真实地反映射手的动态过程；射手的残差信号基本符合正态分布，且具有色噪声的频率特性。