自适应迭代学习控制

2025-01-31

自适应迭代学习控制（精选7篇）

自适应迭代学习控制篇1

0 引言

迭代学习控制不需要依赖动态系统的精确数学模型就可以使系统实际输出完全跟踪期望输出。少量的先验知识就可以使迭代学习控制得以实现,并且重复运行次数越多控制精度越高,因此在很多具有重复运行特性的被控对象中得到了应用[1,2,3]。PID迭代学习控制以其良好的可靠性、鲁棒性已成为迭代学习控制领域中最为常用的控制策略。为了使PID学习律达到更好的控制效果,常将其与其他一些智能控制方法相结合,以期对非线性、不确定的复杂控制系统达到更好的控制效果。模糊控制对解决不确定性高的复杂非线性系统表现出诸多优越之处,因此很多学者将模糊控制和迭代学习控制相结合,以求达到较为理想的控制效果。张航等[4]采用自整定模糊控制器作为迭代学习律,并对单臂机械手进行了仿真研究,取得了较好的控制效果;张丽萍等[5]运用并行分配补偿方法确定T-S模型的迭代学习控制器结构,并给出了误差收敛条件;蒋思中等[6]利用作为模糊控制器的输入,生成PD增益矩阵的调整因子,以达到提高收敛速度的目的。

本文提出一种模糊自适应PID迭代学习控制算法,利用模糊控制器对PID学习律的参数进行实时整定,可以提高系统的控制精度,使控制系统具有较好的稳定性和较强的鲁棒性,同时可以提高系统的收敛速度。

1 问题的提出

考虑如下形式的线性系统:

其中t∈[0,T],x(t)∈Rm为系统的状态,u(t)∈Rn为系统的输入,y(t)∈Rr为系统的输出,A、B、C为适当维数的矩阵。

迭代学习控制是对给定控制系统在时间区间t∈[0,T]内,利用被控对象的期望输出yd(t),通过迭代学习算法寻找优化的控制信号uk(t),使得系统响应yk(t)相对于y0(t)有所改善,寻找uk(t)的过程就是系统学习的过程,且使得当k趋于∞时,满足yk(t)与yd(t)的偏差为零,即系统实现实际输出完全跟踪期望输出[7]。

在第次运行时,系统的动态方程可表示为:

闭环PID型学习律如下:

式中Γ、Φ、Ψ为学习增益矩阵,ek+1(t)为实际输出与期望输出的偏差:

可见PID学习律中的参数是一成不变的,目前所需要的就是寻找到一种方法来对式中Γ、Φ、Ψ三个参数进行实时校正,实现系统的动态学习过程,以提高系统的控制精度和收敛速度。

2 模糊自适应PID迭代学习控制算法

2.1 控制系统基本结构

模糊自适应PID迭代学习控制算法结构如图1所示,控制系统采用闭环学习控制算法,利用传统的经验PID参数作为参考,通过模糊整定单元对经验PID参数进行实时校正,生成精确度更高的模糊PID学习律,提高系统收敛速度,加强抗干扰能力,最终实现系统的完全跟踪,即。

2.2 模糊自适应PID学习律

模糊控制器的输入选用系统输出误差e和输出误差的变化量ec,而模糊控制器的输出则为PID学习律的增益整定量,采用PID学习律时,Kp用于提高系统响应速度,调节系统控制精度,过大则导致系统不稳定,过小则导致系统响应速度缓慢;Ki用于消除系统稳态误差,过大则会导致积分饱和,超调量较大,过小则无法起到相应作用;Kd有助于提高系统的动态性能,模糊控制器生成的整定参数应以此为原则而确定,其控制规则形式如下:

IF ek(t)is E and dek(t)is EC THEN Kpk(t)is Kp。

IF ek(t)is E and dek(t)is EC THEN Kik(t)is Ki。

IF ek(t)is E and dek(t)is EC THEN Kdk(t)is Kd。

定义模糊控制器的输入、输出量的模糊子集均为NB,NM,NS,ZE,PS,PM,PB,即负大,负中,负小,零,正小,正中,正大,其模糊论域均为[-E,E],同时可以建立如表1所示控制规则表:

通过模糊规则对其进行判定,并利用重心法对其进行解模糊处理,以Kp为例,可得模糊控制器输出的的PID参数分别为:

式中μKp为输出量的隶属度函数,Kpq为输出的模糊化变量,Kp为解模糊后的输出量。

由此可以得到模糊PID迭代学习律为:

其中Γ、Φ、Ψ为PID学习律的固定参数。

3 仿真研究

利用仿真对算法的有效性进行验证,选取如下形式LTI系统作为被控对象:

期望跟踪轨迹选用:。各模糊控制器的输入输出论域均选择在[-3,3]范围内,经验PID学习律参数为:Γ=30,Φ=2,Ψ=3。学习周期为T=6s。各变量均采用如图2形式的隶属度函数,其中ZE选用三角形隶属度函数,其余均选择梯形隶属度函数。

控制系统仿真框图如图3所示:

通过仿真可以得到如图4所示的第一次迭代学习下的跟踪轨迹,图5则为系统跟踪误差。

从仿真结果可以看出,系统具有较好的跟踪性能,并随学习次数的增加系统误差不断趋近于零,使得控制系统具有良好的动态性能,达到高精度趋近期望轨迹的控制目的,由此可见模糊自适应PID迭代学习控制系统对于轨迹跟踪具有较好的控制效果。

4 结语

现实问题中常常会遇到需要高精度轨迹跟踪的控制问题,针对这一问题本文提出一种模糊自适应PID迭代学习控制算法,利用模糊推理整定PID学习律的参数,以达到精确控制和快速收敛的目的,并通过仿真验证了方法的有效性,对同类控制对象的控制系统设计具有一定的参考价值。

参考文献

[1]王仲鸿,王强,张东霞,姜齐荣.电力系统暂态稳定问题和迭代学习控制的研究[J].电力系统自动化,1999,23(8):6～10.

[2]徐敏,林辉.基于迭代学习控制理论的励磁控制器设计[J].电力自动化设备,2006,26(3):69～72.

[3]王岩,付永领.模糊滑模迭代学习控制算法在液压系统中应用[J].北京航空航天大学学报,2007,33(1):86～89.

[4]张航,罗大庸,黄浩江,罗熊.机器人模糊迭代学习控制及其仿真研究[J].自动化技术与应用,2002,(2):3～5.

[5]张丽萍,杨富文.基于T-S模型的迭代学习控制算法及其在机器人点位控制中的仿真研究[J].系统仿真学报,2005,17(1):166～169.

[6]蒋思中,朱芳来,王心开,王改云.模糊增益PD迭代学习算法及其应用[J].电光与控制,2009,16(8):72～74.

[7]李士勇.模糊控制·神经控制和智能控制论(第2版)[M].哈尔滨:哈尔滨工业大学出版社,1998:495～496.

自适应迭代学习控制篇2

目前针对装载机的轨迹跟踪控制问题已经取得一定程度的进展。柳波,等基于模糊自整定PID控制实现了装载机自动铲掘时的轨迹控制问题[1]。龚捷,等在考虑作业阻力补偿的情况下,基于计算力矩法设计装载机铲掘作业的轨迹跟踪控制器,实现铲斗轨迹跟踪控制[2]。

自适应迭代学习控制整合自适应控制和迭代学习控制的优点,基于PD反馈控制器增益自适应调节,迭代学习控制器修正控制力矩,该算法具有快速的渐进收敛性能[3]。拉格朗日法与牛顿-欧拉法是机构动力学建模的两种常见方法,封闭形式的拉格朗日方程是控制器设计和综合的有效工具,而递推形式的牛顿-欧拉方程则是实时控制、求解逆动力学问题的有效方法[4]。采用牛顿-欧拉法推导了装载机工作装置的二自由度动力学模型,在matlab环境下设计自适应迭代学习控制器并进行稳定性证明,对装载机工作装置的轨迹跟踪问题进行仿真控制,以期达到高精度追踪。

1 基于牛顿-欧拉法的动力学方程

将装载机工作装置机器人化处理,视为平面二连杆操作臂机构,建立如图1所示的D-H坐标系。设定动臂与车架铰接点为O0,铲斗与动臂铰接点为O1,铲斗斗尖位置为O2,分别以O0、O1、O2为原点建立坐标系O0x0y0z0、O1x1y1z1、O2x2y2z2。

图1中,θ1、θ2分别为动臂和铲斗相对于车架和动臂的转角;l1=lO0O1、l2=lO1O2分别为动臂和铲斗长度;m1、m2分别为铲斗和动臂的质量;假定动臂和铲斗的质心处于连杆O0O1、O1O2的中心位置。

1.1 牛顿-欧拉动力学方程

根据旋转关节的运动推导关节力矩的过程可细分为外推法和内推法两个步骤。外推法采用牛顿欧拉方程,按连杆编号从小到大的顺序逐次计算连杆的速度和加速度;内推法则按连杆编号从大到小的顺序计算连杆间的力、力矩[5]。运动学方程如下。

外推法求解过程:i=0,1,…,n。

内推法求解过程:i=n,n-1,…,1。

式中:iωi、 $^{i} ω_{\dot{i}}$ 、ivi、 $^{i} b_{\dot{i}}$ 为连杆i的角速度、角加速度、线速度、线加速度相对于坐标系{i}的表达式;i+1iR和 $_{i + 1}^{i}$ R为坐标系{i}和{i+1}之间的姿态变换矩阵,且i+1iR $_{i + 1}^{i}$ R=I;i $\hat{Ζ}$ i为Zi轴方向上的单位矢量;iPi+1为节点i到节点i+1的矢量;iPCi为节点i到连杆i质心位置处的矢量;CiIi为连杆i质心处的惯性张量;mi为连杆i的质量;为连杆i质心上的惯性力和力矩相对于坐标系{i}的表达式;为连杆i-1作用在连杆i上的力和力矩相对于坐标系{i}的表达式;τi为关节i上的旋转力矩。

1.2 装载机工作装置的动力学计算

以图1中装载机工作装置的二连杆机器人操作臂为研究对象,推导其牛顿-欧拉动力学方程。假设装载机铲斗斗尖处不受力的作用f3=0,n3=0;车架固定不旋转有 $ω_{0} = 0, \overset{\cdot}{ω_{0}} = 0$ ;计及重力因素[6]的影响,在式(1)中令 $v_{\overset{\cdot}{i + 1}}_{0} = G - G y_{0} ‚ g = 9.8 m / s^{- 2}$ 。相邻坐标系之间的姿态变换矩阵为(其中ci+1=cos(θi+1)),(si+1=sin(θi+1))。

$\begin{array}{l} _{i + 1}^{i} R = [\begin{matrix} c_{i + 1} & - s_{i + 1} & 0 \\ s_{i + 1} & c_{i + 1} & 0 \\ 0 & 0 & 1 \end{matrix}], \\ _{i}^{i + 1} R = [\begin{matrix} c_{i + 1} & s_{i + 1} & 0 \\ - s_{i + 1} & c_{i + 1} & 0 \\ 0 & 0 & 1 \end{matrix}], i = 0, 1 (3) \end{array}$

连杆在坐标系{i}中的惯性张量为:

根据平行移轴定理[7]得到以刚体质心{C}为原点的坐标系与{i}坐标系间的转换式

