Segmentation of Breast Cancer on DCE-MRI Images With MRF Energy and Fuzzy Speed Function
-
摘要: 乳腺癌灶的精确分割是乳腺癌计算机辅助诊断的重要前提. 在动态对比增强核磁共振成像(Dynamic contrast-enhanced magnetic resonance imaging, DCE-MRI)的图像中, 乳腺癌灶具有对比度低、边界模糊及亮度不均匀等特点, 传统的活动轮廓模型方法很难取得准确的分割结果. 本文提出一种结合马尔科夫随机场(Markov random field, MRF)能量和模糊速度函数的活动轮廓模型的半自动分割方法来完成乳腺癌灶的分割, 相对于专业医生的手动分割, 本文方法具有速度快、可重复性高和分割结果相对客观等优点. 首先, 计算乳腺DCE-MRI图像的MRF能量, 以增强目标区域与周围背景的差异. 其次, 在能量图中计算每个像素点的后验概率, 建立基于后验概率驱动的活动轮廓模型区域项. 最后, 结合Gabor纹理特征、DCE-MRI时域特征和灰度特征构建模糊速度函数, 将其引入到活动轮廓模型中作为边缘检测项. 在乳腺癌灶边界处, 该速度函数趋向于零, 活动轮廓曲线停止演变, 完成对乳腺癌灶的分割. 实验结果表明, 所提出的方法有助于乳腺癌灶在DCE-MRI图像中的准确分割.
-
关键词:
- 乳腺癌 /
- 动态对比增强核磁共振成像 /
- 马尔科夫随机场能量 /
- 活动轮廓模型 /
- 模糊聚类
Abstract: Accurate segmentation of breast cancer is an important step for computer-aided diagnosis. In images obtained from dynamic contrast-enhanced magnetic resonance imaging technique, the traditional active contour model method is difficult to obtain accurate segmentation results due to the low contrast, blurred boundary and intensity inhomogeneous of the breast cancer images. In this paper, a semi-automatic segmentation of active contour model combining Markov random field (MRF) energy and fuzzy velocity function is proposed to perform the segmentation of breast cancer in dynamic contrast-enhanced magnetic resonance imaging (DCE-MRI) images. This method has the advantages of fast speed, objectivity and the ability to reproduce the segmentation result compared to the manual segmentation of a professional doctor. First, the MRF energy of DCE-MRI is calculated to enhance the difference between the target area and the background. Second, the posterior probability of each pixel is calculated in the energy map. Then, region term of active contour model based on the posterior probability is developed. Finally, a fuzzy speed function, which derived by combining the image intensity, time domain characteristics of DCE-MRI and the Gabor texture feature, is introduced into the active contour model as edge function. At the boundary of breast cancer, the edge function approaches zero and the evolution of the contour curve will stop. The experimental results showed that the proposed segmentation method can accurately segment breast cancer in the images of DCE-MRI. -
预设性能控制(Prescribed performance control, PPC)是近些年提出的一种约束控制方法[1−3], 该控制方法在保证系统状态既满足稳定性要求又始终处于预设的性能函数之内的同时, 还可通过人为设定已知的边界函数刻画系统的动态性能, 从而解决实际中存在的状态约束及过程动态性能优化问题.
常用的预设性能控制方法由Bechlioulis等[4−5]提出, 其核心在于对于误差变量的非线性变换, 即通过双曲正切函数等手段将受约束状态转换映射为无约束形式, 然后采用常规的状态反馈、反步法等方法开展控制系统设计, 由于其易于与其他方法结合, 因此在近些年获得大量的关注. 在理论研究方面, 已用于解决有限时间控制问题[6]、固定时间控制问题[7]、控制方向未知问题[8] 以及容错控制问题[9] 等; 在实际应用方面, 已开始研究应用于高超声速飞行器姿态控制问题[10−11]、卫星姿态控制问题[12]、工业机电系统[13] 等. 预设性能控制方法存在的问题在于其边界函数选取有限制, 要求单调递减且终值不能为零. 而且, 不同的稳定性需求下预设性能控制需要重新设计, 方法也不同, 缺乏统一的设计框架.
通道控制[14]通过将状态与边界性能函数差的倒数作为控制增益, 同样可满足预设边界. 目前, 该方法已扩展到任意阶系统, 并实现了有界和渐近稳定性. 比如, 针对任意相对阶的不确定非线性多输入多输出系统的跟踪控制问题, 文献[15] 提出一种低复杂度的无模型控制器, 实现跟踪误差在预设性能漏斗通道内演化. 由于非最小相位不确定线性系统的零动态可能存在不稳定部分, 文献[16] 中的通道控制方法通过增加相对度来完全消除内部动力学中不稳定部分, 实现有界稳定性. 除此之外, 近年来以控制障碍函数[17]和障碍李雅普诺夫函数[18]为代表的安全控制方法, 通过将系统状态限制在安全区域内来保证系统在满足约束的同时趋于稳定, 其也能够实现预设性能的效果. 但上述文献的方法同样难以通过同一种控制方法实现不同的稳定性需求.
因此, 本文受预设性能控制与通道控制的启发, 针对一类不确定非线性系统提出一种性能驱动控制 (Performance-driven control, PDC). 不同于传统预设性能控制与通道控制以系统稳定为目标, 该方法通过基于Metzler矩阵的正系统理论和切换控制保证状态量与边界性能函数误差始终为正, 从而保证状态量始终在预设的性能函数范围内. 该方法的特色是: 1)系统的状态始终保证在预设的边界性能函数范围内. 2)系统稳定性由边界性能函数决定, 在不同的边界性能函数下通过同一种控制框架实现有界和渐近多种稳定性. 3)该方法采用的是系统状态与边界性能函数的差, 而非传统预设性能控制与通道控制中采用的系统状态与边界性能函数的商. 一方面放宽边界性能函数的选取范围, 允许边界终值为零; 另一方面保证系统状态在受突发扰动穿越边界后仍能返回预设的边界性能函数内. 由于系统稳态、动态等性能完全由边界性能函数决定, 因此称本文方法为“性能驱动控制”. 本文方法通过将状态量约束在预设性能函数之内, 可以有效实现期望的系统动态性能, 而且面对突发扰动可保证良好的鲁棒性, 对于航天器、飞行器等性能需求较高的实际物理对象控制系统设计具有较好的参考与应用价值.
本文后续内容安排如下: 第1节详细阐述问题及基本假设; 第2节给出性能驱动控制方法的主要理论结果并进一步开展方法讨论; 第3节仿真验证了方法有效性; 第4节是结束语.
1. 问题阐述
本文考虑一类如下形式的非线性不确定系统:
$$\left\{ \begin{aligned} {{\dot x}_1}\left( t \right) &= {x_2}\left( t \right)\\ {{\dot x}_2}\left( t \right) &= F\left( {x,\;t} \right) + B\left( {x,\;t} \right)u\left( t \right) + d\left( t \right) \end{aligned}\right. $$ (1) 其中状态$ x_1(t),\;x_2(t)\in \mathbf{R}^m $, 控制输入$ u(t)\in \mathbf{R}^m $, 不确定性$ d(t)\in \mathbf{R}^m $. 非线性时变函数$ F\left( {x,\;t} \right): {{\mathop{\mathbf{R}}\nolimits} ^{2m}} \times {{\mathop{\mathbf{R}}\nolimits} _ + } \to {{\mathop{ \mathbf{R}}\nolimits} ^m} $和$ B\left( {x,\;t} \right):{{\mathop{\mathbf{R}}\nolimits} ^{2m}} \times {{\mathop{ \mathbf{R}}\nolimits}_+ } \to {{\mathop{\mathbf{R}}\nolimits} ^{m \times m}} $在状态$ x $上满足局部Lipschitz条件, 在时间$ t $上是时变的分段连续有界函数. 本文中考虑$ x_1(t) $和$ x_2(t) $是可测的, 且$ F(x,\;t) $和$ B(x,\;t) $为已知.
首先给出如下的性能函数定义.
定义 1. 考虑函数$ \rho(t) = [\rho_1(t),\; \cdots , \;\rho_m(t)]^{\rm{T}}\in \mathbf{R}^{m} $满足如下性质: 1) $ \rho(t) $ 是光滑的连续有界函数; 2) $ \mathop {\lim }_{t \to 0} \rho_i (t) = {\rho _{i,\;0}},\;{\rm{ }}\mathop {\lim }_{t \to \infty } \rho_i (t) = {\rho_{i,\;\infty} } $, 其中$ \rho_{i,\;0}>\rho_{i,\;\infty}\ge0 $ 是常数, 分别表示函数的初值与终值. 那么, $ \rho(t) $ 称为性能函数.
定义1中, $ \rho_{i,\;0}、\rho_{i,\;\infty} $是影响性能函数特性的关键参数, 在后文详细介绍. 另外, 注意定义1与经典文献[4]定义是不同的, 其不要求终值$ \rho_{i,\;\infty} $ 必须为正, 非负也可. 这个区别使文献[4]中的预设性能控制方法不再适用. 其次, 对文中将出现的部分符号进行定义.
定义 2. 对于任意向量, 定义符号$ \succ、 \prec、 \succeq $、$ \preceq$表示向量中所有元素均大于、小于、不小于和不大于.
针对上述系统(1), 本文控制目标是设计控制律$ u(t) $使 $ \underline{\rho}\left( t \right)\preceq x_1(t) \preceq{\bar \rho} \left( t \right) $ 成立, 即保证状态$ x_1(t) $中每一个元素满足预设的性能约束, 其中函数$ \bar{\rho}(t),\; \underline{\rho}(t)\in \mathbf{R}^{m} $ 是边界性能函数.
假设 1. 边界性能函数$ \bar{\rho}(t) $和$ \underline{\rho}(t) $ 满足关系: $ \bar{\rho}(t)=-\underline{\rho}(t) = \rho(t) $, 其中$ \rho(t) = [\rho_1(t),\, \cdots , \,\rho_m(t)]^{\rm{T}} \in \mathbf{R}^{m} $是给定的性能函数.
一般也可将函数$ \bar{\rho}(t),\; \underline{\rho}(t) $称为边界函数, 表明其约束状态$ x_1(t) $ 的运动边界. 假设1表明本文主要考虑的是对称边界来约束系统的动态性能, 而这也是预设性能控制中常用的假设条件.
本文的控制目标并未明显提及系统的稳定性, 这与经典预设性能控制不同. 本文希望通过以上对状态的约束, 由性能函数决定系统稳定性. 比如说, 如果性能函数$ \rho(t) $有界, 那么状态$ x_1(t) $也是有界稳定的; 如果性能函数$ \rho(t) $渐近收敛, 那么状态$ x_1(t) $也是渐近稳定的. 这是本文提出的控制方法的巧妙之处, 也是与预设性能控制的最大区别之处.
假设 2. $ B(x,\;t) $在任意时刻和状态保持非奇异.
假设 3. 不确定性$ d(t)=[d_1,\; \cdots,\; d_m]^{\rm{T}} $是有界的, 即满足$ 0 \le |{d_i}| \le {D_i},\;i = 1,\;2,\;\cdots,\;m $, 其中$ D_i $ 是已知正常数, 表示不确定性的上界.
注 1. 控制系数$ B(x,\;t) $是方阵且非奇异实际已保证非线性系统(1) 是完全能控的. 假设2在机器人领域、航天器和飞行器等工程系统控制中广泛适用, 与文献[19]中的全驱概念是一致的, 具体介绍可参看文献[19]. 假设3 认为系统受到的不确定性是有界且匹配的, 这也是一个鲁棒控制中常用的假设条件.
2. 性能驱动控制
2.1 控制器设计及稳定性分析
首先给出一些必要的定义和引理.
定义 3. 若方阵$ S\in \mathbf{R}^{k \times k} $中任意非对角元素为非负数, 即$ {S_{i,\;j}} \ge 0,\;{\rm{ }}i \ne j \in \left[ {1,\;k} \right] $, 则该矩阵为Metzler 矩阵[20].
引理 1. 考虑非齐次线性系统$ \dot z\left( t \right) = Sz\left( t \right) + \theta \left( t \right) $, 其中状态向量$ z(t)\in \mathbf{R}^k $, $ \theta(t):{{\mathop{\mathbf{R}}\nolimits} _ + } \to {\mathop{\mathbf{R}}\nolimits} _ + ^k $. 系统满足性质: 1)方阵$ S $为Metzler 矩阵; 2)系统初始状态$ z(t_0)\succeq 0 $时, 方程的解满足$ z(t)\succeq 0,\; \forall t\ge t_0 $ [21].
引理1给出正系统性质, 即系统状态始终为非负. 注意正系统并不一定要求系统是稳定的, 因此矩阵$ S $不一定是Hurwitz 的. 然后, 给出本文提出控制器的具体形式. 定义$ \bar{e}_1(t) $和$ \underline{e}_1(t) $ 分别为状态$ x_1(t) $与上边界性能函数和下边界性能函数的误差, 即
$$ \begin{align} \bar{e}_1(t) = \bar{\rho}(t)-x_1(t) \end{align} $$ (2) $$ \begin{align} \underline{e}_1(t) = x_1(t)-\underline{\rho}(t) \end{align} $$ (3) 根据以上定义, 本文的控制目标可以转变为$ \bar{e}_1(t),\;\underline{e}_1(t)\succeq0,\;\forall t\ge0 $. 进一步, 定义变量
$$ \begin{align} \xi(t) = {x_2}(t)- {K_1}{x_1} (t) \end{align} $$ (4) 其中$K_1 = {\rm diag}\{K_{1,\;i}\}\in \mathbf{R}^{m\times m}\;(i=1,\;\cdots,\;m)$是对角阵, $K_{1,\;i} $是常数. 系统(1)的相对阶为2, 通过式(4)中对$ \xi $的定义, 使得$ u $ 可直接作用到$ \xi $ 的动态中, 从而实现系统的降阶. 为方便起见, 下文中省略时间$ t $和状态$ x $在函数中的显式表达. 定义如下的新误差变量
$$ \begin{align} \begin{aligned} {\bar e_2} ={\dot {\bar e}_1} - {K_1}{\bar e_1} \end{aligned} \end{align} $$ (5) 对变量$ {\bar e_2} $求导并代入式(1)可得
$$ \begin{align} \begin{aligned} {\dot {\bar e}_2} = \bar h +K_1x_2 - F - Bu - d \end{aligned} \end{align} $$ (6) 其中函数$ \bar h $表达式为 $ \bar h = \ddot {\bar \rho} - {K_1}\dot {\bar \rho} $.
同样地, 针对误差$ \underline{e}_1(t) $ 定义如下新变量
$$ \begin{align} \underline{e}_2 = \dot{\underline{e}}_1 - K_1\underline{e}_1 \end{align} $$ (7) 对$ \underline{e}_2 $求导并代入式(1)得到
$$ \begin{align} \dot{\underline{e}}_2 = - \underline{h} -K_1x_2 + F + Bu + d \end{align} $$ (8) 其中函数$ \underline{h} = \ddot {\underline \rho} - {K_1}\dot {\underline \rho} $. 可以发现, 如果定义函数$ h = \ddot{\rho}-K_1\dot{\rho} $, 那么有关系$ \bar{h} = -\underline{h} = h $.
根据式(2) ~ (4), 误差变量$ \bar{e}_2 $和$ \underline{e}_2 $ 可以写为
$$ \begin{align} \bar{e}_2 = \bar{g} - \xi \end{align} $$ (9) $$ \begin{align} \underline{e}_2 = \xi - \underline{g} \end{align} $$ (10) 其中函数$ \bar{g},\;\underline{g} $表达式为 $ \bar{g} = -\underline{g} =\dot{\rho} - K_1\rho $. 可以看出, 式(9)和(10) 与式(2)和 (3)类似, 反映新变量$ \xi $与新的性能边界函数$ \bar{g},\;\underline{g} $之间的距离. 基于以上定义, 本文提出如下的性能驱动控制PDC方法
$$ \begin{split} u(t) =\;& {B^{ - 1}(x,\;t)}( {K_1}{x_2}(t) - F(x,\;t)\; + \\ &{K_2}\xi(t)- {K_2}{\mathop{\rm sgn}} \left( \xi(t) \right)g(\rho(t)) \;- \\ & {\mathop{\rm sgn}} \left( \xi(t) \right)D- {\mathop{\rm sgn}} \left( \xi(t) \right)\psi(t)) \end{split} $$ (11) 其中 $K_2 = {\rm diag}\{K_{2,\;i}\}\in \mathbf{R}^{m\times m}\;(i=1,\;\cdots, \, m)$是对角阵, $K_{2,\;i} $是常数; 向量$ D = [D_1,\;\cdots,\;D_m]^{\rm{T}} $; $ \psi(t) $是需要设计的函数, 如下述定理1的条件1)所示; $ {\rm sgn}(\cdot) $ 表示对角化的符号函数, 即对于任意向量$ y = [y_1, \cdots,\;y_m]^{\rm{T}} $, 有$ {\rm sgn}(y) = {\rm diag}\{{\rm sgn}(y_1),\,\cdots,\,{\rm sgn}(y_m)\} \in \mathbf{R}^{m \times m} $. 函数$ g(\rho(t)) $是显含性能函数$ \rho(t) $ 的函数, 表达式为
$$ \begin{align} g(\rho(t)) = \dot{\rho}(t) - K_1\rho(t) \end{align} $$ (12) 下面, 首先给出如下定理来分析性能驱动控制对预设性能的满足情况.
定理 1. 针对符合假设1 ~ 3 的非线性系统(1), 设计性能驱动控制律(11), 当系统满足以下条件时, 状态$ x_1 $满足$ \underline{\rho} \preceq {x_1} \preceq \bar \rho $.
1) 函数$ \psi \succeq - h + H,\;{\rm{ }}\forall t \ge 0 $, 其中$ H= [{H _1},\;\cdots, {H _m}]^{\rm{T}}\succ 0 $, $H_i\;(i = 1,\;2,\;\cdots,\;m)$为正常数;
2) 选取合适的矩阵$ {K_1} $满足$ g \succeq 0,\;\forall t \ge 0 $;
3) $ {\bar e_2}\left( 0 \right) \succeq 0,\; {\underline e_2}\left( 0 \right) \succeq 0 $, $ {\bar e_1}\left( 0 \right) \succeq 0,\; {\underline e_1}\left( 0 \right) \succeq 0 $.
证明. 首先观察变量$ \xi $的变化. 不失一般性, 认为变量初值$ \xi(0) $ 由非负和负两种元素组成. 假设非负元素有$ l $个, 那么负元素有$ m-l $ 个. 定义如下变量
$$\left\{ \begin{aligned} {\xi _s} =\;& {\left[ {{\xi _{{s_1}}},\; \cdots ,\;{\xi _{{s_i}}},\; \cdots ,\;{\xi _{{s_l}}}} \right]^{\rm{T}}} \in {{\mathop{\mathbf{R}}\nolimits} ^l},\;\\ &{s_i} \in \{ 1,\;2,\; \cdots ,\;m\},\;i = 1,\; \cdots ,\;l \\ {\xi _q}=\; & {\left[ {{\xi _{{q_1}}},\; \cdots ,\;{\xi _{{q_j}}},\; \cdots ,\;{\xi _{{q_{m - l}}}}} \right]^{\rm{T}}} \in {{\mathop{\mathbf{R}}\nolimits} ^{m - l}},\;\\ &{q_j} \in \{ 1,\;2,\; \cdots ,\;m\},\;j = 1,\; \cdots ,\;m - l \end{aligned}\right. $$ (13) 其中标量$ 0 \le l \le m $. 因此向量初值$ {\xi _s}\left( 0 \right) $和$ {\xi _q}\left( 0 \right) $ 分别满足$ {\xi _s}\left( 0 \right) \succeq 0,\;{\xi _q}\left( 0 \right) \prec 0 $.
可以看出, 通过定义$ \xi_s $和$ \xi_q $实现对变量$ \xi\in \mathbf{R}^{m} $在$ t=0 $时刻的划分. 对于系统(1), 存在时刻$ t_1 $使得这种划分在$ t\in[0,\; t_1] $保持不变. 因此引入划分矩阵$ I_s\in \mathbf{R}^{l\times m} $ 和$ I_q \in \mathbf{R}^{(m-l)\times m} $, 表示对任意$ m $ 维向量的划分, 其中矩阵$ I_s $ 和$ I_q $ 表示从单位矩阵$ I_{m \times m} $ 中挑出对应于$ s_i $ 和$ q_j $ 的行向量组成的矩阵. 比如, 对于任意三维向量$ y = [y_1\;\;y_2\;\;y_3]^{\rm{T}} $, 如果$ s_1=1,\;s_2=3 $ 和$ q_1 =2 $, 则表示向量 $ y $ 划分为两个向量$ y_s = [y_1\;\;y_3]^{\rm{T}} $ 和$ y_q = [y_2] $, 此时划分矩阵$ I_s = \begin{bmatrix} 1&0&0\\0&0&1 \end{bmatrix} $, $ I_q \;=\; \begin{bmatrix} 0&1&0\ \end{bmatrix} $. 对于划分矩阵$ I_s $和$ I_q $, 容易验证如下性质: $ y_s = I_s y,\;y_q = I_q y,\; y = I_s^{\rm{T}}y_s+I_q^{\rm{T}}y_q $. 本文后续所出现的所有划分矩阵$ I_s $ 和$ I_q $都是由变量$ \xi $中元素的正负性所唯一决定的. 因此, 状态变量$ \xi $, 函数$ g,\;\psi $ 和向量$ D $ 可重新表示为
$$ \begin{align} \xi = I_s^{\rm{T}}{\xi _s} + I_q^{\rm{T}}{\xi _q} \end{align} $$ (14) $$ \begin{align} g = I_s^{\rm{T}}{g_s} + I_q^{\rm{T}}{g_q} \end{align} $$ (15) $$ \begin{align} \psi = I_s^{\rm{T}}{\psi _s} + I_q^{\rm{T}}{\psi _q} \end{align} $$ (16) $$ \begin{align} D = I_s^{\rm{T}}{D_s} + I_q^{\rm{T}}{D_q} \end{align} $$ (17) 其中$ g_s \in \mathbf{R}^l,\;g_q\in \mathbf{R}^{m-l} $, $ \psi_s\in \mathbf{R}^l,\;\psi_q\in \mathbf{R}^{m-l} $ 和$ D_s \in \mathbf{R}^l,\;D_q\in \mathbf{R}^{m-l} $分别是基于$ s $和$ q $对函数$ g,\;\psi $ 和向量$ D $ 的划分.
考虑$ {\rm sgn}(\cdot) $函数的定义, 经计算可进一步得到
$$ \begin{align} {\mathop{\rm sgn}} \left( \xi \right)g = I_s^{\rm{T}}{g_s} - I_q^{\rm{T}}{g_q} \end{align} $$ (18) $$ \begin{align} {\mathop{\rm sgn}} \left( \xi \right)\psi = I_s^{\rm{T}}{\psi _s} - I_q^{\rm{T}}{\psi _q} \end{align} $$ (19) $$ \begin{align} {\mathop{\rm sgn}} \left( \xi \right)D = I_s^{\rm{T}}{D_s} - I_q^{\rm{T}}{D_q} \end{align} $$ (20) 针对控制律(11), 计算可得
$$ \begin{split} {I_s}Bu =\;& {I_s}{K_2}I_s^{\rm{T}}{\xi _s} + {I_s}{K_2}I_q^{\rm{T}}{\xi _q} - {I_s}{K_2}I_s^{\rm{T}}{g_s} \;+ \\ & {I_s}{K_2}I_q^{\rm{T}}{g_q} - {I_s}I_s^{\rm{T}}{\psi _s} + {I_s}I_q^{\rm{T}}{\psi _q} \;- \\ & {I_s}\left( {F-K_1x_2} \right) - {I_s}I_s^{\rm{T}}{D_s} + {I_s}I_q^{\rm{T}}{D_q}{\rm{ }} \end{split} $$ (21) $$ \begin{split} {I_q}Bu =\;& {I_q}{K_2}I_s^{\rm{T}}{\xi _s} + {I_q}{K_2}I_q^{\rm{T}}{\xi _q} - {I_q}{K_2}I_s^{\rm{T}}{g_s}\; + \\ & {I_q}{K_2}I_q^{\rm{T}}{g_q} - {I_q}I_s^{\rm{T}}{\psi _s} + {I_q}I_q^{\rm{T}}{\psi _q}\; - \\ & {I_q}\left( {F-K_1x_2} \right) - {I_q}I_s^{\rm{T}}{D_s} + {I_q}I_q^{\rm{T}}{D_q}{\rm{ }} \end{split} $$ (22) 注意有$ {I_s}I_q^{\rm{T}} = {{\boldsymbol{0}}_{l \times \left( {m - l} \right)}} $ 和$ {I_q}I_s^{\rm{T}} = {{\boldsymbol{0}}_{\left( {m - l} \right) \times l}} $, 因此利用式(16)得到
$$ \begin{align} {I_s}\psi = {I_s}I_s^{\rm{T}}{\psi _s} + {I_s}I_q^{\rm{T}}{\psi _q} = {I_s}I_s^{\rm{T}}{\psi _s} \end{align} $$ (23) $$ \begin{align} {I_q}\psi = {I_q}I_s^{\rm{T}}{\psi _s} + {I_q}I_q^{\rm{T}}{\psi _q} = {I_q}I_q^{\rm{T}}{\psi _s} \end{align} $$ (24) 由于$ K_2 $是对角矩阵, 可得$ {I_s}{K_2}I_q^{\rm{T}} = {{\boldsymbol{0}}_{l \times \left( {m - l} \right)}} $ 和$ {I_q}{K_2}I_s^{\rm{T}} = {{\boldsymbol{0}}_{\left( {m - l} \right) \times l}} $, 代入式(21) 和(22) 得
$$ \begin{align} {I_s}Bu = - {I_s}{K_2}I_s^{\rm{T}}{\bar e_{2,\;s}} - {I_s}\left( {F -K_1x_2+ \psi } \right) - {I_s}D \end{align} $$ (25) $$ \begin{align} {I_q}Bu = {I_q}{K_2}I_q^{\rm{T}}{\underline{e}_{2,\;q}} - {I_q}\left( {F -K_1x_2 - \psi } \right) + {I_q}D \end{align} $$ (26) 其中向量$ {\bar e_{2,\;s}} = {[ {{{\bar e}_{2,\;{s_1}}},\; \cdots ,\;{{\bar e}_{2,\;s_l}}} ]^{\rm{T}}} $ 和$ {\underline{e}_{2,\;q}} = [ {{\underline{e}}_{2,\;{q_1}}}, \cdots , \;{{\underline{e}}_{2,\;q_{m-l}}} ]^{\rm{T}} $ 定义为
$$ \begin{align} {\bar e_{2,\;s}} = {g_s} - {\xi _s} = {\bar g_s} - {\xi _s} \end{align} $$ (27) $$ \begin{align} {\underline{e}_{2,\;q}} = {\xi _q} + {{g}_q} ={\xi _q} - {\underline{g}_q} \end{align} $$ (28) 其中$ \bar{g}_s = -\underline{g}_s = g_s $, $ \bar{g}_q = -\underline{g}_q = g_q $, 根据条件2)可知$ g_s \succeq 0,\; g_q \succeq 0 $.
代入式(6), (8), (25) 和(26) 得到$ {\dot {\bar e}_{2,\;s}} $和$ {\dot {\underline e}_{2,\;q}} $, 并考虑$ \bar{h} = -\underline{h} = h $可得
$$ \begin{split} {{\dot {\bar e}}_{2,\;s}} = \;&{I_s}{{\dot {\bar e}}_2} = {I_s}{K_2}I_s^{\rm{T}}{{\bar e}_{2,\;s}} + {I_s}\psi\; + \\ &{I_s}h + {I_s}\left( {D - d} \right) \end{split} $$ (29) $$ \begin{split} {{\dot {\underline e}}_{2,\;q}} = \;&{I_q}{{\dot {\underline{e}}}_2} = {I_q}{K_2}I_q^{\rm{T}}{{\underline{e}}_{2,\;q}} + {I_q}\psi \;+\\ &{I_q}h + {I_q}\left( {D + d} \right) \end{split} $$ (30) 因此, 系统方程可表示为
$$ \begin{split} \left[ {\begin{array}{*{20}{l}} {{{\dot {\bar e}}_{2,\;s}}}\\ {{{\dot {\underline e}}_{2,\;q}}} \end{array}} \right] = \;&\left[ {\begin{array}{*{20}{c}} {{I_s}{K_2}I_s^{\rm{T}}}&{{{\boldsymbol{0}}_{l \times \left( {m - l} \right)}}}\\ {{{\boldsymbol{0}}_{\left( {m - l} \right) \times l}}}&{{I_q}{K_2}I_q^{\rm{T}}} \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {{{\bar e}_{2,\;s}}}\\ {{{{\underline e}}_{2,\;q}}} \end{array}} \right] +\\ &\left[ {\begin{array}{*{20}{l}} {{v_s}}\\ {{v_q}} \end{array}} \right] \\[-1pt]\end{split} $$ (31) 其中向量$ {v_s} $和向量$ {v_q} $表示为
$$ \begin{align} {v_s} = {I_s}\left( {\psi + h + D - d} \right) \end{align} $$ (32) $$ \begin{align} {v_q} = {I_q}\left( {\psi + h + D + d} \right) \end{align} $$ (33) 根据假设3可知$ {D_i} - |{d_i}|\ge0 $, 同时结合条件1)可得对$ \forall t \ge 0 $均满足
$$ \begin{align} {v_s} = {I_s}\left( {\psi + h + D - d} \right) \succ 0 \end{align} $$ (34) $$ \begin{align} {v_q} = {I_q}\left( {\psi + h + D + d} \right) \succ 0 \end{align} $$ (35) 因为$ {K_2} \in {{\mathop{\mathbf{R}}\nolimits} ^{m \times m}} $为对角阵, 且满足Metzler 矩阵性质, 由此矩阵$ \left[ {\begin{array}{*{20}{c}} {{I_s}{K_2}I_s^{\rm{T}}}&{{{\boldsymbol{0}}_{l \times \left( {m - l} \right)}}}\\ {{{\boldsymbol{0}}_{\left( {m - l} \right) \times l}}}&{{I_q}{K_2}I_q^{\rm{T}}} \end{array}} \right] \in \mathbf{R}^{m \times m} $也属于Metzler矩阵.
由条件3)可知$ {\bar e_2}\left( 0 \right) \succeq 0 $和$ {\underline e_2}\left( 0 \right) \succeq 0 $, 由引理1可得
$$ \begin{align} \left[ \begin{array}{l} {{\bar e}_{2,\;s}}\\ {{{\underline e}}_{2,\;q}} \end{array} \right] \succeq 0,{\rm{ \;\;}}\forall t \ge 0 \end{align} $$ (36) 注意上式仅保证$ \xi_s\preceq \bar{g}_s $和$ \xi_q \succeq \underline{g}_q $, 因此需要将其扩展到$ \underline{g} \preceq \xi \preceq \bar{g} $.
如果定义$ {{\underline e}_{2,\;s}} = \xi_s - \underline{g}_s $和$ \bar{e}_{2,\;q} = \bar{g}_q - \xi_q $, 那么有 $ {\bar e_{2,\;s}} - {{\underline e}_{2,\;s}} = - 2{\xi _s} $和$ {\bar{e}_{2,\;q}} - {{\underline e}_{2,\;q}} = - 2{\xi _q} $, 变换得到$ {\underline e_{2,\;s}} = {\bar e_{2,\;s}} + 2{\xi _s} $和$ {\bar e_{2,\;q}} = {\underline e_{2,\;q}} - 2{\xi _q} $. 注意根据$ s $和$ q $的划分可得$ \xi_s \succeq 0,\; \xi_q \prec 0 $, $ t\in [0,\; t_1] $, 其中时刻$ t_1 $ 表示$ s $ 和$ q $的划分发生改变的时刻, 此时首次$ \xi_s $有元素小于0 或$ \xi_q $有元素大于0. 那么, 由式(36)可得$ {\underline e_{2,\;s}}\succeq 0,\; {\bar e_{2,\;q}}\succeq 0 $, 进而有$ \bar{e}_2\succeq 0,\; \underline{e}_2 \succeq 0 $, 对于$ t\in [0,\; t_1] $. 在时刻$ t_1 $后, 假设状态变量$ \xi $按照$ \bar{s} $和$ \bar{q} $ 划分, 即非负向量和负向量被重新分组为$ {\xi _{\bar s}} \in {\mathop{\mathbf{R}}\nolimits} ^{\bar l} $, $ {\xi _{\bar q}} \in {{\mathop{\mathbf{R}}\nolimits} ^{m - \bar l}} $, 其中$ 0 \le \bar{l}\le m $. 此时依据与上面过程类似的推导, 可得$ \bar{e}_2\succeq 0,\; \underline{e}_2\succeq 0 $, 对于$ t\in (t_1,\; t_2] $, 其中$ t_2 $ 是下一次$ \bar{s} $ 和$ \bar{q} $ 划分发生改变的时刻. 以此类推, 我们可得
$$ \begin{align} \left[ \begin{array}{l} {{\bar e}_{2}}\\ {{{\underline e}}_{2}} \end{array} \right] \succeq 0,{\rm{ \;\;}}\forall t \ge 0 \end{align} $$ (37) 因此可得$ \underline{g} \preceq \xi \preceq \bar{g},\;\forall t \ge 0 $. 根据式(5)和(7)可得
$$ \begin{align} \left[ {\begin{array}{*{20}{l}} {{{\bar e}_{2}}}\\ {{{{\underline e}}_{2}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{l}} {{{\dot {\bar e}}_{1}}}\\ {{{\dot {\underline e}}_{1}}} \end{array}} \right] - \left[ {\begin{array}{*{20}{c}} {{K_1}}&{{{\boldsymbol{0}}_{m \times {m} }}}\\ {{{\boldsymbol{0}}_{ {m} \times m}}}&{{K_1}} \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {{{\bar e}_{1}}}\\ {{{{\underline e}}_{1}}} \end{array}} \right] \end{align} $$ (38) 因为$ K_1 $为对角阵, 且满足Metzler矩阵性质, 由此矩阵$ \left[ {\begin{array}{*{20}{c}} {{K_1}}&{{{\boldsymbol{0}}_{m \times {m} }}}\\ {{{\boldsymbol{0}}_{ {m} \times l}}}&{{K_1}} \end{array}} \right] \in {{\mathop{\mathbf{R}}\nolimits} ^{2m \times 2m}} $也属于Metzler 矩阵. 根据引理1可得, 当初值满足$ \bar e_1(0)\succeq 0, \underline e_1(0)\succeq0 $时, 下式成立
$$ \begin{align} \left[ {\begin{array}{*{20}{l}} {{{\bar e}_{1}}}\\ {{{{\underline e}}_{1}}} \end{array}} \right] \succeq 0,\;\forall t \ge 0 \end{align} $$ (39) 即$ \forall t \ge 0 $, 有$ {\bar e_{1}} = {\bar \rho} - {x_{1}} \succeq 0 $ 和$ {{\underline e}_{1}} = {x_{1}} - {{\underline \rho}} \succeq 0 $均成立, 那么$ {\underline \rho} \preceq {x_{1}} \preceq {\bar \rho} $, $ \forall t \ge 0 $ 成立.
□ 以上定理保证状态始终在预设的边界性能函数中, 下面分析整个系统的稳定性. 通过适当选择性能函数$ \rho $, 本文所提出的PDC控制律可以使系统具有不同稳定性. 考虑到PPC控制常用到如下性能函数
$$ \begin{align} {\rho _i}\left( t \right) = \left( {{\rho _{i,\;0}} - {\rho _{i,\;\infty }}} \right){{\rm e}^{ - {l_i}t}} + {\rho _{i,\;\infty }} \end{align} $$ (40) 其中$ {\rho _{i,\,0}} > \rho_{i,\,\infty}\ge0,\, l_i>0 $, 式(40)属于标准的幂次函数. $ \rho_{i,\,\infty} > 0 $ 时, 性能函数$ \rho $最后趋近$ \rho_{i,\;\infty} $; $ \rho_{i,\;\infty} =0 $时, 性能函数$ \rho $最后趋近零. 基于性能函数这一特点, 分析所提出的控制方法所实现的稳定性, 如下所述:
定理 2. 针对符合假设1 ~ 3的非线性系统(1), 采用提出的性能驱动控制律(11), 且满足定理1中的3个条件, 那么有如下结论:
1) $K_{2,\;i} < 0\;(i=1,\;\cdots,\;m)$时, 状态$ x_1 $的稳定性由如式(40)的性能函数$ \rho $ 决定, 即如果$ \rho_{i,\;\infty} > 0 $, 那么状态$ x_1 $ 是李雅普诺夫意义下稳定的; 如果$ \rho_{i,\;\infty} =0 $, 那么状态$ x_1 $ 是渐近稳定的.
2) $K_{2,\;i} > 0\;(i=1,\;\cdots,\;m)$时, 如果选取$ K_1 $使得 $ h_i\le0 $, 那么不论$ \rho_{i,\;\infty} >0 $或$ \rho_{i,\;\infty} =0 $, 状态$ x_1 $必然是渐近稳定的.
证明. 由于闭环系统包含符号函数, 是非光滑的, 因此需要首先分析闭环系统的解的存在性. 首先, 引入Filippov解的概念[22], 其允许在不连续点定义系统的解. 然后结合文献[23]中的命题3, 可分析得到系统存在一个绝对连续的Filippov 解. 同时, 由于$ F(x,\;t)、B(x,\;t) $连续且在状态$ x $上满足局部Lipschitz条件, 根据文献[24] 可知, 该解是唯一的.
注意$ K_1 $是对角阵, 那么根据定理1中的条件2), 可知$ g \succeq 0 $, 进一步结合函数$ g $的定义(12), 可得$ K_{1,\;i}\le \dot{\rho}_i(t)/\rho_i(t),\; \forall t\ge 0 $. 由式(40) 可知, $ K_{1,\;i} $ 应该取负值, 即$K_{1,\;i} < 0\;(i=1,\;\cdots,\;m)$. 基于以上分析, 下面将分别对定理2的结论1)和2)进行证明.
1) $K_{2,\;i} < 0\;(i=1,\;\cdots,\;m)$
如果边界性能函数$ \bar{\rho} = -\underline{\rho} = \rho $的终值$ \rho_{i,\;\infty}>0 $, 那么对于任意常数$ \delta $使得: 当$ |\rho(0)|\le\delta $ 时, 存在常数$ \gamma $使得函数$ \rho $ 满足$ |\rho|\le \gamma $. 考虑到定理1的结论$ {\underline \rho} \preceq {x_{1}} \preceq {\bar \rho} $, $ \forall t \ge 0 $ 成立, 那么也有: 当$ |x_1(0)|\le\delta $ 时, 存在常数$ \gamma $使得状态$ x $满足$ |x_1|\le \gamma $, 因此状态$ x_1 $是李雅普诺夫意义下稳定的. 同样, 根据渐近稳定性的定义, 也可得到: 如果性能函数$ \rho $的终值$ \rho_{i,\;\infty}= 0 $, 那么状态$ x_1 $ 是渐近稳定的.
2) $K_{2,\;i} > 0\;(i=1,\;\cdots,\;m)$
对变量$ \xi(t) $求导并代入控制律(11)可得
$$ \begin{align} \begin{aligned} \dot{\xi} = K_2\xi-K_2{\rm sgn}(\xi)g - {\rm sgn}(\xi)D - {\rm sgn }(\xi)\psi +d \end{aligned} \end{align} $$ (41) 考虑Lyapunov函数$V_{\xi} = (1/2) \xi^{\rm{T}}\xi$, 对其求导得
$$ \begin{split} \dot{V}_{\xi} =\;& \xi^{\rm{T}}K_2\xi-\xi^{\rm{T}}K_2{\rm sgn}(\xi)g - \xi^{\rm{T}}{\rm sgn}(\xi)D \;- \\ & \xi^{\rm{T}}{\rm sgn}(\xi)\psi + \xi^{\rm{T}}d \end{split} $$ (42) 注意到$ K_2 $是对角阵
$$ \begin{split} \dot{V}_{\xi}=& \sum_{i=1}^{m}K_{2,\;i}|\xi_i|(|\xi_i|-g_i) - \sum_{i=1}^{m}|\xi_i|(D_{i}-{\rm sgn}(\xi_i)d_i) \,- \\ & \sum_{i=1}^{m}|\xi_i|\psi_i = \sum_{i=1}^{m}K_{2,\;i}|\xi_i|(|\xi_i|-g_i) \;- \\ &\sum_{i=1}^{m}|\xi_i|(\varepsilon_i+|d_i|-{\rm sgn}(\xi_i)d_i)-\sum_{i=1}^{m}|\xi_i|\psi_i\\[-1pt] \end{split} $$ (43) 其中$\varepsilon_i={D_i} - |{d_i}|\ge0\,(i = 1,\,2,\,\cdots,\,m)$. 考虑到$ \underline{g} \preceq \xi \preceq \bar{g} $, 那么有$ |\xi_i|\le {g_i} $. 而且, $ |d_i|\ge {\rm sgn}(\xi_i)d_i $ 必然成立. 注意此时$ K_{2,\;i}>0 $, 因此
$$ \begin{split} \dot{V}_{\xi}\le\;& -\sum_{i=1}^{m}|\xi_i|\varepsilon_i-\sum_{i=1}^{m}|\xi_i|\psi_i \le -\sum_{i=1}^{m}|\xi_i|\psi_i \le \\ &\sum_{i=1}^{m}|\xi_i|(h_i-H_i) \le -\sum_{i=1}^{m}|\xi_i|H_i \le -H_{\min}\sqrt{2V_{\xi}} \end{split} $$ (44) 其中$ H_{\min}=\min\{H_1,\;\cdots,\;H_m\}>0 $. 因此, 由文献[25]可知, 状态$ \xi $ 是有限时间收敛到0. 考虑到$ \xi = {x_2}- {K_1}{x_1} = \dot{x}_1- {K_1}{x_1} $以及$ K_{1,\;i}<0\;(i=1,\;\cdots,\;m) $, 那么状态$ x_1 $是渐近收敛的.
□ 注 2. 如式(40)所示, 当$ \rho_{i,\;\infty}>0 $时, $ \rho_i $的运动与李雅普诺夫稳定中系统状态的运动形式是一致的[26]. 而当$ \rho_{i,\;\infty}=0 $时, $ \rho_i $的运动与渐近稳定中系统状态的运动形式是一致的. 因此从定理2的证明可知, 性能函数$ \rho $决定系统的稳定性, 通过选取不同的$ \rho $可实现不同的稳定性.
注 3. 传统预设性能控制和通道控制方法中采用系统状态与边界性能函数的商和$ \ln $函数等形式, 导致突发受扰时如果系统状态等于边界时必然引发系统发散. 而如式(2)、(3)、(9)、(10) 所示, PDC方法的设计是基于系统状态与边界性能函数的差, 因此当$ K_{2,\;i}<0 $时, 由于突发受扰导致系统状态穿越边界性能函数并不会引起系统发散, 且能保证系统状态最终返回到预设边界性能函数内. 限于篇幅, 给出简略分析. 假设$ \hat{t}_1 $ 时刻, 突发扰动导致$ \xi $穿越边界函数$ \bar g $, 即$ \bar e_{2,\;s}(\hat t_1)\preceq0 $. 由式(31) 可得$ \dot {\bar e}_{2,\;s}(t) \succeq I_s K_2 I^{\rm{T}}_s \bar e_{2,\;s}(t)+C $, 其中常数向量$ C\succ0 $ 为向量$ v_s $的下界. 进而可得
$$ \begin{split} \bar e_{2,\;s}(t)\succeq\;& {\rm e}^{I_s K_2 I^{\rm{T}}_s (t-\hat{t}_1)} \int_{\hat{t}_1}^{t} {\rm e}^{-I_s K_2 I^{\rm{T}}_s \tau}C {\rm{d}}\tau\;+\\ &{\rm e}^{I_s K_2 I^{\rm{T}}_s (t-\hat t_1)}\bar e_{2,\;s}(\hat t_1)=\sigma_1+\sigma_2\end{split} $$ 其中
$$ \begin{split} \sigma_1=\;&{\rm e}^{I_s K_2 I^{\rm{T}}_s (t-\hat t_1)}(\bar e_{2,\;s}(\hat t_1)\;+\\ &{\rm e}^{-I_s K_2 I^{\rm{T}}_s \hat t_1}{(I_s K_2 I^{\rm{T}}_s)}^{-1}C)\\ \sigma_2=\;&-{\rm e}^{-I_s K_2 I^{\rm{T}}_s \hat t_1}{(I_s K_2 I^{\rm{T}}_s)}^{-1}C \end{split} $$ 由于$ I_s\succeq0 $, $ K_2\preceq0 $, $ C\succ0 $, 因此$ \sigma_2\succeq0 $. 而由指数函数性质可知, 随着$ t $增大$ \sigma_1 $将趋近于零, 所以一定存在$ \hat t_2>\hat t_1 $, 使得$ \bar e_{2,\;s}(\hat t_2)\succeq0 $, 此时便又回到定理1所描述的情况. 类似上述分析可知PDC方法允许短时间内出现$ \bar e_1,\; \underline e_1,\; \bar e_2,\; \underline e_2\preceq0 $, 即系统状态穿越性能函数. 因此, 这也是PDC方法相比于传统预设性能和通道控制的优势之一, 后续的仿真结果也验证了这一结论.
注 4. 从定理2的结论2)部分证明可以看出, 本文提出的PDC方法在$ K_{2,\;i}>0 $时蜕化为一种时变增益滑模控制. 后文仿真部分的情况3, 控制输入$ u $出现的抖振现象也符合滑模控制的特点. 但本文方法在$ K_{2,\;i}<0 $时与滑模控制并无关系. 这说明本文方法既涵盖滑模控制, 又能保证边界性能函数的满足.
2.2 性能驱动控制中的“驱动”的含义解释
下面对本文所提出的性能驱动控制中的“驱动”的含义进行解释与讨论. PDC 方法控制的是误差$ \bar{e}_2,\;\underline{e}_2 $ 和$ \bar{e}_1,\; \underline{e}_1 $, 因此在控制状态$ x_1 $ 时, 通过保证 $ \bar{e}_2,\;\underline{e}_2 \succeq 0 $和$ \bar{e}_1,\; \underline{e}_1 \succeq 0 $ 使得系统被驱动压缩到零. 下面的图1更为形象地展示了这种性能驱动的含义.
控制目标是设计控制器将$ x_1 $约束在区间$ [-\rho,\;\rho] $并驱使其收敛到零. 图1中, $ \bar{g}(0) $和$ \underline{g}(0) $与横轴相交于$ \rho_0 $和$ -\rho_0 $, 而$ \bar{g}(\infty) $和$ \underline{g}(\infty) $与横轴相交于$ \rho_{\infty} $和$ -\rho_{\infty} $. 如果终值$ \rho_{\infty} $设置为零, 那么$ \bar{g}(\infty) $和$ \underline{g}(\infty) $最后重合在红色虚线$ x_2-K_1x_1 = 0 $. 首先选取$ A $点的初值$ (x_1,\;x_2) $进行分析. 首先, $ A $点一定在区域$ \bar{g}(0) $和$ \underline{g}(0) $之间, 意味着$ \bar{e}_{2}\succeq0 $和$ \underline{e}_{2}(0)\succeq0 $成立. 然后, 随着边界函数的减小, 区域逐渐被压缩, 而$ A $点保持在此区域内. $ A $点可能最终运动到$ B $点, 也可能随着$ \xi $的变号移动到$ A' $点. 但系统仍会从$ A' $点运动到$ B' $点, 即最终状态必然在$ \bar{g}(\infty) $和$ \underline{g}(\infty) $区域内. 因此, 整个过程中状态$ x_1 $始终在$ [-\rho,\;\rho] $内. 上图清晰地展现了状态是通过控制器被驱动压缩到区间内, 稳定性取决于边界性能函数. 传统PPC方法的稳定性依赖于经过非线性变换后采用的不同控制设计方法, 与选取的边界函数无直接关系, 而这也是本文方法与传统PPC的主要不同.
3. 数值仿真
考虑卫星刚体姿态模型[27]
$$ \left\{\begin{aligned} {{\dot \Omega}}\left( t \right) &= F_\Omega(t)+\omega(t)\\ {{\dot \omega}}\left( t \right) &= F_\omega\left( {t} \right) + B\left( t \right)u(t) + d(t) \end{aligned}\right. $$ (45) 其中$ \Omega( t ) = {[\phi( t ) {\rm{\;\; }}\theta( t ) {\rm{\;\; }}\psi( t ) ]^{\rm{T}}} $, $ \phi( t ) $、$ \theta( t ) $、$ \psi( t ) $表示滚动角、俯仰角和偏航角; $ \omega(t)=[\omega_x(t),\;\omega_y(t), $ $ \omega_z(t)]^{\rm{T}} $ 表示角速度, 其中$ \omega_x(t) $、$ \omega_y(t) $、$ \omega_z(t) $分别表示滚转、俯仰和偏航通道的角速度; 控制力矩$ u(t)={[{u_x}( t ){\rm{\; }}{u_y}( t ){\rm{ \;}}{u_z}( t )]^{\rm{T}}} $, 其中$ u_x(t)、u_y(t) $、$ u_z(t) $表示三通道的控制力矩; $ B $为控制系数矩阵; 向量$ d( t ) $表示受到的干扰; 非线性函数$ F_\Omega(t)=[{\omega_0\psi(t)} {\omega_0}{\rm{\;\; }}{-\omega_0\phi(t)}]^{\rm{T}} $, $F_\omega(t) = {[{(1/J_x)(J_y - J_z)}}$ ${\omega_y\omega_z}{\rm{\; \;}}(1/ J_y) (J_z- J_x)\omega_z\omega_x{\rm{\;\; }}(1/J_z)(J_x-J_y)\omega_x{{\omega_y}]^{\rm{T}}}$, ${\omega _0} = 0.001\, 2 {\rm rad /s}$为轨道角速度. 定义$ {x_1}( t )\; = $ $ \Omega( t ),\; {x_2}( t )\;= F_\Omega(t)+\omega(t) $, 则卫星姿态系统为
$$\left\{ \begin{aligned} {{\dot x}_1}\left( t \right) &= {x_2}\left( t \right)\\ {{\dot x}_2}\left( t \right) &= F\left( {t} \right) + Bu\left( t \right) + d \end{aligned}\right. $$ (46) 其中$F(t)=\frac{\partial F_\Omega}{\partial \Omega^{\rm{T}}}\frac{\partial \Omega}{\partial t}+F_\omega(t)$. 仿真中其他参数取值如下: 转动惯量矩阵$ J = B^{-1} = {\rm diag}\{18.73,\; 20.77,\;23.63\} $, 初始姿态角$\Omega (0) = [0.011\,2\;\,- 0.044\, 5\,\;0.023\;5]^{\rm{T}}\; {\rm rad}$, 初始姿态角速度$\omega(0) = [-0.041 \, 6\;\,0.048\,4\;\,-0.055\,6]^{\rm{T}} {\rm rad/s}$, 干扰
$$ d(t) = \left[ {\begin{array}{*{20}{c}} {{A_0}\left( {3\cos \left( {{\omega _0}t} \right) + 1} \right)}\\ {{A_0}\left( {3\cos \left( {{\omega _0}t} \right) + 1.5\sin \left( {{\omega _0}t} \right)} \right)}\\ {{A_0}\left( {3\sin \left( {{\omega _0}t} \right) + 1} \right)} \end{array}} \right] $$ 其中${A_0} = 1.5 \times {10^{- 5}}\;{\rm (N\cdot m)}$.
为验证提出的方法的有效性, 本文首先选择如下两种情况来分别验证有界收敛和渐近收敛.
情况1 (有界收敛). 性能函数$ \rho \left( t \right)\; =\; (4\; - 0.2){{\rm e}^{ - 2t}}\; +\; 0.2\;=\;\bar\rho \left( t \right)=-\underline\rho \left( t \right) $. $ {K_1} = {\rm{diag}}\{-5, \;-5, -5\} $和$ {K_2} = {\rm{diag}}\{-2.5,\;-2.5,\; -2.5\} $.
情况2 (渐近收敛). 性能函数$ \rho \left( t \right) \;=\; (4\; - 0.02){{\rm e}^{ - 2t}} $, $ {K_1} = {\rm{diag}}\{-5,\, -5,\,-5\} $, $ {K_2} = {\rm{diag}}\{-2.5, -2.5,\; -2.5\} $.
以上两种情况中, 根据定理1选择函数$ \psi \left( t \right) = - h\left( t \right) + [1\times10^{-5}{\rm{\;\; }}1\times10^{-5}{\rm{\;\; }}1\times10^{-5}]^{\rm{T}} $, 对$ \forall t \ge 0 $ 满足条件$ \psi \left( t \right) \succeq - h(t) + H $. 情况1和情况2下的仿真结果分别如图2和图3所示. 可以看出, 状态$ {x_1}\left( t \right) $, 也就是俯仰偏航和滚转三个姿态角被成功地约束到规定的区域$ [{\underline{\rho}}\left( t \right),\;\bar \rho \left( t \right)] $; 而且, 注意$ K_{2,\;i}<0 (i=1,\;2,\;3) $, 系统稳定性与所选取边界性能函数一致, 从而验证了所提出控制方案的有效性. 图2(b)和图3(b)中控制输入$ u(t) $效果平滑. 同时, 从图2(c)和图3(c)中可以看出, 变量$ \xi \left( t \right) $也在边界驱动下有界收敛或渐近收敛至零. 图2(d)和图3(d)给出系统相平面的变化, 其所展示出的系统稳定性也与上文中所描述的一致.
设置以下情况分析$K_{2,\;i} > 0\;(i=1,\;2,\;3)$ 的情形, 进而讨论验证定理2所阐述的参数对稳定性的影响.
情况3. 性能函数$ \rho ( t ) = (4 - 0.2){{\rm e}^{ - 2t}}= \bar\rho ( t )= -\underline\rho ( t ) $, $ {K_1} = {\rm{diag}}\{-5,\,-5,\,-5\} $, $ {K_2} = {\rm{diag}}\{2.5, \; 2.5, 2.5\} $.
仿真结果如图4所示. 图4(a)显示出三个通道的姿态角均渐近收敛到零, 符合定理 2中的结论. 通过图4(b)观察到, 控制信号$ u(t) $产生抖动, 原因在于矩阵$ K_2 $ 此时是非Hurwitz的, 而在情况1中则是符合Hurwitz条件的. 当矩阵$ K_2 $是非Hurwitz 时, 系统(46) 是发散的, 因此通过频繁的切换从而将系统驱动压缩到零. 这种现象与滑模控制有些类似, 从图4(d)上能更为清晰地看出. 在变量$ \xi(t) $收敛至零后, 状态$ x(t) $ 沿$ x_2(t)-K_1x_1(t)=0 $“滑动”至零.
除上述情况外, 为分析本文的性能驱动控制方法相对于现有传统控制方法的优势,下面以滚转通道为例, 将本文方法与通道控制[15]以及传统预设性能控制[28]进行对比, 以进一步验证本文方法的有效性.
情况4. 性能函数$ \rho \left( t \right) = (4 - 0.02){{\rm e}^{ - 2t}} + 0.02 $, $ {K_1} \;=\; {\rm{diag}}\{-5,\; -5,\; -5\}, $ $ {K_2}\; =\; {\rm{diag}}\{-2.5, \;\; -2.5, -2.5\} $与预设性能控制方法以及通道控制方法进行对比. 首先, 给出预设性能控制方法的控制律设计: $ u_{PPC}(t)\;=\;B^{-1}(x,\;t)(K_{p2}(x_2-x_{2c})+\dot{x}_{2c}\;- F(x,\;t)) $, 其中$ x_{2c} = K_{p1} \varepsilon'(t)-(\underline{\rho}(t)\dot{\bar{\rho}}(t)- \dot{\underline{\rho}}(t)\bar{\rho}(t)\;- x_1(\dot{\bar{\rho}}(t)- \dot{\underline{\rho}}(t)))/(\bar{\rho}(t)\;-\;\underline{\rho}(t))$, $ \varepsilon'(t)\;=\;\ln(\vartheta(t)/(1\;- \vartheta(t)))$, $\vartheta(t)=(x_1-\underline{\rho}(t))/(\bar{\rho}(t)-\underline{\rho}(t)) $, $ K_{p1}=K_{p2}= \rm{diag}\{-0.8,\;-0.8,\;-0.8\} $. 其次, 给出通道控制方法的控制律设计: $ u_{FC}(t)=-\alpha(\varphi(t)\Vert e(t)\Vert)e(t) $, 其中$ \alpha(s)=1/(1-s^2),\; \varphi(t)=1/g(t),\; e(t)=\xi(t) $. 除此之外, 为探究三种方法在突加扰动下的鲁棒性, 突出本文提出方法的优势, 将该情况滚转通道的干扰在$ t\in(3,\;3.5)\;{\mathrm{s}} $ 时设置为$ d_1(t)=0.25\cos(\omega_0t)+0.5 $.
由图5可知, 在突加扰动情况下, 由于预设性能控制和通道控制均采用的是误差与性能函数的比值, 因此导致其在触碰到性能函数边界时, 系统很快就会发散, 控制输入$ u_x $ 也发散. 此时, 性能驱动控制方法用系统状态与边界性能函数做差的优势便体现出来. 如图5(a)的子图所示, 由于$ K_{2,\;i}<0 $时是Hurwitz的, 所以系统在扰动结束后会趋向于收敛, 当系统状态重新返回到边界性能函数之内后, 便又回到定理1和定理2所证明的情况中.
4. 结束语
本文针对一类不确定非线性系统提出一种保证系统状态满足预设边界性能函数的新型性能驱动控制方法, 其优势在于通过设计控制器使得系统的稳定性取决于边界性能函数的选取, 而不改变控制设计. 同时, 卫星姿态系统的数值仿真例子验证了方法的有效性. 未来将尝试把提出方法扩展到控制矩阵非方阵的情况, 并进一步放宽不确定性的约束假设条件, 从而适用于更为一般性的不确定非线性系统.
-
图 2 乳腺DCE-MRI图像的MRF能量和后验概率
Fig. 2 MRF energy and posterior probability of breast DCE-MRI image
图 6 肿块形乳腺癌灶分割结果((a)医生手动分割; (b)结合二值水平集和形态学的分割方法; (c)基于边界的ACM; (d)基于局部高斯能量的ACM; (e)基于区域的ACM; (f)基于MRF的区域ACM; (g) U-net; (h) V-net; (i)本文方法)
Fig. 6 Tumor-shaped breast lesion segmentation results ((a) Expert manual segmentation; (b) Segmentation method based on binary level set and morphological operation; (c) Edge-based ACM; (d) ACM based on local Gaussian energy; (e) Region-based ACM; (f) Regional ACM based on MRF; (g) U-net; (h) V-net; (i) The method in this paper)
图 7 非肿块形乳腺癌灶分割结果((a)医生手动分割; (b)结合二值水平集和形态学的分割方法; (c)基于边界的ACM; (d)基于局部高斯能量的ACM; (e)基于区域的ACM; (f)基于MRF的区域ACM; (g) U-net; (h) V-net; (i)本文方法)
Fig. 7 Non tumor-shaped breast lesion segmentation results ((a) Expert manual segmentation; (b) Segmentation method based on binary level set and morphological operation; (c) Edge-based ACM; (d) ACM based on local gaussian energy; (e) Region-based ACM; (f) Regional ACM based on MRF; (g) U-net; (h) V-net; (i) The method in this paper)
表 1 V-net和U-net训练集结果
Table 1 Results of V-net and U-net in training sets
分割
算法$R_{\rm{TP}}$
(均值±标准差)$R_{\rm{FP}}$
(均值±标准差)$D_{\rm{S}}$
(均值±标准差)U-net 0.9706±0.1189 0.6673±0.6662 0.6252±0.1456 V-net 0.9408±0.0637 0.0732±0.1287 0.8836±0.0917 
下载: 导出CSV
表 2 50个临床样本分割结果
Table 2 Segmentation results of 50 clinical samples
分割算法 $R_{\rm{TP}}$ (均值±标准差)$R_{\rm{FP}}$ (均值±标准差)$D_{\rm{S}}$ (均值±标准差)结合二值水平集和形态学的分割方法 0.9759±0.0276 0.9100±1.6139 0.6910±0.2537 基于边界的ACM 0.9676±0.0551 0.6270±0.4238 0.5146±0.2272 基于局部高斯能量的ACM 0.9548±0.0818 0.8368±1.0217 0.6147±0.2040 基于区域的ACM 0.9892±0.0135 1.0604±1.0193 0.5857±0.2357 基于MRF的区域ACM 0.8709±0.1178 0.2018±0.3271 0.7524±0.1438 U-net 0.9389±0.1156 0.5988±0.3772 0.6143±0.1456 V-net 0.8762±0.1408 0.1131±0.0961 0.7920±0.1402 本文方法 0.9297±0.0413 0.0337±0.0222 0.8996±0.0410
下载: 导出CSV
表 3 算法的时间复杂度
Table 3 Time complexity of the algorithm
分割算法 时间复杂度 平均执行时间 (s) 结合二值水平集和形态学的分割方法 ${\rm O}(n)$ 3.3500 基于边界的ACM ${\rm O}(n)$ 3.2550 基于局部高斯能量的ACM ${\rm O}(n)$ 3.6790 基于区域的ACM ${\rm O}(n)$ 1.1188 基于MRF的区域ACM ${\rm O}(n^2)$ 2.7925 U-net ${\rm O}(n)$ 1.6102 V-net ${\rm O}(n)$ 3.3669 本文方法 ${\rm O}(n^2)$ 4.6575
下载: 导出CSV
-
[1] Fitzmaurice C, Akinyemiju T F, Al Lami F H, Alam T, Alizadeh-Navaei R, Allen C, et al. Global, regional, and national cancer incidence, mortality, years of life lost, years lived with disability, and disability-adjusted life-years for 29 cancer groups, 1990 to 2016: A systematic analysis for the global burden of disease study. JAMA Oncology, 2018, 4(11): 1553−1568 doi: 10.1001/jamaoncol.2018.2706 [2] Uematsu T, Kasami M, Watanabe J. Is evaluation of the presence of prepectoral edema on T2-weighted with fat-suppression 3 T breast MRI a simple and readily available noninvasive technique for estimation of prognosis in patients with breast cancer? Breast Cancer, 2014, 21(6): 684−692 [3] Cheon H, Kim H J, Lee S M, Cho S H, Shin K M, Kim G C, et al. Preoperative MRI features associated with lymphovascular invasion in node-negative invasive breast cancer: A propensity-matched analysis. Journal of Magnetic Resonance Imaging, 2017, 46(4): 1037−1044 doi: 10.1002/jmri.25710 [4] Fasihi M S, Mikhael W B. Overview of current biomedical image segmentation methods. In: Proceedings of the 2016 International Conference on Computational Science and Computational Intelligence. Las Vegas, USA: IEEE, 2016. 803−808 [5] van Dijk N P, Maute K, Langelaar M, van Keulen F. Level-set methods for structural topology optimization: A review. Structural and Multidisciplinary Optimization, 2013, 48(3): 437−472 doi: 10.1007/s00158-013-0912-y [6] 马超, 刘亚淑, 骆功宁, 王宽全. 基于级联随机森林与活动轮廓的3D MR图像分割. 自动化学报, 2019, 45(5): 1004−1014Ma Chao, Liu Ya-Shu, Luo Gong-Ning, Wang Kuan-Quan. Combining concatenated random forests and active contour for the 3D MR images segmentation. Acta Automatica Sinica, 2019, 45(5): 1004−1014 [7] 孙文燕, 董恩清, 曹祝楼, 郑强. 一种基于模糊主动轮廓的鲁棒局部分割方法. 自动化学报, 2017, 43(4): 611−621Sun Wen-Yan, Dong En-Qing, Cao Zhu-Lou, Zheng Qiang. A robust local segmentation method based on fuzzy-energy based active contour. Acta Automatica Sinica, 2017, 43(4): 611−621 [8] Zhang K H, Zhang L, Lam K M, Zhang D. A level set approach to image segmentation with intensity inhomogeneity. IEEE Transactions on Cybernetics, 2016, 46(2): 546−557 doi: 10.1109/TCYB.2015.2409119 [9] Rouhi R, Jafari M. Classification of benign and malignant breast tumors based on hybrid level set segmentation. Expert Systems with Applications, 2016, 46: 45−59 doi: 10.1016/j.eswa.2015.10.011 [10] 张旭梅, 范虹, 乔柱. 融合全局和局部信息的水平集乳腺MR图像分割. 计算机应用研究, 2015, 32(1): 307−311 doi: 10.3969/j.issn.1001-3695.2015.01.072Zhang Xu-Mei, Fan Hong, Qiao Zhu. Level set method combined global and local information and its application to breast MRI. Application Research of Computers, 2015, 32(1): 307−311 doi: 10.3969/j.issn.1001-3695.2015.01.072 [11] 虞红伟, 厉力华, 徐伟栋, 刘伟, 张娟, 邵国良. 一种基于水平集的多尺度乳腺肿块分割方法. 仪器仪表学报, 2010, 31(6): 1418−1423Yu Hong-Wei, Li Li-Hua, Xu Wei-Dong, Liu Wei, Zhang Juan, Shao Guo-Liang. Multi-scale approach for breast mass segmentation using level set method. Chinese Journal of Scientific Instrument, 2010, 31(6): 1418−1423 [12] 范虹, 朱艳春, 王芳梅, 张旭梅. 多分辨率水平集算法的乳腺MR图像分割. 物理学报, 2014, 63(11): Article No.118701Fan Hong, Zhu Yan-Chun, Wang Fang-Mei, Zhang Xu-Mei. Segmentation of breast MR images based on multiresolution level set algorithm. Acta Physica Sinica, 2014, 63(11): Article No.118701 [13] Kuo H C, Giger M L, Reiser I, Boone J M, Lindfors K K, Yang K, et al. Level set segmentation of breast masses in contrast-enhanced dedicated breast CT and evaluation of stopping criteria. Journal of Digital Imaging, 2014, 27(2): 237−247 doi: 10.1007/s10278-013-9652-1 [14] Rodtook A, Makhanov S S. Multi-feature gradient vector flow snakes for adaptive segmentation of the ultrasound images of breast cancer. Journal of Visual Communication and Image Representation, 2013, 24(8): 1414−1430 doi: 10.1016/j.jvcir.2013.09.009 [15] Xie Q Y, Pan X, Zhu H Q. A multiphase level set clustering approach using MRF-based Student's-t mixture model. In: Proceedings of the 8th International Conference on Wireless Communications and Signal Processing. Yangzhou, China: IEEE, 2016. 1−5 [16] Yang X, Gao X B, Tao D C, Li X L, Li J. An efficient MRF embedded level set method for image segmentation. IEEE Transactions on Image Processing, 2015, 24(1): 9−21 doi: 10.1109/TIP.2014.2372615 [17] Kallenberg M, Petersen K, Nielsen M, Ng A Y, Diao P F, Igel C, et al. Unsupervised deep learning applied to breast density segmentation and mammographic risk scoring. IEEE Transactions on Medical Imaging, 2016, 35(5): 1322−1331 doi: 10.1109/TMI.2016.2532122 [18] Rouhi R, Jafari M, Kasaei S, Keshavarzian P. Benign and malignant breast tumors classification based on region growing and CNN segmentation. Expert Systems with Applications, 2015, 42(3): 990−1002 doi: 10.1016/j.eswa.2014.09.020 [19] Ashraf A B, Gavenonis S, Daye D, Mies C, Feldman M, Rosen M, et al. A multichannel Markov random field approach for automated segmentation of breast cancer tumor in DCE-MRI data using kinetic observation model. In: Proceedings of the 14th International Conference on Medical Image Computing and Computer-Assisted Intervention. Toronto, Canada: Springer, 2011. 546−553 [20] Achuthan A, Rajeswari M, Ramachandram D, Aziz M E, Shuaib I L. Wavelet energy-guided level set-based active contour: A segmentation method to segment highly similar regions. Computers in Biology and Medicine, 2010, 40(7): 608−620 doi: 10.1016/j.compbiomed.2010.04.005 [21] Shahvaran Z, Kazemi K, Helfroush M S, Jafarian N, Noorizadeh N. Variational level set combined with Markov random field modeling for simultaneous intensity non-uniformity correction and segmentation of MR images. Journal of Neuroscience Methods, 2012, 209(2): 280−289 doi: 10.1016/j.jneumeth.2012.06.012 [22] 宋艳涛, 纪则轩, 孙权森. 基于图像片马尔科夫随机场的脑MR图像分割算法. 自动化学报, 2014, 40(8): 1754−1763Song Yan-Tao, Ji Ze-Xuan, Sun Quan-Sen. Brain MR image segmentation algorithm based on Markov random field with image patch. Acta Automatica Sinica, 2014, 40(8): 1754−1763 [23] Li B, Chen K, Peng G M, Guo Y X, Tian L F, Ou S X, et al. Segmentation of ground glass opacity pulmonary nodules using an integrated active contour model with wavelet energy-based adaptive local energy and posterior probability-based speed function. Materials Express, 2016, 6(4): 317−327 doi: 10.1166/mex.2016.1311 [24] 张喆, 韩德强, 杨艺. 基于证据马尔可夫随机场模型的图像分割. 控制与决策, 2017, 32(9): 1607−1613Zhang Zhe, Han De-Qiang, Yang Yi. Image segmentation based on evidential Markov random field model. Control and Decision, 2017, 32(9): 1607−1613 [25] Li C M, Kao C Y, Gore J C, Ding Z H. Minimization of region-scalable fitting energy for image segmentation. IEEE Transactions on Image Processing, 2008, 17(10): 1940−1949 doi: 10.1109/TIP.2008.2002304 [26] Khan S, Hussain M, Aboalsamh H, Mathkour H, Bebis G, Zakariah M. Optimized Gabor features for mass classification in mammography. Applied Soft Computing, 2016, 44: 267−280 doi: 10.1016/j.asoc.2016.04.012 [27] Kwedlo W. A new random approach for initialization of the multiple restart EM algorithm for Gaussian model-based clustering. Pattern Analysis and Applications, 2015, 18(4): 757−770 doi: 10.1007/s10044-014-0441-3 [28] 郑强, 董恩清. 一种新的基于二值水平集和形态学的局部分割方法. 电子与信息学报, 2012, 34(2): 375−381Zheng Qiang, Dong En-Qing. A new local segmentation method based on binary level set and morphological operation. Journal of Electronics and Information Technology, 2012, 34(2): 375−381 [29] Xu J, Monaco J P, Madabhushi A. Markov random field driven region-based active contour model (MaRACel): Application to medical image segmentation. In: Proceedings of the 13th International Conference on Medical Image Computing and Computer-Assisted Intervention. Beijing, China: Springer, 2010. 197−204 期刊类型引用(14)
1. 党红宇，唐利明，徐雅雅. 基于Bayesian最小误判的活动轮廓模型. 内江师范学院学报. 2024(04): 53-64 . 
百度学术2. 孙淑婷，刘铖枨，周广茵，韩锐，陈立超，羊月褀，许玥. 图像分割算法在医学图像中的应用综述. 现代仪器与医疗. 2024(02): 59-68 . 百度学术3. 蒲秋梅，殷帅，李正茂，赵丽娜. U型卷积网络在乳腺医学图像分割中的研究综述. 计算机科学与探索. 2024(06): 1383-1403 . 百度学术4. 刘侠，吕志伟，李博，王波，王狄. 基于多尺度残差双域注意力网络的乳腺动态对比度增强磁共振成像肿瘤分割方法. 电子与信息学报. 2023(05): 1774-1785 . 百度学术5. 阮旺，郝国生，王霞，胡晓婷，杨子豪. 面向目标识别的特征融合模糊模型及其应用. 计算机科学. 2023(S1): 505-511 . 百度学术6. 聂生东，魏传令，张小兵. 基于磁共振影像的乳腺癌计算机辅助检测方法综述. 上海理工大学学报. 2022(01): 1-10+41 . 百度学术7. 罗圣钦，陈金怡，李洪均. 基于注意力机制的多尺度残差UNet实现乳腺癌灶分割. 计算机应用. 2022(03): 818-824 . 百度学术8. 徐胜军，周盈希，孟月波，刘光辉，史亚. 基于多节点拓扑重叠测度高阶MRF模型的图像分割. 自动化学报. 2022(05): 1353-1369 . 本站查看9. 王小鹏，魏统艺，房超，朱生阳. 自适应非局部空间约束与K-L信息的模糊C-均值噪声图像分割算法. 控制理论与应用. 2022(07): 1261-1271 . 百度学术10. 王超，王永顺，狄凡. 快速自动模糊C-均值聚类彩色图像分割算法. 激光与光电子学进展. 2022(22): 82-88 . 百度学术11. 邱成健，刘青山，宋余庆，刘哲. 基于循环显著性校准网络的胰腺分割方法. 自动化学报. 2022(11): 2703-2717 . 本站查看12. 秦传波，宋子玉，曾军英，田联房，李芳. 联合多尺度和注意力-残差的深度监督乳腺癌分割. 光学精密工程. 2021(04): 877-895 . 百度学术13. 胡高珍，徐胜军，孟月波，刘光辉，冯峰，段中兴. 基于边缘约束局部区域MRF的图像分割方法. 计算机工程. 2021(06): 253-261+270 . 百度学术14. 刘壮盛，李昌林，衣利磊，李智，陈业航，李荣岗，罗学毛，龙晚生，冯宝. MRI纹理特征分析预测浸润性乳腺癌脉管浸润. 中国医学影像技术. 2020(11): 1637-1642 . 百度学术其他类型引用(8)
-
下载:
计量
- 文章访问数: 1564
- HTML全文浏览量: 322
- PDF下载量: 186
- 被引次数: 22