${\begin{cases} ^{i} Ι_{x x} =^{C_{i}} Ι_{x x} + m (y_{c}^{2} + z_{c}^{2}) \\ ^{i} Ι_{y y} =^{C_{i}} Ι_{y y} + m (x_{c}^{2} + z_{c}^{2}),^{i} Ι_{z z} = \\ ^{C_{i}} Ι_{z z} + m (x_{c}^{2} + y_{c}^{2}) \\ ^{i} Ι_{x y} =^{C_{i}} Ι_{x y} - m x_{c} y_{c},^{i} Ι_{x z} = \\ ^{C_{i}} Ι_{x z} - m x_{c} z_{c},^{i} Ι_{y z} =^{C_{i}} Ι_{y z} - m y_{c} z_{c} \end{cases} (6)$

式(6)中pc=[xcyczc]T表示连杆质心在坐标系{i}中的位置矢量。

将式(3)—式(6)代入式(1)和式(2)计算装载机工作装置的动力学方程如下(式子太长未列出)。

$\begin{array}{l} ^{C_{1}} Ι_{1} = [\begin{matrix} 0 & 0 & 0 \\ 0 & \frac{m_{1} l_{1}^{2}}{12} & 0 \\ 0 & 0 & \frac{m_{1} l_{1}^{2}}{12} \end{matrix}], \\ ^{C_{2}} Ι_{2} = [\begin{matrix} 0 & 0 & 0 \\ 0 & \frac{m_{2} l_{2}^{2}}{12} & 0 \\ 0 & 0 & \frac{m_{2} l_{2}^{2}}{12} \end{matrix}] (7) \end{array}$

设 $\vec{i}$ 、 $\vec{j}$ 、 $\vec{k}$ 分别为坐标系O0-x0y0z0中的单位向量

${\begin{cases} ^{1} ω_{1} = \dot{θ}_{1} \vec{k} \\ ^{1} \dot{ω}_{1} = \ddot{θ}_{1} \vec{k} \\ ^{1} \dot{v}_{1} = g s_{1} \vec{i} + g c_{1} \vec{j} \\ ^{1} \dot{v}_{C_{1}} = (- \frac{l_{1}}{2} \ddot{θ}_{1} + g s_{1}) \vec{i} + (\frac{l_{1}}{2} \ddot{θ}_{1} + g c_{1}) \vec{j} \\ ^{2} ω_{2} = (\dot{θ}_{1} + \dot{θ}_{2}) \vec{k} \\ ^{1} \dot{ω}_{2} = (\ddot{θ}_{1} + \ddot{θ}_{2}) \vec{k} \\ ^{2} \dot{v}_{2} = (l_{1} \ddot{θ}_{1} s_{2} - l_{1} \dot{θ}_{1}^{2} c_{2} + g s_{12}) \vec{i} + (l_{1} \ddot{θ}_{1} c_{2} + \\ l_{1} \dot{θ}_{1}^{2} s_{2} + g c_{12}) \vec{j} \\ ^{2} \dot{v} C_{2} = (l_{1} \ddot{θ}_{1} s_{2} - l_{1} \dot{θ}_{1}^{2} c_{2} + g s_{12} - \frac{l_{2}}{2} (\dot{θ}_{1} + \dot{θ}_{2})^{2}) \vec{i} + \\ (l_{1} \ddot{θ}_{1} c_{2} + l_{1} \dot{θ}_{1}^{2} s_{2} + g c_{12} + \frac{l_{2}}{2} (\ddot{θ}_{1} + \ddot{θ}_{2})) \vec{j} \\ ^{1} Ν_{1} = \frac{m_{1} l_{1}^{2}}{12} \ddot{θ}_{1} \vec{k} \\ ^{1} F_{1} = (- \frac{m_{1} l_{1}}{2} \ddot{θ}_{1} + m_{1} g s_{1}) \vec{i} + (\frac{m_{1} l_{1}}{2} \ddot{θ}_{1} + m_{1} g c_{1}) \vec{j} \\ ^{2} Ν_{2} = \frac{m_{2} l_{2}^{2}}{12} (\ddot{θ}_{1} + \ddot{θ}_{2}) \vec{k} \\ ^{2} F_{2} = m_{2} (l_{1} \ddot{θ}_{1} s_{2} - l_{1} \dot{θ}_{1}^{2} c_{2} + g s_{12} - \frac{l_{2}}{2} (\dot{θ}_{1} + \dot{θ}_{2})^{2}) \vec{i} + \\ m_{2} (l_{1} \ddot{θ}_{1} c_{2} + l_{1} \dot{θ}_{1}^{2} s_{2} + g c_{12} + \frac{l_{2}}{2} (\ddot{θ}_{1} + \ddot{θ}_{2})) \vec{j} \end{cases} (8)$

${\begin{cases} τ_{1} = (\frac{m_{1} l_{1}^{2}}{3} + \frac{m_{2} l_{2}^{2}}{3} + m_{2} l_{1}^{2} + m_{2} l_{1} l_{2} c_{2}) \ddot{θ_{1}} + \\ (\frac{m_{2} l_{2}^{2}}{3} + \frac{m_{2} l_{1} l_{2} c_{2}}{2}) θ_{2}^{\ddot{2}} - \frac{m_{2} l_{1} l_{2} s_{2}}{2} \overset{\cdot}{θ_{2}}^{2} - \\ m_{2} l_{1} l_{2} s_{2} \overset{\cdot}{θ_{1}} \overset{\cdot}{θ_{2}} + \frac{m_{2} g l_{2}}{2} c_{12} + (\frac{m_{1}}{2} + m_{2}) g l_{1} c_{1} \\ τ_{2} = (\frac{m_{2} l_{2}^{2}}{3} + \frac{m_{2} l_{1} l_{2} c_{2}}{2}) \ddot{θ_{1}} + \frac{m_{2} l_{2}^{2}}{3} \ddot{θ_{2}} + \\ \frac{m_{2} l_{2} s_{2}}{2} \overset{\cdot}{θ_{1}}^{2} + \frac{m_{2} l_{2} g c_{12}}{2} \end{cases} (9)$

2 工作装置轨迹的自适应迭代学习控制

2.1 问题描述

考虑实际中存在不确定性、摩擦及外部干扰的影响,式(9)描述的装载机工作装置动力学方程可写为

$Μ (q_{i} (t)) \ddot{q_{i}} (t) + C (q_{i} (t), \ddot{q_{i}} (t)) \ddot{q_{i}} (t) + G (q_{i} (t)) = τ_{i} (t) + d_{i} (10)$

式(10)中:qi(t)=θi(t),i=1,2,t∈[0,T]为时间变量,i∈Z+为迭代次数, $q_{i} \in R^{n}, \overset{\cdot}{q_{i}} \in R^{n}, \ddot{q_{i}} \in R^{n}$ 表示第i次迭代的关节角度、角速度、角加速度量。M(qi)∈Rn×n为惯性矩阵, $C (q_{i}, \overset{\cdot}{q_{i}}) \overset{\cdot}{q_{i}} \in R^{n}$ 表示科氏力和离心力量,G(qi)∈Rn为重力项,τi(t)∈Rn为作用在关节位置处的控制力或力矩,di∈Rn为未建模动态和扰动。

2.2 收敛性分析

控制目的是∀t∈[0,T]及∀i∈Z+,设计有界的自适应迭代学习控制器τi(t)来保证 $q_{i} \in R^{n}, \overset{\cdot}{q_{i}} \in R^{n}, \overset{\cdot \cdot}{q_{i}} \in R^{n}$ 有界,且当i→∞时, $q_{i} (t), \overset{\cdot}{q_{i}} (t)$ 收敛到理想的参考位置 $q_{d} (t), \overset{\cdot}{q_{d}} (t)$ 。为实现控制目的,给出以下基本假设:

(A1) qd(t)为可实现的关节参考轨迹角度。

(A2) 对于∀t∈[0,T]及 $\forall i \in Ζ_{+} ‚ q_{i} (t), \overset{\cdot}{q_{i}} (t), \overset{\cdot \cdot}{q_{i}} (t) ‚ d_{i} (t)$ 有界。

(A3) ∀i∈Z+,初始条件满足:

$\overset{\cdot}{q_{d}} (0) - \overset{\cdot}{q_{i}} (0) = q_{d} (0) - q_{i} (0) = 0$ 。

且满足如下4个基本特征:

(B1) M(qi)∈Rn×n为有界正定对称矩阵;

$(B_{2}) \dot{Μ} (q_{i}) - 2 C (q_{i}, \overset{\cdot}{q_{i}}) \in R^{n \times n}$ 为对称矩阵,且 $ω^{Τ} (\dot{Μ} (q_{i}) - 2 C (q_{i}, \overset{\cdot}{q_{i}})) ω = 0, \forall ω \in R^{n}$ ;

$(B_{3}) G (q_{i}) + C (q_{i}, \overset{\cdot}{q_{i}}) \overset{\cdot}{q_{d}} (t) = ψ (q_{i}, \overset{\cdot}{q_{i}}) ξ^{Τ} (t), ψ (q_{i}, \overset{\cdot}{q_{i}}) \in R^{n \times (m - 1)}$ 为已知矩阵,ξ(t)∈Rm-1为未知向量;

$(B_{4}) \forall q_{i}, \overset{\cdot}{q_{i}} ‚ t \in [0, Τ] ‚ i \in Ζ_{+} ‚ | | C (q_{i}, \overset{\cdot}{q_{i}}) | | \leq k_{c} | | \overset{\cdot}{q_{i}} | | ‚ | | G (q_{i}) | | \leq k_{g}, k_{c}$ 、kg均为正实数。

定理定义关节位置跟踪误差和速度误差分别为: $e_{i} (t) = q_{d} (t) - q_{i} (t), \overset{\cdot}{e_{i}} (t) = \dot{q}_{d} (t) - \dot{q}_{i} (t)$ ,采用控制律

${\begin{cases} τ_{i} (t) = Κ_{Ρ} e_{i} (t) + Κ_{D} \overset{\cdot}{e_{i}} (t) + φ (q_{i}, \overset{\cdot}{q_{i}}, \overset{\cdot}{e_{i}} (t)) \overset{⌢}{θ_{i}} (t) (11 a) \\ \overset{⌢}{θ_{i}} (t) = \overset{⌢}{θ}_{i - 1} (t) + Γ φ^{Τ} (q_{i}, \overset{\cdot}{q_{i}}, \overset{\cdot}{e_{i}} (t)) \overset{\cdot}{e_{i}} (t) (11 b) \end{cases}$

式(11)中: $\overset{⌢}{θ}_{i - 1} (0) = 0 ‚ φ (q_{i}, \overset{\cdot}{q_{i}}, \overset{\cdot}{e_{i}}) \in R^{n \times n}$ ,且 $φ (q_{i}, \overset{\cdot}{q_{i}}, \overset{\cdot}{e_{i}}) = [ψ (q_{i}, \dot{q}_{i}) sgn (\overset{\cdot}{e_{i}})]$ ,矩阵KP、KD、Γ均为对称正定矩阵,则ei(t)、 $\overset{\cdot}{e_{i}} (t)$ 有界,且

$\lim_{i \to \infty} e_{i} (t) = \lim_{i \to \infty} \overset{\cdot}{e_{i}} (t) = 0, \forall t \in [0, Τ]$ 。

证明第i次迭代时,构造Lyapunov函数

$W_{i} (t) = V_{i} (t) + \frac{1}{2} \int_{0}^{t} \overset{⌢}{θ_{i}} (t) Γ^{- 1} \overset{⌢}{θ_{i}} (t) d τ (12)$

式(12)中,设定 $\overset{⌢}{θ}_{i} (t) = θ (t) - \overset{⌢}{θ}_{i} (t) ‚ θ (t) = [ξ^{Τ} (t) β]^{Τ}, \hat{θ}_{i} (t) = [\hat{ξ}_{i}^{Τ} (t) \hat{β}_{i} (t)]^{Τ}$ 为θ(t)的估计值,由假设(A1)、和基本特征(B1)知,

$| | Μ (q_{i}) \overset{\cdot \cdot}{q}_{d} - d_{i} | | \leq β ‚ β > 0$

$\begin{array}{l} V_{i} (e_{\dot{i}} (t), e_{i} (t)) = \frac{1}{2} e_{i}^{\dot{Τ}} (t) Μ (q_{i}) e_{\dot{i}} (t) + \\ \frac{1}{2} e_{i} (t)^{Τ} Κ_{p} e_{i} (t) (13) \end{array}$

由于 $\bar{θ_{i}} = - θ + \overset{⌢}{θ}_{i} - \overset{⌢}{θ}_{i + 1} + θ = - \overset{⌢}{θ}_{i} + \overset{⌢}{θ}_{i + 1}$ ,

即 $\overset{⌢}{θ}_{i + 1} = \bar{θ}_{i} + \overset{⌢}{θ}_{i}$ ,则

$\begin{array}{l} \tilde{θ}_{i} (t) Γ^{- 1} \tilde{θ}_{i} (t) - \tilde{θ}_{i - 1} (t) Γ^{- 1} \tilde{θ}_{i - 1} (t) = \\ - 2 \bar{θ}_{i}^{Τ} (t) Γ^{- 1} \tilde{θ}_{i} (t) - \bar{θ}_{i}^{Τ} (t) Γ^{- 1} \bar{θ}_{i} (t) (14) \end{array}$

2.2.1 Wi(t)的非递增性证明

$Δ W_{i} = W_{i} - W_{i - 1} = V_{i} - V_{i - 1} - \frac{1}{2} \int_{0}^{t} (\bar{θ}_{i}^{Τ} Γ^{- 1} \bar{θ}_{i} +$

$2 \bar{θ}_{i - 1}^{Τ} (t) Γ^{- 1} \bar{θ}_{i - 1} (t)) d τ$ (15)

式(15)中 $\bar{θ}_{i} = \overset{⌢}{θ}_{i} - \overset{⌢}{θ}_{i - 1}$ 。由于 $\int_{0}^{t} \dot{V} (t) d τ = V_{i} (t) - V_{i} (0)$ ,即 $V_{i} (t) = V_{i} (0) + \int_{0}^{t} \dot{V} (t) d τ$ ,

又由于

${\begin{cases} V_{\dot{i}} (t) = e_{i}^{\dot{Τ}} (t) Μ (q_{i}) e_{\overset{\cdot \cdot}{i}} (t) + \frac{1}{2} e_{i}^{\dot{Τ}} (t) \times \\ \dot{Μ} (q_{i}) e_{\dot{i}} (t) + e_{i}^{\dot{Τ}} (t) Κ_{p} e_{i} \\ V_{i} (e_{\dot{i}} (t), e_{i} (t)) = \\ V_{i} (e_{\dot{i}} (0), e_{i} (0)) + \\ \int_{0}^{t} [e_{i}^{\dot{Τ}} Μ (q_{i}) e_{\overset{\cdot \cdot}{i}} + \frac{1}{2} e_{i}^{\dot{Τ}} \dot{Μ} (q_{i}) e_{\dot{i}} + e_{i}^{\dot{Τ}} Κ_{p} e_{\dot{i}}] d τ \end{cases} (16)$

由式(10)和特性(B2)和特性(B3)得

${\begin{cases} e_{i}^{\dot{Τ}} Μ e_{\overset{\cdot \cdot}{i}} = e_{i}^{\dot{Τ}} Μ (q_{\overset{\cdot \cdot}{d}} - q_{\overset{\cdot \cdot}{i}}) = e_{i}^{\dot{Τ}} Μ q_{\overset{\cdot \cdot}{d}} - \\ e_{i}^{\dot{Τ}} (- C q_{\dot{i}} - G + τ_{i} + d_{i}) \\ \frac{1}{2} e_{i}^{\dot{Τ}} \dot{Μ} (q_{i}) e_{\dot{i}} = e_{i}^{\dot{Τ}} C e_{\dot{i}} = e_{i}^{\dot{Τ}} C (q_{\dot{d}} - q_{\dot{i}}) = \\ e_{i}^{\dot{Τ}} C q_{\dot{d}} - e_{i}^{\dot{Τ}} C q_{\dot{i}} \end{cases} (17)$

由已知得

$\begin{array}{l} e_{i}^{\dot{Τ}} (Μ (q_{i}) q_{\overset{\cdot \cdot}{d}} - d_{i}) \leq \\ β | | e_{i}^{\dot{Τ}} | | = e_{i}^{\dot{Τ}} β sgn (e_{\dot{i}}) ψ (q_{i}, q_{\dot{i}}) ξ^{Τ} + β sgn (e_{\dot{i}}) = \\ [ψ (q_{i}, q_{\dot{i}}) sgn (e_{\dot{i}})] [ξ^{Τ} β]^{Τ} = ϕ (q_{i}, q_{\dot{i}}, e_{\dot{i}}) θ (18) \end{array}$

则

$\begin{array}{l} V_{i} (e_{\dot{i}} (t), e_{i} (t)) = \\ V_{i} (e_{\dot{i}} (0), e_{i} (0)) + \int_{0}^{t} e_{i}^{\dot{Τ}} (Μ (q_{i}) q_{\overset{\cdot \cdot}{d}} - d_{i} + \\ C (q_{i}, q_{\dot{i}}) q_{\dot{d}} + G (q_{i}) + Κ_{Ρ} e_{i} - τ_{i}) d τ \leq V_{i} (e_{\dot{i}} (0), e_{i} (0)) + \int_{0}^{t} e_{i}^{\dot{Τ}} (ψ (q_{i}, q_{\dot{i}}) ξ^{Τ} + Κ_{Ρ} e_{i} + β sgn (e_{\dot{i}}) - \\ τ_{i}) d τ \leq V_{i} (e_{\dot{i}} (0), e_{i} (0)) + \\ \int_{0}^{t} e_{i}^{\dot{Τ}} (ϕ (q_{i}, q_{\dot{i}}, e_{i}) θ + Κ_{Ρ} e_{i} - τ_{i}) d τ (19) \end{array}$

将控制律式(11a)代入式(19),得

$\begin{array}{l} V_{i} (e_{\dot{i}} (t), e_{i} (t)) \leq V_{i} (e_{\dot{i}} (0), e_{i} (0)) + \\ \int_{0}^{t} e_{i}^{\dot{Τ}} (ϕ (q_{i}, q_{\dot{i}}, e_{i}) \tilde{θ}_{i} - Κ_{D} e_{\dot{i}}) d τ (20) \end{array}$

据式(11b)得: $\bar{θ}_{i}^{Τ} (t) = (Γ ϕ^{Τ} e_{\dot{i}})^{Τ} = e_{i}^{\dot{Τ}} ϕ Γ$ ,有

由假设(A2)得 $V_{i} (\overset{\cdot}{e_{i}} (0), e_{i} (0)) = 0$ ,将式(21)代入式(15),有

$\begin{array}{l} Δ W_{i} = - V_{i - 1} + V_{i} - \frac{1}{2} \int_{0}^{t} (\bar{θ}_{i}^{Τ} Γ^{- 1} \bar{θ}_{i} + 2 \bar{θ}_{i}^{Τ} Γ^{- 1} \tilde{θ}_{i}) d τ \leq \\ - V_{i - 1} + \int_{0}^{t} e_{i}^{\dot{Τ}} (ϕ \tilde{θ}_{i} - Κ_{D} e_{\dot{i}}) d τ - \frac{1}{2} \int_{0}^{t} (e_{i}^{\dot{Τ}} ϕ Γ ϕ^{Τ} e_{\dot{i}} + 2 e_{i}^{\dot{Τ}} ϕ \tilde{θ}_{i}) d τ \leq - V_{i - 1} - \frac{1}{2} \int_{0}^{t} e_{i}^{\dot{Τ}} (ϕ Γ ϕ^{Τ} + \\ 2 Κ_{D}) e_{\dot{i}} d τ \leq 0 (22) \end{array}$

因Vi-1、Γ、KD均为正定阵,ΔWi≤0,因此Wi为非递增数列,得结论如下:若W0有界,那么Wi必定有界。

2.2.2 W0(t)的连续有界性证明

由式(12)和式(20)得

$\begin{array}{l} W_{\dot{0}} (t) \leq e_{0}^{\dot{Τ}} (ϕ (q_{0}, q_{\dot{0}}, q_{\overset{\cdot \cdot}{0}}) \tilde{θ}_{0} - Κ_{D} e_{\dot{0}}) + \\ \frac{1}{2} \tilde{θ}_{0}^{Τ} Γ^{- 1} \tilde{θ}_{0} (23) \end{array}$

由于 $\overset{⌢}{θ}_{- 1} (t) = 0$ ,根据控制律式(11b)得

$\overset{⌢}{θ}_{0} (t) = Γ φ^{Τ} (q_{0}, \overset{\cdot}{q_{0}}, \overset{\cdot}{e_{0}}) \overset{\cdot}{e_{0}} (t)$ ,有

$\begin{array}{l} W_{\dot{0}} (t) \leq - e_{0}^{\dot{Τ}} Κ_{D} e_{\dot{0}} + (\hat{θ}_{0}^{Τ} + \frac{1}{2} \tilde{θ}_{0}^{Τ}) Γ^{- 1} \tilde{θ}_{0} = \\ - e_{0}^{\dot{Τ}} Κ_{D} e_{\dot{0}} - \frac{1}{2} \tilde{θ}_{0}^{Τ} Γ^{- 1} \tilde{θ}_{0} + θ^{Τ} Γ^{- 1} \tilde{θ}_{0} (24) \end{array}$

因为 $θ^{Τ} Γ^{- 1} \overset{⌢}{θ}_{0} \leq Κ | | Γ^{- 1} \overset{⌢}{θ}_{0} | |^{2} + \frac{1}{4 Κ} | | θ | |^{2}, Κ > 0$ ,可得

$\overset{\cdot}{W_{0}} (t) \leq - β_{1} | | \overset{\cdot}{e_{0}} | |^{2} - β_{2} | | \overset{⌢}{θ}_{0} | |^{2} + \frac{1}{4 Κ} | | θ | |^{2} (25)$

$\begin{array}{l} 式 (25) 中 ∶ β_{1} = λ_{\min} (Κ_{D}), \\ β_{2} = \frac{1}{2} λ_{\min} (Γ^{- 1}) - Κ λ_{\max}^{2} (Γ^{- 1}), Κ \leq \frac{λ_{\min} (Γ^{- 1})}{2 λ_{\max}^{2} (Γ^{- 1})} ‚ 因此有 \end{array}$

$\overset{\cdot}{W_{0}} (t) \leq \frac{1}{4 Κ} | | θ | |^{2} ‚ \forall t \in [0, Τ] (26)$

因为θ(t)连续有界,推出W0(t)亦连续有界。

2.2.3 Wi(t)的连续有界性证明

Wi(t)可表示为 $W_{i} (t) = W_{0} + \sum_{j = 1}^{i} Δ W_{j}$ 。

则由式(15)可得

$\begin{array}{l} W_{i} \leq W_{0} - \sum_{j = 1}^{i} V_{j - 1} \leq W_{0} - \sum_{j = 1}^{i} V_{j - 1} \leq W_{0} - \\ \frac{1}{2} \sum_{j = 1}^{i} e_{j - 1}^{Τ} Κ_{Ρ} e_{j - 1} - \frac{1}{2} \sum_{j = 1}^{i} e_{j - 1}^{\dot{Τ}} Κ_{Ρ} e_{\overset{\cdot}{j - 1}} (27) \end{array}$

$\begin{array}{l} 那么 (\sum_{j = 1}^{i} e_{j - 1}^{Τ} Κ_{Ρ} e_{j - 1} + \sum_{j = 1}^{i} e_{j - 1}^{\dot{Τ}} Κ_{Ρ} e_{\overset{\cdot}{j - 1}}) \leq \\ 2 (W_{0} - W_{i}) \leq 2 W_{0}, 即 W_{i} (t) 有界。可得结论 ∶ \end{array}$

$\lim_{i \to \infty} e_{i} (t) = \lim_{i \to \infty} \overset{\cdot}{e_{i}} (t) = 0, \forall t \in [0, Τ]$ ,定理得证。

3 仿真研究

对于式(10)所示的装载机工作装置的动力学方程及式(11)所设计的控制器,仿真研究如下。设置动臂、铲斗质量分别为m1=1 793.080 7 kg、m2=2 706.447 8 kg;动臂、铲斗长度l1=2.668 2 m、l2=1.262 1 m。对式(10)描述的系统,可写为

$[\begin{matrix} m_{11} & m_{12} \\ m_{21} & m_{22} \end{matrix}] [\begin{array}{l} \overset{\cdot \cdot}{θ_{1}} \\ \overset{\cdot \cdot}{θ_{2}} \end{array}] + [\begin{matrix} c_{11} & c_{12} \\ c_{21} & c_{22} \end{matrix}] [\begin{array}{l} \overset{\cdot}{θ_{1}} \\ \overset{\cdot}{θ_{2}} \end{array}] + [\begin{matrix} G_{1} \\ G_{2} \end{matrix}] = [\begin{matrix} τ_{1} \\ τ_{2} \end{matrix}] (28)$

式(28)中, $m_{11} = \frac{m_{1} l_{1}^{2}}{3} + \frac{m_{2} l_{2}^{2}}{3} + m_{2} l_{1}^{2} + m_{2} l_{1} l_{2} c_{2}$ ; $m_{12} = \frac{m_{2} l_{2}^{2}}{3} + \frac{m_{2} l_{1} l_{2} c_{2}}{2}$ ; $m_{21} = \frac{m_{2} l_{2}^{2}}{3} + \frac{m_{2} l_{1} l_{2} c_{2}}{2};$

在matlab环境下编制m函数[8],创建装载机工作装置的自适应迭代学习simulink主程序框图如图2所示。其中,ctrl为控制器子程序,adapt为自适应律子程序,plant为被控对象子程序,input为指令输入子程序。

设定式(11)的PD控制器参数KP=KD=diag[100,100],自适应律参数Γ=diag[150,150,150,150,150],迭代次数为5,仿真时间设定为1 s。干扰项di(t)=[dmsint dmsint],dm为幅值为1的随机信号。对装载机工作装置动臂、铲斗分别施加q1d=sin(2πt),q2d=cos(2πt)的关节位置指令信号,取被控对象初始状态为x(0)=[0 2π 1 0]T,得到动臂、铲斗的跟踪控制仿真结果如图3—图6所示。

图3显示5次迭代时动臂和铲斗的关节位置追踪情况,期望值与迭代值之间存在误差,且随着迭代次数增加,误差越来越小;图4显示5次迭代后动臂和铲斗有极好的速度追踪效果;

图5、图6分别显示了动臂和铲斗的关节位置跟踪误差和速度跟踪误差图,位置和速度跟踪误差随迭代次数逐渐减小。

4 结论

在运用牛顿-欧拉法对装载机工作装置进行动力学分析的基础上,针对其不确定性和外部干扰性, 采用自适应迭代学习控制策略对其进行轨迹跟踪仿真实验。仿真结果表明,该控制器能保证系统在一定时间内稳定地减小跟踪误差,追踪性能良好,充分验证了该算法的有效性和可行性。

摘要：在理想情况下用牛顿-欧拉法建立装载机工作装置的二自由度动力学模型,考虑实际情况中存在的建模误差和扰动干扰,设计自适应迭代学习控制器对其进行轨迹跟踪控制研究,基于Lyapunov函数证明了跟踪误差的稳定性和收敛性,在matlab环境下对所设计的控制器进行仿真研究。仿真结果表明,所设计的控制器对装载机工作装置有极好追踪效果,验证了该算法的可行性与有效性。

关键词：装载机,工作装置,牛顿-欧拉方程,自适应迭代学习控制,轨迹跟踪

参考文献

[1]柳波,黄杰,等.滑移装载机自动铲掘轨迹控制研究.武汉理工大学学报,2010;32(23):107—115

[2]龚捷,鲍金峰,等.基于计算力矩法的装载机工作装置轨迹控制.机械工程学报,2010;46(13):141—146

[3]刘松国.六自由度串联机器人运动优化与轨迹跟踪控制研究.浙江大学学位论文,2009

[4]刘美华,常文森,等.机械手的显式递推牛顿-欧拉动力学及其分散自适应控制的并行方案.国防科技大学学报,1988;10(4):45—55

[5]田志祥,吴洪涛.闭环双臂空间机器人的动力学建模与仿真.华南理工大学学报(自然科学版),2011;39(8):42—47

[6]贠超,等译.机器人学导论.北京:机械工业出版社,2006

[7]张继元.机构分析与综合的解.北京:人民交通出版社,2007

自适应迭代学习控制篇3

在汽车工业领域, 许多电线端子上需要安装防水密封的圆柱形橡胶栓, 其对橡胶栓的外形质量要求较高, 圆柱面上不能有明显的缺陷。传统的橡胶栓的检测主要依靠人工分拣, 费时费力, 而且受人为因素的影响, 检测效果也不理想。因此, 利用数字图像处理技术对橡胶栓进行在线检测以保证橡胶栓的质量, 具有重要的意义。

橡胶栓是通过注塑机成形的, 常见的缺陷主要是圆周上的飞边和缺料的凹坑等, 而CCD一次拍摄只能取到一半的圆柱面, 就可能会漏过橡胶栓背面的缺陷。因此在在线检测中, 课题组巧妙地采用了一组螺杆传送机构[1], 该机构不仅推动橡胶栓往前运动通过CCD, 同时使橡胶栓随螺杆的旋转时本身也回转, 这样在通过CCD下方时隔一定时间连续拍摄两次或三次, 基本上可实现整个柱面的检测。

由于橡胶栓在螺杆上是自由的, 因此在CCD拍摄时其位置角度是随机的, 而在图像识别中, 采用的方法往往是将被检测图像与标准模板图像进行比对来识别缺陷, 因此在利用图像处理技术对橡胶栓的表面缺陷进行识别时, 首先必须对不同位置的橡胶栓图像进行快速自适应旋转, 使图像中橡胶栓的中心轴与标准模板中橡胶栓的中心轴一致。

图像的自适应旋转是数字图像处理中的一个重点和难点, 目前, 这方面的文献报道较少。图像自适应旋转采用的方法主要有惯性主轴法[2]和圆形编码匹配[3]的方法。前者先检测出图像的边缘, 再利用边缘信息计算图像的质心和惯性主轴角度, 然后绕质心旋转, 所以该方法受限于边缘检测的精度, 并且容易受光照等外界条件的影响。后者将图像在圆周方向等分, 并建立方向编码来解决图像旋转匹配的问题, 计算量很大, 速度较慢, 难以满足在线检测的图像处理实时性要求。

本研究提出一种基于迭代匹配的快速图像自适应旋转算法, 将最小距离匹配的方法引入图像旋转, 然后通过迭代计算, 精确快速地实现橡胶栓图像的自适应旋转。

1 最小距离匹配

1.1最小距离匹配的原理

最小距离匹配又称模板匹配[4], 是一种重要的图像识别方法。最小距离匹配主要用于判断两个物体的相似度, 越相似则距离越小, 反之越大。

在利用最小距离匹配前, 需要预先将n个标准对象的特征存储为n个标准模板, 然后将具体的待检测样品与已知的标准模板进行比较, 计算它们之间的距离, 若待测样品与第k (k∈[0, n]) 个模板之间的距离最短, 则判定待测样品与第k个标准对象是同一类的物体。

两个对象间的距离可以用如下的方法计算:

例如, 假设一个标准对象A可以抽象为一个m维的特征向量:

XA=[xA1xA2 … xAm]T (1)

同理, 任意的待检测对象B也可以抽象为一个m维的特征向量:

XB=[xB1xB2 … xBm]T (2)

则, 待检测对象与标准对象间的距离可表示为:

$d_{A B} = d (X_{A}, X_{B}) = [\sum_{i = 1}^{m} (x_{A i} - x_{B i})^{2}]^{\frac{1}{2}}$ (3)

由于上式涉及到乘方、开方等复杂运算, 可将距离公式简化为:

$D_{A B} = D (X_{A}, X_{B}) = \sum_{i = 1}^{m} | x_{A i} - x_{B i} |$ (4)

1.2特征向量的选取

利用最小距离匹配需要将标准对象抽象为一个特征向量, 然后存储为模板。对实际的待检测样品, 也需要用同样的方法将其抽象为一个特征向量, 然后才可以计算它与模板间的距离。因此, 特征向量的设计具有重大的意义。

目前, 对特征向量的选取没有固定和统一的方法;相反, 随着研究的深入, 选取方法越来越多, 如利用图像的几何矩构成特征向量[5]、利用图像的灰度信息构成特征向量[6]等。本研究根据橡胶栓的实际特点和检测要求, 以橡胶栓在图像中所占比重为特征值来构成图像的特征向量。如图1所示, 将整幅图像等分为M×N个小块, 每块的像素数为X×Y。对图像中的任一点像素, 令:

则图像中第m行第n列小块的面积比重为:

$\frac{S_{m, n}}{S} = \frac{\sum_{i = m \cdot X}^{(m + 1) \cdot X - 1} \sum_{j = n \cdot Y}^{(n + 1) \cdot Y - 1} g (i, j)}{X \cdot Y}$ (6)

对任意的m、n, 令k=M×n+m, 于是, 面积比重亦可表示为:

$x_{k} = \frac{S_{k}}{S} = \frac{S_{Μ \cdot n + m}}{S} = \frac{S_{m, n}}{S}$ (7)

其中, m∈[0, M-1], m∈[0, N-1], k∈[0, MN-1]。所以, 该对象的特征向量就可以表示为:

X=[x1x2 … xMN-1]T (8)

上述方法得到了对象的k维特征向量, 在检测过程中, 需要根据实际的需要, 选择适当维数的特征向量。特征向量的维数越多, 那么匹配的精度越高, 但计算的时间复杂度更大;反之, 维数越少, 则精度降低, 但计算速度更快。

1.3最小距离匹配在自适应旋转中的应用

根据最小距离匹配的原理, 在图像旋转中, 假设被存储为模板的标准对象的旋转角度为0°, 则当图像的旋转角度越接近0°, 则它与模板间的距离越小;反之, 而当旋转角度增大时, 它与模板间的距离也随之增大[7,8]。

在橡胶栓检测过程中, 假设图像中橡胶栓中心轴水平时橡胶栓的旋转角度为0°, 并将该图像存储为模板, 那么当被检测的橡胶栓图像的旋转角度从-90°～+90°时, 它与模板间的距离变化关系如图2所示。

从图2可见, 在图像的旋转角度从-90°～+90°的范围内, 用最小距离匹配的方法来实现橡胶栓图像的自适应旋转是可行的。

对旋转角度不在-90°～+90°范围内的橡胶栓图像, 可根据橡胶栓的具体特点, 在实际检测过程中, 先利用数字图像处理技术将橡胶栓图像水平镜像, 然后再利用最小距离匹配的方法来实现图像的自适应旋转。

2 迭代匹配过程

2.1原始迭代方法

根据最小距离匹配原则, 在-90°～+90°的范围内, 待检测对象旋转角度偏离0°越大, 它与模板间的距离越大, 所以, 在实际的检测过程中, 可根据检测精度要求, 指定一个正数C, 当图像与模板间的距离小于指定正数C时, 即可认为图像已经旋转到位。对当前的待检测图像, 可以先计算距离D0, 然后分别计算顺时针和逆时针转动θ角 (θ>0, 下文顺时针用正数表示, 逆时针用负数表示) 后的偏离距离D1和D2, 则D0、D1、D2的4种位置关系如图3所示。

根据D0、D1、D2的位置关系就可以进一步确定之后图像旋转的方向和角度, 直到图像与模板的距离小于指定正数C。具体的迭代步骤如下:

(1) 判断D0是否小于指定正数C, 若是, 则跳到第 (8) 步;

(2) 计算当前位置的距离D0和当前位置正反转动θ角 (θ>0) 后的位置1、位置2的距离D1、D2;

(3) 若D2>D0>D1, 则继续;若D2<D0<D1, 则θ=-θ, 继续;若D0最小, 则跳到第 (6) 步;

(4) 判断D0是否小于指定正数C, 若是, 则跳到第 (8) 步;否则, 转动θ角, 计算当前位置距离D0′;

(5) 若D0′<D0, 则令D0=D0′, 跳到第 (4) 步;否则令θ=-θ/2, 跳到第 (4) 步;

(6) 若D2>D1, θ=θ/2;否则, θ=-θ/2;

(7) 转动θ角, 计算当前位置距离, 并赋值给D0, 若D0小于指定正数C, 跳到第 (8) 步;否则令θ=-θ/2, 跳到第 (7) 步;

(8) 结束。

2.2简化的迭代方法

原始的迭代算法通过不断地调整旋转方向和旋转角度实现自适应旋转, 可确保旋转的精度, 但如果迭代次数非常多, 势必影响迭代运算的速度。所以, 本研究下面提出一种简化的迭代算法, 可有效提高运算速度。

(1) 计算当前的距离D0;

(2) 将图像旋转θ角 (θ>0) ;

(3) 计算新位置的距离D1;

(4) 若D1<D0, 沿原方向不断转动θ角;若D1>D0, 则沿反方向不断转动θ角;

(5) 不断旋转θ角, 直到新位置的距离D1又变大;

(6) 反方向转动θ/2角。

(7) 结束。

简化的迭代算法的实现过程如图4所示。仅经过7步的迭代运算即可将图像旋转到位。

简化的迭代算法比原始的迭代算法的迭代步骤少, 因此, 它具有更快的运算速度。但由图4可见, 简化的迭代算法最后得到的图像与模板的距离并不一定是最小的, 它与最小距离可能有一点偏差, 因此为保证旋转的精度, 在计算过程中需要合理地选取θ值, 保证最后的图像与模板的距离小于指定正数C。

3 实验与结果

为验证基于最小距离匹配的橡胶栓图像自适应旋转算法的有效性, 在实际检测条件下, 本研究分别采用传统的惯性主轴法和原始迭代法、简化迭代法对具有不同旋转角度的橡胶栓图像进行旋转匹配。

实验中相关参数如下:

(1) CCD摄像头的分辨率:640×480;

(2) CPU:Intel Pentium (R) 4, 主频2.80 GHz;

(3) 内存:512 MB;

(4) 光源:LED环形灯;

(5) 模板:旋转0°的标准图像, 100维特征向量。

在实验中, 笔者分别对旋转角度为0°、10°、30°、45°、60°和80°的待匹配的图像进行对比实验, 如图5所示 ( (a) ～ (e) 为不同旋转角度的待检测图像, (f) 为自适应旋转后旋转角度接近0°的图像) 。

若实验中取正数θ=0.5, 初始角度θ=5°, 则两种迭代法匹配所需要的时间、迭代次数及最后距离如表1所示。

由表1可知, 相对于传统的惯性主轴法, 本研究中的原始迭代法能有效地完成橡胶栓图像的自适应旋转, 且效率没有明显降低;而简化迭代法不仅能有效地完成了橡胶栓图像的自适应旋转, 而且计算速度也有显著的提高。但是, 相对于原始迭代法, 简化迭代法的自适应旋转匹配较快, 但精度稍低 (最后距离大) 。因此, 在实际工业检测中, 可根据检测的精度要求选择适当的迭代匹配方法。

4 结束语

本研究将模式识别中的最小距离匹配引入图像的自适应旋转, 并利用迭代法实现匹配过程, 在保证匹配精度的同时, 也具有较快的运算速度, 因而在橡胶栓的检测中得到了很好的旋转匹配效果, 为后续橡胶栓图像表面缺陷识别提供了基础。另外, 笔者提出了用原始的迭代方法和简化的迭代方法分别实现自适应旋转匹配, 两者的侧重点有所不同, 可以根据实际的要求选择合适的算法。

参考文献

[1]魏燕定, 周豪, 陈波扬.在线橡胶栓外形质量检查机:中国, 101349550[P].2009-01-21.

[2]李金泉, 徐正华, 胡晓飞, 等.染色体质心、惯性主轴及旋转算法的研究[J].信号处理, 2001, 17 (3) :274-277.

[3]FARHAN U, SHUN-ICHI K.Using orientation codes for ro-tation-invariant template matching[J].Pattern Recogni-tion, 2004, 37 (2) :201-209.

[4]MCQUEEN M P C.A generalization of template matchingfor recognition of real objects[J].Pattern Recognition, 1981, 13 (2) :139-145.

[5]刘明, 叶正麟, 陈作平.基于二维特征向量的快速分形编码方法[J].计算机工程与应用, 2007, 43 (8) :82-84.

[6]陈玲, 沈红标, 李咸伟, 等.改进的图像纹理检索方法在矿石识别中的应用[J].中国图象图形学报, 2006, 11 (11) :1700-1703.

[7]RAFAEL C G, RICHARD E W.Digital Image Processing[M].2nd ed.London:Prentice Hall, 2007.

自适应迭代学习控制篇4

由于图像采集终端的性能影响、图像传输信道的干扰问题以及图像存储介质和展示介质的粒噪声等问题,不同类型和程度的噪声广泛存在于各种图像当中。这些噪声对基于图像的应用以及人类对图像信息的理解造成了困难,因此降噪是图像处理领域的重要问题。图像检索、人脸识别、增强现实以及图像数据融合等多个课题都与降噪的效果有密切关系。随着Web时代图像成为大数据越来越重要的数据形式[1],图像降噪成为网络数据处理领域不容忽视的问题。

随着图像处理应用范围的拓展和程度的深入,降噪技术引起了越来越多的关注。边界判别噪声检测方法是一种基于选择性中值滤波的降噪算法[2]。该方法通过扫描每个像素点的邻域,逐步缩小被损坏的邻域范围,从而确定噪声点。该方法可通过迭代使降噪性能逐步达到很高的水平,然而算法的复杂度和时间开销也是较大的。DT-CWT是一种基于离散小波变换的降噪方法,该方法通过重尾分布的二元分布模型提高了降噪的性能[3],但算法的收敛条件问题没有很好地解决。M3TM是一种快速的图像降噪方法,由于优化了马尔科夫树模型的求解流程,该算法取得了不错的实时应用效果[4]。不幸的是,该算法应用于医学影像以外的很多图像降噪的性能不佳。

以上各类算法多通过迭代提高降噪的效果,在这类算法中,迭代的终止条件对算法的时间开销和降噪的性能有至关重要的影响。在文献[2]中,迭代的初始区域选择为21×21像素,当迭代区域为3×3像素仍然有被损坏时,直接认定为噪声像素,这种方式使迭代的终止无弹性可言。伦敦大学帝国理工学院的Chourmouzios教授和希腊信息技术学院的Maria博士提出了一种图像降噪算法自动终止策略[5],然而这种策略依然不是根据图像自适应的。华南理工大学的肖宿等人提出了一种基于稀疏正则优化的降噪方法,在该方法中,通过引入迭代软阈值实现了自动终止[6]。然而该方法实现自动终止迭代所付出的额外计算开销较大。波尔多大学的Sutour等人提出了一种基于非本地均值的图像、视频通用降噪方法,该方法邻域像素加权实现迭代自终止[7]。实验结果反映出这种方法存在一定的过早结束迭代问题。

本文提出一种根据具体图像自适应停止迭代的图像降噪终止策略。该方法通过判断每一步迭代邻域内被还原的像素数目知否最小化而决定是否停止迭代,实现了降噪后的图像具有最大峰值信噪比。

1基本降噪方法

为了使关于本文提出迭代终止策略的讨论有所依托,在本文中采用基于选择性中值滤波技术的方法作为降噪算法。

1.1选择性中值滤波

设I是含有随机脉冲噪声的图像,Iij代表坐标为(i,j)处的灰度值。设ω为像素点(i,j)处的滑动窗,其尺寸满足由式(1)定义:

其中σ是大于零的整数,则选择性中值过滤器可由式(2)定义:

其中代表滑动窗ωi,j的均值,δ是选择阈值。如果用滑动窗内参考像素与目标像素均值的倒数作权重,则可以由滑动窗的加权标准差来计算阈值。加权后的均值Mw(i,j)定义如式(3):

定义加权标准差s(i,j)如式(4):

由上述步骤,可得到最终的降噪规则如式(5)所示:

其中α是给定的参数。通过迭代计算,就可以求得降噪之后的图像。

1.2传统迭代求参方法

通常式(5)的参数α在初始化时都被设置为一个较大的值,然后通过迭代求解参数,参数设定与求解如式(6)[8]:

在这种情况下,算法一般都无法自动求得调出迭代的条件,只能依靠人工预先设置的经验参数。

为了提高算法的性能,还有学者提出了其他迭代方法,如式(7)[9]:

在上述迭代方法中,依然没有提出循环终止的自适应策略,仍需依赖人工经验。

2自适应迭代终止策略

2.1降噪细节分析

为了引出文本的自适应迭代终止策略,需要精确分析基于选择性中值滤波器的随机脉冲噪声去除过程。

首先,如式(2)所定义的阈值,首先会被设置为一个较大的数值,只有被探测到的噪声像素被清除出图像。经过初始化之后,阈值会随着迭代次数的增多而逐渐下降。随着阈值的降低,降噪算法就会逐渐表现出去除非噪声像素的危险倾向,这时就需要终止迭代过程。

为了更形象地说明这一现象,图1给出了过度迭代和迭代不足的示例。在图1中,原图和含噪声图的尺寸都为20×20像素。

在图1中,(a)表示无噪声的原图,(b)是含噪声的图像,(c)是降噪过程中,迭代不足的情况,(d)表示降噪过程中,迭代过度的情况。(c)表明当迭代不足时,图像中会残留未被处理的噪声像素;(d)表明当迭代过度时,图像中会有非噪声像素被误清除。因此,合适的迭代次数,应当在(c)和(d)之间。

2.2 PNSR最大化迭代终止条件

遗憾的是,在实际应用中,人们无法知道哪些是噪声像素,哪些是非噪声像素,因此无法通过类似图1的观察来决定迭代是否终止。所幸在每一轮降噪过程中很容易确定被处理的像素点数目。在基于选择性中值滤波的降噪过程中,每一步被处理的像素点数目如图2所示。

图2表明,降噪算法处理迭代点的总数目是一个先减后增的函数。并且图2明显表现出当迭代中处理的像素点数目达到最小值时,降噪算法即将开始清除非噪声像素。

为了更准确地建立迭代终止判据,我们需要降噪的质量评价标准。为此引入了恢复图像的峰值信噪比。根据以往大量实验的结果,峰值信噪比在迭代中处理的像素点数目有递减变为递增之后的第二轮达到最大值[10],并且通常峰值信噪比越大,图像的质量越高。

由于迭代中处理像素点数目与峰值信噪比的关系,以及峰值信噪比与降噪质量之间的关系已经确立,因而可以得到本文提出的迭代终止策略的完整流程,具体算法如下:

图像降噪的自适应迭代终止算法

3实验

在实验部分,我们以FSWA[11]和BDND[12]为降噪算法,测试了使用本文提出的迭代终止条件的降噪效果,并与算法本身迭代终止方法所得到的降噪效果作了对比。实验的数据为经典图像处理测试材料Barbara和Lena图像,以及数据堂提供的图像处理数据集。

3.1图像降噪示例

首先以Barbara为实验图像,对比了基于本文提出迭代终止策略的FSWA和BDND算法降噪效果。实验结果如图3所示。

在图3中,(a)是Barbara的原始图像,(b)是含有40%噪声的图像,(c)是通过FSWA降噪的图像,(d)是通过BDND降噪的图像。图3表明基于本问题提出迭代终止策略的降噪算法具有明显的图像还原能力,特别是基于本文提出策略的FSWA算法,具有优秀的降噪能力。

3.2实验数据分析

为了定量分析本文提出迭代停止策略对降噪算法性能的影响,我们在人工添加20%高斯白噪声的图像数据集中,对比了FSWA算法和BDND算法在原有迭代模式下和在基于本文提出迭代终止策略下的降噪性能。着重对比了不同迭代终止策略下,相同图像集的平均未处理噪声像素数、平均误处理非噪声像素数和平均总有效处理像素数,实验结果如表1所示。

如表1所示,对两种降噪算法,本文提出的迭代终止策略都能提高它们的降噪效果。表1解释了FSWA的问题是有可能让迭代提前结束,BDND的问题是有可能过度迭代。

为了更直观地对比不同迭代终止策略对降噪性能的影响,我们对比算法应用不同迭代模式降噪后的图像峰值信噪比,实验结果如表2所示。

表2表明本文提出的迭代终止策略能切实提高图像降噪后的峰值信噪比,从而提高图像降噪的效果。

4结语

降噪是图像处理领域的基础问题,特别是随着物联网的发展,越来越多的图像被廉价终端所采集。并且由于图像传输信道与存储介质甚至展示设备的影响,噪声在各类图像数据集里是广泛存在的。现有大部分图像降噪算法是通过对像素邻域的迭代处理实现图像还原的。这一类方法面临的一个挑战就是如何根据具体图像的情况,自适应地终止迭代。

针对这一问题,本文提出了一种与图像相关的降噪算法的迭代终止策略。该策略通过对比每一轮降噪所处理的像素点数来发现处理像素数目曲线的变化趋势。当处理数目逐渐减少时,认为噪声像素点也在逐渐减少;当处理像素数量又开始增多时,说明算法已经开始误处理非噪声像素了,此时需要终止迭代。通过本文提出的迭代终止策略,算法可在尽量不处理非噪声像素的前提下最大化清除噪声像素,从而提高了降噪后图像的峰值信噪比,达到了更高的降噪效果。

摘要：降噪是图像处理领域的重要问题。针对现有大多数迭代类图像降噪算法终止条件判断复杂或不恰当终止的现象,提出一种与图像相关的迭代终止判据,不需要预设经验参数,而通过对图像数据的学习实现自适应的方法。该方法通过判断每一步被降噪的像素点的数量是否最小化来决定是否终止迭代,从而使处理后的图像具有最大峰值信噪比。实验结果表明,该方法在各种噪声条件图片的降噪中都实现了峰值信噪比最大化,是一种高效的、与图像相关的方法。

关键词：降噪,迭代终止条件,自适应,峰值信噪比,选择性中值滤波器

参考文献

[1]侯刚.基于内容的图像检索中特征表示与检索策略研究[D].吉林大学,2014.

[2]Wan Y,Chen Q,Kang Y.Robust impulse noise variance estimation based on image histogram[J].IEEE Signal Processing Letters,2010,17(5):485-488.

[3]Hilla P R,Achima A M,Bulla D R,et al.Dual-Tree Complex Wavelet Coefficient Magnitude Modelling Using the Bivariate Cauchy-Rayleigh Distribution for Image Denoising[J].Singnal Processing,2014,105(12):464-472.

[4]Hanzouli H,Lapuyade-Lahorgue J,Monfrini E,et al.PET/CT image denoising and segmentation based on a multi observation and a multi scale Markov tree model[C]//IEEE Nuclear Science Symposium and Medical Imaging Conference,Seoul,Korea,2013:1-4.

[5]Tsiotsios Chourmouzios,Petrou Maria.On the Choice of the Parameters for a Nisotropic Diffusion in Image Processing[J].Pattern Reconintion,2012,46(5):1369-1381.

[6]肖宿,韩国强.基于稀疏正则优化的图像复原算法[J].计算机应用,2012,32(1):261-264.

[7]Sutour C,Deledalle C A,Aujol J F.Adaptive Regularization of the NLMeans:Application to Image and Video Denoising[J].IEEE Transactions on Image Processing,2014,23(8):3506-3521.

[8]Pesnel S,Akkoul S,Lédée R,et al.Use of an image restoration process to improve spatial resolution in bioluminescence imaging[J].Molecular Imaging,2011,10(6):446-452.

[9]Wei C,Yang C,Lin C.A latent semantic indexing-based approach to multilingual document clustering[J].Decision Support Systems,2012,45(3):606-620.

[10]邱宇,基于双边滤波的图像去噪及锐化技术研究[D].重庆大学,2011.

[11]Varghese Justin,Ghouse Mohamed,Subash Saudia,et al.Efficient adaptive fuzzy-based switching weighted average filter for the restoration of impulse corrupted digital images[J].IET Image Processing,2014,8(4):199-206.

自适应迭代学习控制篇5

1 材料与方法

1.1 模型制作

分离一段长约30 cm的新鲜猪结肠, 清洗后将肠壁外翻, 于黏膜面钳夹小块黏膜组织, 用医用丝线结扎其根部形成模拟息肉, 息肉以一定间隔连续排列, 共30个, 其中息肉高度1.0~3.0 mm 10个, 3.1~5.0 mm 10个, 5.1~10.0 mm 10个, 见图1。然后将黏膜翻入恢复原肠壁结构, 注入空气使肠管充分膨胀后结扎肠管。扫描前将结肠模型固定于塑料容器内, 并向容器内注满0.1%碘水溶液。

图1模拟息肉。将离体猪结肠肠壁外翻, 于黏膜面用医用丝线结扎其根部形成连续线性排列的模拟息肉

1.2 CT检查

采用GE Discovery CT 750HD宝石能谱CT机扫描结肠息肉模型, 分别以不同管电流 (10、20、30、40、50、60、70、80、90、100、120、140、160、180、200、240、320 m A) 重复扫描。其他扫描参数:管电压120 k V, 球管旋转时间0.6 s, 螺距0.984∶1, 层厚0.625 mm, 层间距0.625 mm, 矩阵512×512, 重建视野50 cm。

1.3 图像重建

CT扫描结束后, 对不同管电流扫描条件下所得图像均以FBP、50%ASIR、100%ASIR和MBIR算法进行重建, 软组织算法, 重建层厚0.625 mm, 层间距0.625 mm。

1.4 图像分析

将重建后的数据传至GE ADW 4.3工作站进行图像后处理, 获得CT虚拟结肠镜 (CT virtue endoscopy, CTVE) 、结肠虚拟分割 (virtual dissection, VD) 和多平面重组 (MPR) 图像。由2名有CTC诊断经验的放射科副主任医师采用盲法对每组图像质量以4分制[5,6]进行评分, 评分标准见表1, 并记录每组图像的目测息肉检出率。取2名医师评分结果及息肉检出率的平均值进行分析。2名医师分别对每组图像固定选取2个层面测量图像的噪声和对比噪声比 (CNR) 。将圆形感兴趣区放置在容器上方空气内, 面积为1 cm2, 将被测空气密度的标准差 (SD) 作为图像噪声。CNR为息肉与空气的CT值之差与噪声的比值。噪声和CNR结果取2名观察者的平均值, 再分别计算ASIR和MBIR相对于FBP的噪声降低率和CNR提高率, 以MBIR为例, 噪声降低率MBIR= (SDFBP-SDMBIR) /SDFBP×100%;CNR提高率MBIR= (CNRMBIR-CNRFBP) /CNRFBP×100%。

1.5 统计学方法

采用SPSS 17.0软件, 各组图像噪声及CNR进行析因设计的方差分析, 对不同管电流水平的4种重建方法比较采用LSD法, P<0.05表示差异有统计学意义。

2 结果

2.1 FBP、ASIR及MBIR对图像噪声的影响

4种重建模式下的图像噪声均随管电流减小而增加, 其中MBIR增加幅度最小, FBP增加幅度最大 (图2A) 。MBIR及ASIR图像噪声较FBP明显下降 (F=21.860, P<0.05) ;管电流10~70 m A时, MBIR图像噪声均较FBP、50%ASIR及100%ASIR明显下降 (P<0.05) 。随着管电流降低, 50%ASIR及100%ASIR的噪声降低率无明显变化, 而MBIR的噪声降低率逐渐增加 (图3A) , 提示MBIR及ASIR均较FBP降低图像噪声能力强, 在低剂量条件下MBIR较ASIR降低图像噪声能力更强。

2.2 FBP、ASIR及MBIR对图像CNR的影响

4种重建模式下的图像CNR随着管电流减小而减小, 其中100%ASIR减小幅度最大, FBP减小幅度最小 (图2B) 。MBIR及100%ASIR图像CNR较FBP明显提高 (F=4.209, P<0.05) ;10~80 m A管电流时, MBIR图像CNR较FBP、50%ASIR及100%ASIR明显提高 (P<0.05) 。随着管电流降低, 50%ASIR及100%ASIR的CNR提高率无明显变化, 而MBIR的CNR提高率逐渐升高 (图3B) , 提示MBIR及ASIR图像CNR较FBP提高, 在低剂量条件下, MBIR图像CNR较ASIR及FBP提高。

2.3不同重建算法对结肠息肉的检出情况

利用4种重建算法获得CT仿真内镜及虚拟分割图像在所有管电流条件下息肉检出率均为100%, 3 mm以下息肉在重建获得的MPR图像上由于表现为微小隆起而不易与黏膜皱襞的局部隆起相鉴别, 而3 mm以上息肉在MPR图像上均可以显示。

图2不同重建算法的噪声 (A) 及CNR (B) 与管电流变化的关系

图3噪声降低率 (A) 和CNR提高率 (B) 与管电流变化的关系

2.4 FBP、ASIR及MBIR图像质量评价

10 m A管电流时FBP、50%ASIR、100%ASIR及MBIR 4种重建算法的平均主观评分值分别为 (2.00±0.77) 分、 (2.50±0.92) 分、 (2.50±0.50) 分、 (3.00±0.63) 分, 所得图像符合诊断要求的最低管电流及重建算法为10m A并应用MBIR重建算法。MBIR的CTVE及VD图像比FBP、50%ASIR及100%ASIR的图像对黏膜皱襞细节的显示更为清晰, 黏膜皱襞更为光滑。100%ASIR的MPR图像可见“勾边样”伪影, 而FBP、50%ASIR及MBIR的MPR图像中无此伪影, 且MBIR的MPR图像比FBP、50%ASIR及100%ASIR的图像更细腻, 息肉边缘更为光滑锐利, 病变密度更为均匀 (图4) 。管电流10 m Av及以上时, 4种重建算法的息肉检出率均为100%。

图4扫描电压为120 k V、管电流为10 m A时不同重建算法对图像质量的影响。A~D分别为FBP、50%ASIR、100%ASIR和MBIR的图像重建结果, 各组图像的横列从左至右为CTVE图像、VD图像和MPR图像 (轴位) 。MBIR的CTVE、VD图像较FBP、50%ASIR和100%ASIR对于黏膜皱襞细节显示更为清晰, MBIR的MPR图像息肉边缘较其他3种算法更为光滑锐利, 息肉密度更为均匀, 且没有100%ASIR图像中的“勾边样”伪影

3 讨论

3.1 ASIR和MBIR的原理

ASIR算法使用FBP数据建立一个最终结果, 假设邻近体素之间的噪声差异是由统计波动或图像噪声引起的, 对所有投射数据进行多次选择迭代和矫正, 由此重建出高质量和低噪声的图像。ASIR百分比反映了原始FBP图像与本质上完全通过数学模型建立的无噪声图像的混合比例, 例如0%ASIR实际上即FBP图像[7,8]。

MBIR是一个完全的迭代重建算法, 使用多次迭代和多模型, 除建立系统统计模型之外, 还建立了系统光学模型, 体素、X线光子初始位置和探测器几何因素均通过模型进行模拟, 真实地还原了X线从投射到采集的过程, 提高了重建图像的空间分辨率。MBIR用于重建的模型包括焦点模型、X线锥形束模型、三维体素模型、探测器单元模型、系统噪声模型。MBIR对每一个体素进行精确描述, 并考虑体素受到特定大小的X线焦点照射后激活探测器单元的情况;随体素的位置不同, 与X线焦点的距离不同, 受照射后激活探测器单元区域的大小也会不同。MBIR通过迭代方法不断地去除原始数据中的统计噪声和光学模糊效应, 从而还原一个真实的物体。此外, MBIR还加入了去除硬化伪影和金属伪影的技术, 但MBIR重建所需计算量非常大, 重建时间较长, 目前尚未应用于临床[6,9,10]。

3.2 ASIR和MBIR改善图像质量比较

本研究表明, 无论常规剂量还是低剂量扫描条件下, MBIR和ASIR图像的噪声均比FBP图像低, 而在低剂量扫描条件下MBIR的图像噪声比ASIR更低;在各个剂量条件下, MBIR和100%ASIR图像的CNR均比FBP图像高, 在低剂量扫描条件下MBIR图像的CNR比ASIR更高。低剂量扫描时, MBIR的CTVE、VD图像与ASIR及FBP图像相比, 黏膜皱襞细节显示更为清晰, 降低了息肉检出的假阳性率;而MBIR的CTVE及VD图像在息肉显示方面与ASIR及FBP图像差别不大, 可能是由于后处理时重建CTVE及VD图像的软件应用了某些算法, 改变了图像质量;MBIR的MPR图像对息肉边缘显示较ASIR及FBP图像更为光滑, 病变密度更为均匀, 且没有100%ASIR图像中的“勾边样”伪影, 可见MBIR的图像质量比ASIR及FBP的图像质量高。在较高扫描剂量时, 尽管MBIR与100%ASIR相比不能明显降低图像噪声并提高图像CNR, 但100%ASIR的MPR图像中有“勾边样”伪影, 且临床上一般也不采用100%ASIR图像, 由于过高的ASIR (如100%) 噪声谱与常规不同, 高比例的ASIR可能会大幅改变图像的纹理和特性[11]。

本体外实验模型研究发现, 对息肉满足100%检出率的CTC最低剂量扫描方案为扫描电压为120 k V、管电流为10 m A, 最佳图像质量为应用MBIR, 与Yoon等[12]的研究结果一致。尽管MBIR的MPR图像对3 mm以下息肉的检出率不能达到100%, 但MBIR的CTVE及VD图像息肉检出率可达100%, 能满足诊断要求。

3.3 本研究的局限性

本研究的局限性在于: (1) 本研究样本量小, 且模拟息肉均为带蒂息肉, 未包括扁平息肉和无柄息肉, 而扁平息肉及无柄息肉不易与黏膜皱襞局部隆起相鉴别, 故此离体模型研究有可能高估息肉检出率。 (2) 本实验模型中息肉以一定间隔呈线性排列是为了方便对息肉进行观察, 比较息肉大小, 但是可能会降低息肉检出的假阳性率。 (3) 扫描中盛放肠管的容器体积固定, 不能代表不同体重指数的患者。 (4) 本离体模型研究属于理想条件下的研究, 未考虑肠蠕动和血管搏动等的影响, 并且不存在患者肠道清洁不满意的因素, 还需要在后续的体内研究中加以改进。 (5) 本研究未评估2名观察者对图像质量进行主观评价时观察者间的差异。

总之, 本离体模型研究结果表明在CTC中, ASIR和MBIR重建算法与FBP算法相比, 均可以显著降低图像噪声, 提高图像质量。在低剂量扫描条件下, MBIR提高图像质量的能力比ASIR更强, 为低剂量条件下进行CTC筛查结肠疾病提供了可能。

参考文献

[1]Svensson MH, Svensson E, Lasson A, et al.Patient acceptance of CT colonography and conventional colonoscopy:prospective comparative study in patients with or suspected of having colorectal disease.Radiology, 2002, 222 (2) :337-345.

[2]Laghi A, Iafrate F, Rengo M, et al.Colorectal cancer screening:the role of CT colonography.World J Gastroenterol, 2010, 16 (32) :3987-3994.

[3]Gunn ML, Kohr JR.State of the art:technologies for computed tomography dose reduction.Emerg Radiol, 2010, 17 (3) :209-218.

[4]Ulzheimer S, Endt H, Flohr T.Computed tomography—patient dose and dose reduction technologies.Health Phys, 2011, 100 (3) :325-328.

[5]宁培钢, 朱绍成, 史大鹏, 等.高级迭代重建算法降低腹部CT剂量的潜能:体模研究.中国医学影像技术, 2012, 28 (12) :2243-2247.

[6]吴瑶媛, 王万勤, 刘斌, 等.FBP、ASiR和VEO三种重建算法对常规剂量胸部CT图像质量的影响.中国医学影像技术, 2012, 28 (3) :575-578.

[7]王艳, 史大鹏, 朱绍成, 等.比较以自适应性统计迭代重建技术和滤过反投影重建的低剂量腹部CT的图像质量.中国医学影像技术, 2012, 28 (10) :1902-1905.

[8]于晓坤, 孙浩然, 李又洁, 等.ASIR提高低剂量CT结肠成像息肉检出率有效性的体模研究.临床放射学杂志, 2011, 30 (10) :1536-1540.

[9]Chang W, Lee JM, Lee K, et al.Assessment of a model-based, iterative reconstruction algorithm (MBIR) regarding image quality and dose reduction in liver computed tomography.Invest Radiol, 2013, 48 (8) :598-606.

[10]Miéville FA, Gudinchet F, Brunelle F, et al.Iterative reconstruction methods in two different MDCT scanners:physical metrics and 4-alternative forced-choice detectability experiments--a phantom approach.Phys Med, 2013, 29 (1) :99-110.

[11]Silva AC, Lawder HJ, Hara A, et al.Innovations in CT dose reduction strategy:application of the adaptive statistical iterative reconstruction algorithm.Am J Roentgenol, 2010, 194 (1) :191-199.

自适应迭代学习控制篇6

奈奎斯特采样定理要求对信号进行采样时采样频率要大于或者等于信号带宽的两倍。而这样得到的数据具有较大的冗余。压缩感知技术的出现，打破了采样定理的限制。压缩感知技术表明：只要信号是稀疏的或可压缩的，则数据采样可以低于奈奎斯特采样频率，而且原始信号可以由采样值精确的恢复。从而大大减少了数据的处理量。通过挖掘信号的部分先验信息，可以减少观测值方面的工作。本文中利用一部分的大系数的稀疏系数的位置设置了一个阈值矩阵，之后这个矩阵与高斯随机矩阵构成自适应观测矩阵，达到改进重构性能的目的。

2 压缩感知的原理：

标准的压缩感知的数学模型为：

其中θ为未知信号x在稀疏域Φ中的系数，y是观测矩阵Ψ对x进行先行观测后的得到的观测向量。若θ中的非零的向量为K(K<

在此选择阈值迭代算法(IHT)解决此问题。

通过挖掘信号的的部分先验信息，可以减少观测值的长度，当观测值的维数及重构的算法相同时，信号稀疏度越高，重构的精度越高。显然信号的未知支撑集合比原始信号更稀疏。利用自适应的观测矩阵，可以将稀疏域的小系数更接近于零，从而优化重构建效果。

3 观测矩阵

观测矩阵的设计性能直接影响压缩感知的质量，观测通过观测矩阵Ψ对信号x进行映射，或者压缩感知矩阵T对θ进行线性投影，若T满足受限等距特性(RIP)，则此时L0可以由θ重构，从而的到信号x的重构。自适应观测是借助信号x或者稀疏系数θ的某些先验特性，生成相应的观测矩阵，这类矩阵性能上显著由于随机观测矩阵和确定观测矩阵。文中以高斯随机观测矩阵为原始矩阵由θ的稀疏特性对矩阵的特性进行改进，使得改进后的矩阵满足受限等距特性(RIP)的同时，使得θ的小系数更接近于0。

3.1 观测矩阵的构造

设W=supp(θ);K=|W|，W1是信号的已知支撑集部分，W2是信号的未知支撑集部分，K1=|W1|，W=W1+W2，观测矩阵的自适应过程如下：

1)生成M×N的高斯随机矩阵，作为原始矩阵Ψ;

2)将Ψ最大值归一化并且按照已知和未知支撑部分分为两块;

3)选取θ中的大系数的K0个并把他们的位置记录在矩阵A中，记A为{a1,a2，…ak0};得到域阀矩阵，Φ是信号的稀疏域变换矩阵。

∆中的元素的选取规则是：若元素的的行列数不等则为零，若元素行列数相等且列数位于大系数矩阵A中，则此元素等于1，除此之外的元素赋值为极小数10-q其中q<-10，新的自适应观测矩阵是不难看出我们使用θ的部分大系数的位置信息W1，从而使优化过的自适应观测矩阵的行数缩减了K0行，从而减少了K0个观测值。在保证才数据传输量变化不大的前提下，应用后面相应的重构算法，使得应用自适应矩阵的重构效果明显好于应用高斯随机矩阵的重构效果。

3.2 压缩感知矩阵的受限等距特性（RIP)

应用自适应观测矩阵代替原有矩阵于(1)中，可知压缩感知模型变化为：

其中，,通过对压缩感知系数θ进行线性观测后可以得到观测矩阵。且当受限等距常数δk满足：若C∈(cj)j∈Γ,|Γ|≤K时有

则称压缩感知矩阵T满足RIP。

由于T、都是高斯随机矩阵，故可的RIP常数和T的RIP常数关系有有。的RIP常数的RIP常数大小关系为

4 自适应观测矩阵对压缩感知重构性能的改进

下面我们介绍使用的阈值迭代算法。IHT算法源于梯度下降算法，其采用高斯随机矩阵Ψ最为观测矩阵，传感矩阵为T=，从而观测向量值y=Tθ+e(e为观测噪声)。在运用自适应矩阵代后，传感矩阵变为，观测值向量为,IHT的迭代算法如下：

其中分别为第n+1和第n次迭代后得到的重构稀疏系数。因为可得：。IHT算法具有实现速度快，重构效果好的特点。下面对于自适应观测矩阵对IHT算法重构性能方面的改进作具体的分析。

4.1 自适应观测矩阵对IHT重构误差的改进

对于IHT重构的性能判断，其中一条是其重构误差的大小。IHT的重构误差的定义式：

其中β2k与T2k的RIP常数δ2k的关系为，应用自适应观测矩阵后IHT的重构误差为，根据参考资料[4]可知，对于由于所以的到以下结论：

可以看出，使用自适应观测矩阵后的IHT重构误差会明显小于用高斯随机观测矩阵的IHT重构误差。

4.2 仿真结果

对于信号经稀疏分解后的系数的组成往往是有限个的大系数和其余的小系数，而往往这些小系数并不为零，只是相对较小，有些甚至和零相差较大，而这种近视稀疏必然会降低压缩感知的性能。当采用自适应观测矩阵对稀疏系数θ进行观测时，可以看作是压缩感知矩阵进行线性观测。从而的稀疏性要优于θ的稀疏性。在仿真中，稀疏分解采用离散余弦变换(DCT)，重构算法采用IHT算法，信噪比取多帧的平均值。仿真软件为Matlab 2011b，分别对灰度图像rice,peppers和football进行实验。结果表明应用自适应观测矩阵的到的重构图像效果要明显优于使用高斯随机观测矩阵重构的图像。

5 小结

本文利用自适应观测矩阵，对信号进行重构。自适应观测矩阵以高斯随机观测矩阵为原始矩阵，结合信号的稀疏稀疏的部分先验信息，利用小部分的大稀疏系数的位置信息设计了一个域阀矩阵。在运算的同时减少了观测矩阵的行数，是的数据的传输代价增加较小。最后利用DCT变换和IHT算法进行仿真，证明了自适应观测矩阵的重构性能要优于高斯随机矩阵的重构性能。

关键词：自适应观测矩阵,压缩感知,阈值迭代算法

参考文献

[1]John G.Proakis,Dimitris G.Manolakis.数字信号处理[M].电子工业出版社,第四版方艳梅等译,2007年

[2]Katsuhiko Ogata.现代控制工程[M].电子工业出版社,第五版,卢伯英等,2011年

[3]U.Meyer-Baese.数字信号处理的FPGA实现[M].清华大学出版社,第三版,刘凌译,2011年版

自适应迭代学习控制篇7

目前, CT已经成为一种主要的医学放射暴露源。最新的流行病学证据表明, 高剂量头颅CT的辐射会增加患白血病和颅脑肿瘤的风险[1,2,3], 所以寻找一种降低CT剂量的新方法是目前迫切需要解决的问题。CT最常用的重建算法是滤过反投影 (filetered back projection, FBP) 重建算法, 该算法的特点是:当降低放射剂量时, 图像信噪比会随之下降。而自适应统计迭代重建 (adaptive statistical iterative reconstruction, ASIR) 算法是最新的一种技术, 与FBP算法相比, ASIR最大的优势是能在常规FBP算法表现为低信噪比时明显提高图像质量, 也就是说ASIR能在低剂量扫描时提高图像信噪比, 因此ASIR为临床低剂量扫描开辟了一种新途径。

对于ASIR算法, 前人在胸部[4,5,6]、腹部[7,8,9]、心脏[10,11]及大血管[12]方面低剂量扫描进行了大量的研究并报道, 而在颅脑方面的应用报道较少。本研究旨在发现ASIR在颅脑扫描中降低辐射剂量的证据。

1 资料与方法

1.1 研究对象

受试者来自于2011年3月至2013年5月解放军总医院门诊就诊患者, 实验经解放军总医院医学伦理委员会批准。受试者均签署知情同意书。

纳入标准:CT检查示脑内无明确病变, 能顺利并配合完成检查, 签署知情同意书。

排除标准:CT检查结果证实有以下疾病: (1) 脑血管病, 如脑梗塞或脑出血病史等; (2) 脑肿瘤, 如脑原发肿瘤、转移瘤; (3) 脑外肿瘤, 如脑膜瘤等; (4) 脑外伤; (5) 异物; (6) 其他颅内疾病。

1.2 研究方法

1.2.1 数据采集

采用型号为Discovery CT 750 HD宝石能谱CT, 对32例受检者进行头部扫描, 年龄 (74±16.48) 岁。扫描条件120 k V、220 m A, 球管旋转1 s。轴扫, 间隔:5 mm, 层厚:5 mm, 共24张图像。

1.2.2 数据测量与处理

扫描结束后在双侧侧脑室最大层面重建图像, FOV 25 cm、层厚5 mm、间隔20 mm, 共5层, 分别将ASIR值设定为0%、20%、40%、60%、80%、100%重建图像, 窗宽为100 HU、窗位35 HU。在一侧侧脑室内选取圆形感兴趣区 (region of interest, ROI) (面积为9.9 mm2, 周长为7.8 mm, 像素点为30) , 测量不同的ASIR值时侧脑室脑脊液的CT值与标准差 (如图1所示) 。由2名高年资放射诊断科医生对数据进行测量, 每一个数据测量2次后取平均值, 采用双盲的方法记录数据。

1.2.3 统计学分析

应用统计学软件SPSS 13对数据进行统计学分析, 首先检验病灶的CT值和标准差等数据是否符合正态分布, 之后对数据进行t检验。

2 结果

当ASIR值为0%、20%、40%、60%、80%、100%时, 侧脑室脑脊液的CT值SD分别为 (4.706±1.192) 、 (3.466±0.872) 、 (2.8±0.761) 、 (2.284±0.618) 、 (1.784±0.547) 、 (1.425±0.507) , 各组间具有明显的统计学差异, P<0.01, 随ASIR百分比的增加, 图像噪声水平明显降低。

3 讨论

医用辐射已成为人类接受X线照射最大的来源。特别是近年来CT技术的发展, 病灶的检出率大大增加, CT应用的范围也越来越广泛, 导致CT检查次数在不断增加, 从而引起辐射剂量的增加。CT产生的电离辐射具有潜在致癌的风险, 会增加患白血病和颅脑肿瘤等恶性肿瘤发生的几率[1,2], 同时因CT扫描速度比较快, 所以头部CT扫描临床应用也越来越多, 进而带来的是放射剂量的危害也越来越严重。因此, 当前降低CT辐射剂量、提高有效剂量是目前迫切需要解决的问题。

随着社会发展进步和人们生活水平的提高, 人们越来越意识到放射防护的重要性。人们之前用了很多方法尝试降低头部扫描剂量[13,14,15,16], 如降低曝光量、遮挡腺体和眼睛等重要部位、自动球管电流调节、倾斜扫描架避免扫描眼眶等, 这些方法能在一定程度上降低CT扫描剂量。此外, 学者们还尝试在重建技术上降低辐射剂量, 其中, 临床上采用FBP重建算法时是通过降低曝光量而达到降低剂量的目的。但这种方法的局限性在于:该方法会导致图像噪声的增加和图像对比噪声比的下降。这是因为:FBP重建方法虽然速度较快, 但是该方法要求每次投影测量数据都是精确定量和完全的, 而FBP重建方法对噪声和伪影都很敏感, 当辐射剂量降低或投影数据采集不足时, 重建出的图像质量就会很差。因此, 为了保证图像质量, 采用FBP重建方法就不能大幅度地降低CT扫描的辐射剂量。

ASIR重建方法是一种较新的重建方法。该方法使用从FBP算法获取的图像信息作为图像重建的原始结构单元, 通过建立系统噪声模型并且利用迭代方法加以抑制, 从而得到更加清晰的图像。该方法是通过将FBP数据和ASIR数据按不同比例加权融合而达到不同程度的降噪效果。与FBP重建算法相比, AISR算法重建图像时图像噪声更低[17]。本研究就是依据这个原理进行研究的, 研究结果证实了AISR重建算法能有效地降低图像噪声。结果显示, 随着ASIR数据按不同比例加权的增加, 图像对比越来越平滑, 噪声越来越小。也就是说ASIR百分比的增加, 图像噪声水平明显降低。提高ASIR权重值能够降低颅脑CT图像噪声, 因此扫描时采用较高ASIR值重建图像具有降低射线剂量的作用。本研究结果为临床采用ASIR方法降低颅脑扫描剂量提供了理论依据。

本研究中, 随着ASIR算法比例加权的增加, FBP算法比例加权降低, 图像噪声降低。从这个角度上也说明, ASIR算法比FBP算法在降低图像噪声上更优越, 是一种更好的降噪手段, 对于降低辐射剂量更有效。在临床使用过程中, ASIR算法重建需要的时间很短 (据文献报道不会超过5 min[18]) , 不会影响临床的效率, 说明该算法有很好的临床可行性。

本研究采用图像中脑脊液作为测量对象, 这是因为有研究对ASIR算法和FBP算法2组脑组织噪声差异进行对比发现, ASIR算法组脑脊液图像噪声明显下降, 而脑白质、颅骨和空气的噪声对比2组没有差异[19], 原因是脑脊液密度接近于水的密度, 且CT值相对比较恒定, 不容易受其他因素的干扰。这里唯一需要考虑的是脑脊液流动伪影造成的影响, 而双侧侧脑室内脑脊液相对流速较慢, 就有效地避免了这个问题。颅骨和空气CT值变化较大, 同时与其他组织相邻会有界面征, 会给CT值测量带来麻烦。脑白质虽然没有颅骨和空气对于CT值影响那么大, 但脑白质本身在不同部位和纤维束走行都会影响测量结果。综上所述, 选择侧脑室脑脊液来测量CT值结果更科学。

【自适应迭代学习控制】推荐阅读：