一种自适应鲁棒最小体积高光谱解混算法

王天成 刘相振 董泽政 王海波

一种自适应鲁棒最小体积高光谱解混算法

王天成1刘相振1董泽政1王海波1

对高光谱图像解混的目的在于从低空间分辨率的高光谱图像中找到端元与对应的丰度.本文根据解混算法中的最小体积准则,提出了一种自适应鲁棒最小体积高光谱解混算法(Robust minimum volume based algorithm with automatically estimating regularization parameters for hyperspectral unmixing,RMVHU).本算法通过引入负数惩罚正则项,替换了同类算法中的丰度非负性约束(Non-negativity constraint,ANC),使算法对图像中的噪声与异常值具有更强的鲁棒性;采用循环最小化方法,将非凸优化问题分解为凸优化子问题,然后应用交替方向乘子法解决随着像素点个数增大带来的求解困难问题;对于正则项系数,本算法提出了一种自适应调整策略,提高了算法的收敛性,并且通过定性分析,说明了该调整方法的合理性.将算法应用于合成数据与实际数据,实验结果表明,与同类算法相比,本文提出的算法能够取得更为优秀的效果.

高光谱解混,交替方向乘子法,凸优化,最小体积,自适应估参

随着遥感传感器的飞速发展,高光谱图像由于包含丰富的光谱间信息与空间信息,在分类、探测等方面得到了更加广泛的应用.与普通图像不同,在高光谱图像中,每一个像素点(像元)包含了成百上千个波段的反射率信息[1−2].高光谱的高谱间分辨率,为像素级乃至亚像素级的图像分类与探测提供了条件,例如文献[3]就利用高光谱的空谱特性进行了有效的异常探测.在光谱图像中,若某一像元中只存在一种物质,我们将这样的像元称为纯像元[4].但是高光谱图像的空间分辨率较低,图像中存在着大量混合像元[4].混合像元广泛存在于高光谱图像中,是影响遥感分类精度和目标探测效果的重要因素之一.为了进一步利用高光谱数据,我们需要通过分解混合像元得到图像中一系列基本物质(这样的基本物质被称作端元)的光谱信息,同时还要求取端元在混合像元中的占比(丰度),而上述通过原始高光谱图像得到端元与丰度的过程就是高光谱图像的解混过程[5].

近年来,学者们提出了一系列基于线性混合模型的解混算法.在线性混合模型中,像元可以由一系列不相关的端元线性表示,且丰度满足“非负性”(ANC)与“和为1”(ASC)两个约束条件[6].在线性模型的基础上,基于凸面几何学的解混算法被广泛研究.该类算法的主要思想是:根据线性模型的全约束条件,所有数据点均被包含在某一类单形体中,在这些单形体中,由端元作为顶点构成的单形体体积是最小的.同时也有较多学者研究基于非负矩阵的解混算法[7−10]与稀疏解混算法[11−15].相较于其他算法,凸面几何学类解混算法运算速度快,解混精度较高,因此本文将继续深入研究凸面几何学的算法.

在基于凸面几何学的算法中,有一类基于纯像元假设的算法,如 VCA[16]、PPI[17]、NFINDR[18]等.虽然这些算法在图像满足纯像元假设时[5]有不俗的表现,但是若图像中存在大量高度混合的数据,解混精度将大大下降.为解决上述问题,学者们提出了一类基于最小体积变换的算法,该类算法能够在不满足纯像元假设的情况下取得较好的端元提取效果,如 MVES[19]、MVSA[20]、SISAL[21]等.但是由于MVES、MVSA这类算法须严格满足丰度非负约束,观测数据的异常值将会对解混精度造成很大的影响;而SISAL虽然通过惩罚项放宽了约束条件,但是在求解优化问题时,SISAL存在两个问题:1)将目标函数做了二次近似,将导致近似的目标函数由于高阶项的缺失,在某些情况下会大尺度地偏离原目标函数;2)二次项的Hesse矩阵采用对角矩阵近似而非目标函数的Hesse矩阵.上述问题的存在将给解混精度的提高带来一定的限制.

为了改善上述算法出现的问题,本文提出了一种自适应鲁棒最小体积解混算法,该算法有如下优点:

1)将约束条件放宽,通过引入负数惩罚正则项替代丰度非负性约束,获得了更强的抗干扰性能.

2)求解的目标函数根据行列式的形式展开,变换之后的目标函数与原函数完全等价,不再是二次近似,解决了二次逼近带来的误差.

3)为了能够解决惩罚正则项带来的大规模凸优化求解困难问题,应用交替方向乘子法[22](ADMM),获得了理想的计算精度与速度.

4)创造性地提出了一种自适应的正则系数调整方法,提高了算法的收敛性与稳定性,并通过定性的分析说明了此自适应参数调整方法的合理性.

5)对于ADMM中的惩罚系数,应用文献[22]中自适应调节的方法,进一步加快了算法的收敛速度.

文章结构主体如下所述.在第1节中详细介绍线性混合模型,凸面几何学的数学基础以及最小体积解混模型;在第2节中讲述自适应鲁棒最小体积高光谱解混算法(RMVHU)模型,参数的自适应调整方法,RMVHU算法步骤以及算法复杂度与收敛性;在第3节中讲述实验用的数据来源,并将数据应用于算法的实验与分析;在第4节中,总结实验的分析结果,得出结论.

1 几何解混相关工作

1.1 高光谱线性混合模型

通常情况下,在线性混合模型中,高光谱图像中的每个像元都可被近似认为是图像中各个端元的线性组合:

式中,p为端元个数,yyy∈RL为任意像元的L维光谱向量(L为图像波段数),为大小是L×p的端元矩阵为端元的向量表示.为丰度的向量表示,si表示像元中端元所占的比例,为误差项.

线性混合模型一般可分为3种情形:式(1)为无约束的线性混合模型;加上约束条件(2)则为非负约束混合模型;再加上约束条件(3)则为全约束混合模型.线性解混就是提取端元,同时求出各个端元在像元中所占的比例,得到端元丰度的过程.

在研究几何学高光谱解混算法之前,本文将依次介绍相关数学符号的定义、求解凸优化问题的著名算法框架——ADMM,以及仿射集与凸包的数学概念.

本文应用到的符合及其意义如表1所示:

表1 数学符号及其意义Table 1 Mathematical notations and their meaning

1.1.1 交替方向乘子法(ADMM)

考虑如下线性约束优化问题：

式中,µ为与收敛速度相关的常值,ddd为尺度对偶变量(Scaled dual variable)[22].

采用ADMM 的框架,能够在一步步的迭代过程中逼近凸问题的最优解.

1.1.2 仿射集与凸包

仿射包是一个仿射集,所以也能被表示为:

在文献'[23]中,通过“求解式(10)描述的问题,根据向量集与k,得到一组参数(C,d)来表示仿射集.

式中,约束CTC=Ik是为了满足rank(C)=k的秩约束条件,是投影到仿射集上的误差,定义为:

问题(10)有如下闭合解的形式[23]:

凸包中的点x若不能被表示为严格凸组合的形式,即: 若且θ/=eeei,∀i=1,···,p,称 x为凸包的顶点.若L=p−1且是仿射无关的,则称凸包为单形体,向量集为单形体顶点的集合.

仿射集、凸包的概念如下图所示:

图1 顶点个数为3时,凸包与仿射集的概念说明Fig.1 Illustration of convex hull and aきne hull when the number of vertices is three

从图1中可以看到,由三个顶点描述的仿射集是一个平面,而凸包则是在仿射集平面上,由顶点组成的三角形内的区域,这个三角形也被称为二维单形体.是

1.2 最小体积解混模型

一般来说,式(1)中的端元矩阵是列满秩的,即各个端元线性无关.又由于式(3)“和为1”的假设,文献[19]指出,由观测数据集组成的仿射包与由端元组成的仿射包是同一个.因此与都能由同一组参数(C,ddd)来表示,其中:

考虑如下优化问题:

式中,V(β1,···, βp) 代表由β1,···, βp作为顶点构成的单形体的体积,代表光谱数据的低维表示.文献[19]指出,在纯像元假设下,问题(18)的最优解如下:

同时文献[19]指出,在纯像元假设不严格成立的情况下,采用式(19)表示的可行解也接近问题(18)的最优解.根据文献[24],单形体的体积能够被表示为:

将式(20)与式(21)代入式(18),式(18)描述的优化问题能够被等价表示为:

式(23)即为最小体积解混算法模型的数学描述.

2 自适应鲁棒最小体积高光谱解混算法

2.1 RMVHU算法模型

在凸面几何学的最小体积类算法中,我们通常要求单形体能够包围所有的数据点集,即必须满足丰度系数非负性条件.然而,在实际情况下,在图像中很有可能存在异常点,或者由于噪声的影响导致像元分布在单行体之外.若在这些情况下仍强行满足丰度非负性条件,则可能出现为了将这些单形体之外的点包入估计的单形体中,所估计的单形体顶点与真实的顶点相差甚远的情况.因此,为了容忍异常点给端元估计带来的影响,我们需要构造一个负数惩罚正则项替换非负性约束,从而在尽可能保证单形体体积最小的情况下,增强算法对异常值、噪声值的鲁棒性.同时通过引入此正则项,减少在求解非凸优化问题时,陷入局部最小解的概率.

考虑引入负数惩罚正则项:式中,sij为矩阵S的元素.该正则项能够对负系数进行惩罚,而对于非负系数则没有任何影响.

式(24)中的h(sij),可以写成如下形式:

将此正则项代替优化问题(23)中的非负约束项,建立如下优化模型:

将式(25)代入式(26),上述问题等价为如下无约束非凸最优化问题:

式(27)描述的问题即为RMVHU算法的问题模型.但是式(27)中的目标函数非凸,很难求解得到全局最优解,因此,可以采用循环最小化思想,逐行更新矩阵变量H与的值,以将此问题转换为带有绝对值的凸优化问题,随后通过代数余子式展开行列式,且脱去|det(H)|项的绝对值,将凸优化问题转换为两个等价的凸优化子问题.

具体的,考虑更新矩阵变量H与的第i行值,得到RMVHU算法的凸优化子问题p∗与q∗:

式(31)中,Hi,j为H中元素hi,j的代数余子式.

2.2 RMVHU算法子问题的求解

问题(28)与(29)为带有l1范数的最优化问题.下面介绍一种求解此类的方法.

定理 1.问题(33)的最优解与问题(34)的最优解相同:

则目标函数值有:

只有当y0的第i个元素y0(i)=0或者z0的第i个元素z0(i)=0时,能够满足在当前可行解 x0下,目标函数取得最小值.根据此结论,结合约束条件y − z =A x + b ,有:

式中ai为矩阵A的行向量,bi为向量b的第i个元素.所以,问题(34)可以等价为问题(33)的形式.□

利用定理1,可以将凸优化子问题(28)与(29)转换为线性规划问题求解.

但是此种解法仅仅适用于图像像素点较少的情况,原因在于A ∈R2N×p,转换后新增的变量y∈R2N,z∈R2N,变量的个数为p+4N,约束条件的个数为1+6N.当像素点个数达到N=10000的时候,所需解决的就是一个中规模的线性规划问题.而这仅仅是算法中更新矩阵变量H与ggg第i行元素的子步骤,若要完成算法的一个完整的迭代过程,将带来很大的时间开销.因此需要能够快速求解子问题(28)与(29)的新方法.

2010年后,ADMM 在求解大型凸优化问题的领域展现了其强大的能力,也给求解此类凸优化问题带来了契机.

以子问题(28)为例,给出采用ADMM算法解决此问题的步骤.

问题(28)有如下等价形式:

此问题的ADMM形式为:

根据ADMM算法的框架,采用表2所描述的算法流程,依次更新x,z1, z2,d1,d2,求解子问题(28)中的p∗.

表2 RMVHU中求解子问题p∗的算法步骤Table 2 Steps for solving subproblem p∗in RMVHU

下面将介绍如何求解步骤2至步骤4的优化子问题.

步骤2的子问题为:

式(39)的问题为简单的无约束二次规划问题,目标函数梯度为0的极值点即为最优解.令问题(39)目标函数的梯度为0,得到步骤2的解析最优解:

步骤3的子问题为:

该问题是对变量z1解耦的不等式约束最优化问题,因此我们能够直接给出问题(41)的最优解形式.更方便的是,z1其实就是一个实数.式(41)最优解的解析形式表示为:

步骤4描述的问题为:

因此,问题的解可以采用著名的软阈值[25]给出:

式中,soft¡Ω,ν¢即为著名的软阈值函数.soft¡Ω,ν¢表示对矩阵Ω中每个元素Ωi,j作如下映射:

在矩阵的运算中,最为耗时的是求解矩阵的逆.我们可以发现,在问题(28)的求解中,除了采用式(40)更新 xk时需要进行矩阵的逆运算,其余变量的更新只要进行简单的矩阵加减乘除与元素大小的比较操作,因此着重关注式(40)中矩阵Ψ的大小.根据式(30)与(31),A∈R2N×p, c∈Rp,因此Ψ∈Rp×p,p为图像中的端元个数.在高光谱图像中,端元个数p几乎都在几十以内,因此,求解Ψ−1的运算量很小,并且在子问题的求解中,Ψ−1为常值,意味着在迭代过程中只需要计算一次Ψ−1,这意味着从时间与空间开销上考虑,求逆运算对子问题的求解并不会产生很大的影响.求解问题(29)的方法与表2中的步骤类似,下面直接给出应用表2中步骤求解时,问题(29)对应的步骤2至步骤4的最优解表达形式:

关于表2步骤7中的迭代停止条件,当满足如下两种情形中的一种时,停止迭代.

情形1.迭代步数k到达最大迭代步数.

情形2.原始残差与对偶残差(见下文第2.3.2节中定义)小于预设阈值.

特别要说明的是关于两个优化问题的优化值取舍的问题.如果p∗＜ q∗,则采用子问题(28)的解,否则采用子问题(29)的解.

最后,假设RMVHU算法最终求解问题(27)得到的最终解为H∗与ggg∗,则可由式(48)得到端元矩阵A:

式中,C与ddd由式(14)得到.

根据式(21),可得到第n个像元的丰度

式中,p为端元个数.

2.3 RMVHU算法参数的自适应调节策略

2.3.1 正则项系数的自适应调节策略

正则项系数的大小控制着正则项在整个优化函数中的权值,也影响着整个算法的解混效果,当正则项所占的比例较小时,其起的作用相应较小,反之亦然.在RMVHU算法中,若正则项系数λ不变,则在每一步的迭代中,由于 c Tx项每次都会变化,其所起的作用将会随每一次的更新而改变,不利于整体优化问题的收敛.因此有必要固定其在每一步子优化问题中的占比.受文献[26]启发,采用如下策略自适应调整参数λ:

式中, x0为 x在迭代之前的初值.ω为常值比例系数,实验中发现其值在30∼50之间算法效果较好.

下面通过定性的分析,表明采用上述的自适应参数调整法是行之有效的.

以上两点原因,解释了为何在正则项的系数与体积约束项大小成正比,与负数惩罚正则项大小成反比时,RMVHU算法能够取得很好的效果.

2.3.2 惩罚系数µ的自适应调整

根据文献[22],子问题ADMM 形式(38)中的惩罚项系数µ在如下简单自适应调节的策略下,使算法的收敛性得到提高:

式中,rk为原始残差,sk为对偶残差.根据原始残差与对偶残差在文献[22]中的定义,在RMVHU算法中,原始残差与对偶残差能够由以下公式进行计算:

算法的主要流程如下表所示:

表3RMVHU算法步骤Table 3 Procedure for solving RMVHU

在表3描述的算法流程中,在更新矩阵变量H与g 时,一次完整的迭代步骤为:采用步骤3至步骤5将矩阵变量的每一行元素都更新一遍.也就是说,在算法的整个求解过程中,分别最多求解M×(p−1)个子问题(28)与(29).算法在采用表2描述的ADMM算法求解子问题后是相对高效的,因此,RMVHU算法必然能够在有限时间内给出可行解.

2.4 RMVHU算法时间复杂度与收敛性分析

为进一步研究RMVHU算法的效率,本节将详细分析算法的复杂度.由于算法的本质为求解M×(p−1)个子问题(28)与(29),因此,将以求解子问题(28)的时间复杂度为分析的着眼点.

当采用表2所描述的算法求解子问题(28)时,更新xxx需做2Np+4N+p次加法或者乘法运算,更新z1需做p+2次加法或者乘法运算,更新 zz2需做2Np+12N 次加法或者乘法运算,更新d1需做p+2次加法或者乘法运算,更新dd2需做2Np+6N 次加法或者乘法运算,在计算时,需做次加法或者乘法运算,在计算时,需做2Np+4N+p+8次加法或者乘法运算.因此,在表2描述的算法中,循环一次需要做10Np+32N+5p+8次加法或者乘法运算,算法的时间复杂度为相应的为O(Np+N+p),可以发现,时间复杂度与图像像素点个数N 成正比,与端元个数p成正比,当p≪N 时,算法的时间复杂度能够被近似表示为O(N).

当采用定理1的方法,将问题(28)转化为线性规划问题求解时,变量个数为4N+p.采用Karmarkar算法求解时,时间复杂度[27]为O(N3.5).

显然,采用ADMM 算法求解在时间效率上具有显然的优势,尤其当N 很大的时候.

为了分析RMVHU算法的收敛性,本文将从如下两个方面进行阐述.

1)在式(27)描述的优化目标函数中,当固定正则项系数λ不变,采用循环最小化的方法,逐行更新H与g时,所求解的两个子问题都是凸问题,并且ADMM算法在求解凸问题时具有收敛性,因此,式(27)的目标函数值必将逐渐减少,说明在采用循环最小化的方法,执行步骤4至步骤5时,算法是收敛的.

2)当执行完步骤4与步骤5,算法仍未收敛时,将根据式(50)的正则系数自适应调整策略,重新计算λ,这时,根据第2.3.1节中的定性分析,算法将在此种策略下加速收敛.

3 仿真校验

本部分由两块内容组成,首先介绍实验数据的来源以及算法评价准则,随后给出实验结果并进行分析.

3.1 算法实验数据来源与评价准则

算法的实验数据分为合成光谱数据与真实光谱数据.

3.1.1 合成高光谱数据

在USGS光谱库[28]中选择p个光谱向量作为端元.图2显示了库中的3条光谱曲线.

图2 USGS库中不同物质的光谱曲线Fig.2 USGS library spectra of diあerent materials

根据文献[16]的方法,产生像元的丰度信息.像元的丰度信息满足狄利克雷分布,狄利克雷分布表示为:

式中,第i个端元的丰度 αi满足0≤ αi≤ 1,且丰度期望Γ(·)表示Gamma函数.狄利克雷的密度分布不仅能够保证丰度符合非负与和为1约束,还能够根据不同的参数µi,产生不同分布特征的丰度系数.为了产生高度混合的数据,对于丰度αi＞ρ的像元,将其丰度重新匹配,改为1/p.ρ为表示像元混合度的参数且满足ρ∈(0,1),ρ越小,图像混合度越高.

为了使数据更具有真实性,加入独立同分布的零均值高斯加性噪声.信噪比SNR表示为:

式中，y[n]为真实的光谱数据,σ2代表高斯噪声的方差,L为波段维数,N为像素点个数.

在仿真数据中,还应考虑异常点对算法的影响,因此需要生成异常点.假设异常点的个数为k,则重新生成第1至k个像素点的丰度.第i(0≤i≤k)个点的丰度由以下公式重新产生:

式中,δ为常数,k=randperm(p),代表在1与p之间满足均匀分布的随机整数,ss(i)为第i个像元丰度的向量表示,sm(i)代表ss(i)的第m个元素.

3.1.2 实际数据的来源

本文使用的实际数据是美国内华达州Cuprite地区1997年成像的224波段0.4µm 至2.5µm,大小为250×191的子图[29].这组数据在高光谱图像端元提取的研究中广泛应用,具有很好的代表性.由于水汽吸收的干扰和低信噪比的原因,第1∼2,104∼113,148∼167以及221∼224波段的数据被剔除.伪彩色子图如图3所示.

该地区的矿物分布图由Tricorder 3.3软件[30]提供,如图4所示.需要注意的是此物质分布在1995年就已问世,而Cuprite地区的数据是在1997年采集.因此该物质分布图仅作为解混效果好坏的参考.

图3 Cuprite地区高光谱图像伪彩色子图,R:2.109µm,G:2.209µm,B:2.308µmFig.3 The Pseudo color subimage of the AVIRIS Cuprite dataset,R:2.109µm,G:2.209µm,B:2.308µm

图4 Cuprite地区矿物分布图Fig.4 The mineral distribution map of Cuprite

3.1.3 算法评价准则

为了描述解混算法的效果好坏,通常需要研究解混后端元与真实端元的相似度以及解混后估计丰度与真实丰度的相似程度.可以定义光谱角距离(SAD)和误差平方根(RMSE)来评估实验所得的端元和丰度.

对于第i个端元,假设实际的端元的向量表示形式为 ai,采用解混算法得到的端元估计值的向量表示为第i个端元的光谱角距离表示为φi:

对于第i个端元,假设端元对应的真实丰度表示为si,采用解混算法得到与端元对应的丰度估计值表示为N 为图像像素点个数,第i个端元的均方根误差表示为θi:

为全局地评价解混算法的好坏,将采用如下两个平均指标:平均光谱角距离εφ与平均误差平方根εθ进行算法效果的衡量:

式(59)与(60)中,p为图像中端元个数.

3.2 实验及结果分析

实验分为合成数据解混与真实数据解混两组实验.在采用合成数据的实验时,首先研究τ与γ的取值对ADMM 算法收敛的影响,研究第2.3.1节中正则系数自适应调整的策略的有效性,随后在端元数p=3且存在异常点的数据集中,将RMVHU算法与MVES算法进行比较,阐述RMVHU算法的优点;其次将RMVHU算法与被广泛使用的VCA、MVES、MVSA、SISAL四种凸面几何学算法进行比较,分别研究真实端元个数,噪声大小,数据混合度ρ,异常点个数,像元个数以及错估端元个数时对不同算法解混效果的影响.最后,对于真实数据,给出RMVHU算法的解混结果.

3.2.1 合成数据的仿真

实验1.RMVHU惩罚系数的调整实验

本实验将进一步研究第2.3.2节中惩罚系数γ与τ的取值对算法收敛性的影响.实验中,γ取值为1,10,100,200,τ取值为0.5,1.5,2,4.

图5与图6为求解子问题(28)时,目标函数值在不同γ与τ下的变化趋势.从图5可以发现,当γ较小时,函数收敛较快,但波动较大,当γ增大时,目标函数虽然波动较小,但收敛速度变慢.产生上述结果的原因在于,γ较小时,只要原始残差与对偶残差存在差别,那么惩罚系数µ便会频繁的变化,虽然这样的变化能使收敛速度提高,但会带来收敛稳定性的不足;而当γ较大时,原始残差与对偶残差的差异性将很难影响到惩罚系数µ的改变,而合适的µ将对收敛速度产生较大的影响,因此,虽然ADMM算法具有收敛性,但是收敛速度将放缓.在图5中可以发现,γ=10时,函数收敛较快,波动较小.

图5 γ变化时目标函数值Fig.5 The values of object function when γ changes

图6 τ变化时目标函数值Fig.6 The values of object function when τ changes

从图6可以发现,当τ＜1时,函数不收敛,随着τ的增大时,函数收敛速度加快,将很快收敛于极值,但是波动将会增加,其原因在于,当对偶残差与原始残差较大时,若τ＜1,根据调整策略,惩罚系数µ将增大,这将导致对偶残差更大,势必造成算法的不收敛;当τ＞1,则根据调整策略,惩罚系数µ将减小,这将缩小对偶残差与原始残差的大小,使得算法收敛性得到提升,但是过大的τ将使对偶残差波动较大,带来收敛稳定性的不足.在图6中可以发现τ=2时,函数收敛较快,波动较小.

因此在接下来的实验中,取γ=10,τ=2.

实验2.正则项系数的自适应调节对算法的影响

为了进一步研究第2.3.1节的正则项系数的自适应调节策略对算法的影响,将进行相关实验进行对比验证.实验中将比较随着端元个数的变化时,两种方法的解混精度.固定系数的算法取参数λ=0.002,采用第2.3.1节中描述的策略的实验,取参数ω=40.

实验中,模拟数据的像素点个数N=10000,波段维数L=220,端元数p从3至8变化,SNR=30dB,数据混合度ρ=0.8,异常点个数为25,估计端元个数无错估.实验结果如图7与图8所示.

图7 正则项系数固定或自适应变化时,各算法平均光谱角距离Fig.7 The average spectral angle distance of diあerent algorithms when the regularization parameter is fi xed or changes automatically

图8 正则项系数固定或自适应变化时,各算法平均均方根误差Fig.8 The average root mean square error of diあerent algorithms when the regularization parameter is fi xed or changes automatically

图7与图8表明,当正则项系数固定时,实验效果将会随着端元个数的改变发生很大的变化.而采用第2.3.1节策略的算法则能够稳定地估计真实端元与丰度.由此可以说明,采用第2.3.1节策略将使算法具有更强的鲁棒性.

实验3.RMVHU的鲁棒性验证实验

本实验主要通过与一种流行的凸面几何学解混算法MVES进行比较,直观地表明RMVHU算法在存在异常点情况下的鲁棒性.

模拟数据的像素点个数N=10000,波段维数L=220,端元个数p=3(选取的端元光谱曲线见图2),SNR=30dB,数据混合度ρ=0.8,异常点个数为25(像素点总数的0.25%),估计端元个数为3(端元个数无错估).

RMVHU与MVES端元估计结果如图9所示.在图9中,实线代表真实的端元光谱曲线,虚线代表解混算法提取出的端元曲线.可以发现,在少量异常点的干扰下,RMVHU算法的端元结果远远优于MVES算法.其主要原因在于,RMVHU算法引入了正则项,该正则项能够容忍丰度值为负的情形,使得噪声与异常点对最小体积逼近的影响不那么明显,而MVES算法需要严格满足非负的条件,因此MVES计算的单形体顶点需严格包含所有观测数据,使异常点对最小体积逼近的影响尤为突出,导致了端元提取效果的剧烈下降.

图10显示了数据在二维空间的分布以及RMVHU、MVES算法提取的端元在二维空间的分布,此图亦印证了上述缘由.

表4与表5分别记录了RMVHU与MVES算法的光谱角距离,平均光谱角距离与均方根误差,平均均方根误差值.结果最优值在表中加粗.从表4与表5中可以发现,若图像数据中存在异常值,则在进行图像的解混时,RMVHU算法无论是在端元提取精度上,还是丰度反演精度上,都远远优于MVES算法.随后进行的合成数据仿真实验中,将分别研究真实端元个数,噪声大小,数据混合度,有无异常点,像素点个数以及错估端元个数对不同解混算法的影响.

表4 RMVHU与MVES算法的光谱角距离与平均光谱角距离Table 4 Spectral angle distance and average spectral angle distance via RMVHU and MVES

表5 RMVHU与MVES算法的均方根误差与平均均方根误差Table 5 Root mean square error and average root mean square error via RMVHU and MVES

图9 RMVHU与MVES端元估计结果,(a)、(b)、(c)为RMVHU的端元估计结果,(d)、(e)、(f)为MVES的端元估计结果Fig.9 Estimated results of endmembers via two algorithms:RMVHU and MVES((a)、(b)、(c)are the results of RMVHU and(d)、(e)、(f)are the results of MVES)

图10 数据集与端元在p−1(二)维子空间的投影Fig.10 Two dimensional subspace projection of datasets,estimated endmembers and real endmembers

实验4.真实端元个数变化的实验

在每一幅高光谱图像中,端元个数会根据不同的拍摄场景而变化,因此,研究端元个数对解混结果的影响很有必要.

本实验中,生成端元个数p从3∼10变化的8组数据,其他参数设置为:模拟数据的像素点个数N=10000,波段维数L=220,SNR=30dB,数据混合度ρ=0.8,异常点个数为25(像素点总数的0.25%),端元个数无错估.

平均光谱角距离结果如图11所示.平均均方根误差结果如图12所示.

图11与图12表明,相较于其他算法,端元个数的变化对RMVHU算法的影响并不很明显.RMVHU算法的解混效果相对于基于纯像元假设的算法—VCA得到的结果优秀很多.相比于需要严格满足非负性条件的最小体积变化算法:MVES与MVSA,RMVHU提取的端元的精度要优秀很多.对于同样可以不满足非负条件,但是采用二次逼近体积项的SISAL算法来说,由于RMVHU算法子问题中的体积项采用等价的行列式展开公式,因此其求解的精度也在大多数情况下比SISAL算法高.

图11 端元个数变化时,各算法平均光谱角距离Fig.11 The average spectral angle distance of diあerent algorithms when the number of endmembers changes

图12 端元个数变化时,各算法平均均方根误差Fig.12 The average root mean square error of diあerent algorithms when the number of endmembers changes

实验5.抗噪能力的测试实验

设计本实验的目的在于测试在不同大小的噪声下,RMVHU算法的稳定性.

本实验中,信噪比SNR 分别设置为15dB、20dB、25dB、30dB、35dB、40dB,获得6组数据,其他参数设置为:模拟数据的像素点个数N=10000,波段维数L=220,端元个数p=6,数据混合度ρ=0.8,异常点个数为25(像素点总数的0.25%),端元个数无错估.

在不同噪声大小下,平均光谱角距离与平均均方根误差结果如图13与图14所示.

从图13与图14中可以发现,信噪比越高,RMVHU算法的解混效果越好.由于VCA中采用了根据不同噪声大小进行不同的数据降维方法,因此VCA在低信噪比下有较突出的表现.而本文提出的RMVHU算法效果在低信噪比下的端元提取精度仅次于VCA,侧面印证了RMVHU算法正则项及自适应的正则算子调节方法有较强的抗噪能力.相较于其他几何学解混算法,在信噪比大于20dB的情况下,RMVHU算法的表现出色,能够得到最低的平均光谱角距离与平均均方根误差值.RMVHU算法较强的抗噪声能力得到了实验结果的支持.

图13 在不同噪声大小下,各算法平均光谱角距离Fig.13 The average spectral angle distance of diあerent algorithms when SNR changes

图14 在不同噪声大小下,各算法平均均方根误差Fig.14 The average root mean square error of diあerent algorithms when SNR changes

实验6.数据混合度变化的实验

设计本实验的目的在于:验证在不同混合度的场景下,RMVHU算法也有较为突出的表现,即能够胜任不同混合程度的高光谱图像的解混.

在本实验中,数据混合度ρ分别设置为0.65、0.7、0.75、0.8、0.85、0.9、0.95、1,获得 8组数据.其他参数设置如下:模拟数据的像素点个数N=10000,波段维数L=220,真实端元个数p为6,SNR为30dB,异常点个数为25(像素点总数的0.25%),端元个数无错估.

在不同混合度下,平均光谱角距离与平均均方根误差结果如图15与图16所示.

由图15与图16发现,RMVHU在不同的混合度下,都能取得最为优异的结果.同时,由于异常值的存在,使得本应该在混合度较低,即ρ接近1的情况下能够获得更好效果的纯像元提取类算法VCA变得精度很差.但是本文提出的解混算法RMVHU在异常点存在的情况下,对数据的混合度大小并不敏感.同时由于实验4中所述的原因,RMVHU相比于SISAL也更加优秀.本实验的成功意味着该算法在对不同混合度的高光谱图像解混时,能够取得很好的效果.

图15 在不同混合度下,各算法平均光谱角距离Fig.15 The average spectral angle distance of diあerent algorithms when the purity of pixels changes

图16 在不同混合度下,各算法平均均方根误差Fig.16 The average root mean square error of diあerent algorithms when the purity of pixels changes

实验7.高于异常点个数变化的实验

设计本实验的目的在于:验证在异常点个数发生变化时,RMVHU算法有较强的稳定性,并且与其他算法一起,比较异常点的有无及数量对算法解混精度的影响.

实验中, 异常点个数分别设置为像素点总数的 0% (无异常数据),0.25%、0.5%、0.75%、1%、1.25%、1.5%、1.75%、2%,获得9组数据.其他参数设置如下:模拟数据的像素点个数N=10000,波段维数L=220,真实端元个数p为6,SNR为30dB,数据混合度ρ设置为0.8,端元个数无错估.

在异常点个数改变的情形下,平均光谱角距离与平均均方根误差结果如图17与图18所示.

图17 在不同异常点个数下,各算法平均光谱角距离Fig.17 The average spectral angle distance of diあerent algorithms when the number of outliers changes

图18 在不同异常点个数下,各算法平均均方根误差Fig.18 The average root mean square error of diあerent algorithms when the number of outliers changes

由图17与图18可以发现,图像中无异常数据时,基于最小体积变化的算法RMVHU、MVES、MVSA及SISAL算法都能够取得较好端元提取与解混的效果;但是当出现异常数据时,除了RMVHU算法能够取得最为优异的解混结果,SISAL算法勉强能够较为准确地获得图像端元与丰度信息外,其他算法获得的结果都与真实值有较大偏差.同时,相比于SISAL算法的结果随着异常点个数增加而变差,RMVHU算法的解混误差基本不变,维持在最低的水平.本实验进一步说明了RMVHU算法的稳定性.

实验8.关于像素点个数的变化的实验

设计本实验的目的在于研究图像像素点个数对RMVHU算法解混精度的影响以及算法的时间复杂度(即能否在较短的时间内得到较为精确的解混结果).

在本实验中, 设置像素点个数为2000、4000、6000、8000、10000、12000、14000、16000、18000、20000,共得到10组数据.其他参数设置如下:真实端元个数p为6,波段维数L=220,SNR为30dB,数据混合度ρ设置为0.8,异常点个数为像素点总数的0.25%,端元个数无错估.计算机为联想S20工作站.

在不同大小的高光谱图像中,平均光谱角距离与平均均方根误差结果如图19与图20所示,运行时间如图21所示.

图19 在不同像素点个数下,各算法平均光谱角距离Fig.19 The average spectral angle distance of diあerent algorithms when the number of pixels changes

由图19与图20可以发现,相较于其他算法,RMVHU算法解混偏差并不会像SISAL、MVES、MVSA、VCA 等算法随着像素点的个数变化出现巨大的波动,它的表现一直很稳定.当像素点个数超过4000时,RMVHU算法的解混结果是所有算法中最优秀的.此实验验证了本文提出的RMVHU算法能够在各种大小的高光谱图像解混中,取得稳定且优秀的解混效果的说法.

在图21中,可以发现随着像素点个数的增多,RMVHU算法运行的时间线性增长,这也与第2.4节中算法的复杂度分析相吻合.在像素点个数较大时,由于不用求解大规模的线性规划问题,其时间消耗反而要比MVES算法少.所以,在有限时间内,应用RMVHU算法,定能得到理想的解混效果.相较于SISAL与MVSA等算法,虽然其时间消耗较多,但若更注重解混的精度与结果的稳定性,那么在时间上的一些牺牲是值得的.

图20 在不同像素点个数下,各算法平均均方根误差Fig.20 The average root mean square error of diあerent algorithms when the number of pixels changes

图21 在不同像素点个数下,各算法运行时间Fig.21 Computational time of diあerent algorithms when the number of pixels changes

实验6.关于错误估计端元个数的实验

在解混算法的实际应用中,大多数解混算法需要预先估计端元的个数,RMVHU算法也不例外.但是可能由于混合模型的非线性偏差,端元个数的估计有可能发生错误,因此,设计实验研究算法在错误估计端元个数时的表现也很有必要.

在本实验中,真实端元的个数p在3∼8间变化.估计的端元个数为p−1,模拟端元个数错误估计的情况.其他实验参数设置如下:模拟数据的像素点个数N=10000,波段维数L=220,SNR为30dB,数据混合度ρ设置为0.8,异常点个数设置为0.

在对被错估端元个数的高光谱图像解混后,平均光谱角距离与平均均方根误差结果如图22与图23所示.

图22 在错估端元个数时,各算法平均光谱角距离Fig.22 The average spectral angle distance of diあerent algorithms when the number of endmembers is incorrectly estimated

由图22发现在绝大多数端元个数错估的情况下,RMVHU算法相比其他的最小体积变换算法:MVSA、MVES及SISAL,能够取得更为优秀的端元估计结果.

在对丰度值的估计方面,由图23可以发现,RMVHU算法在不同的端元个数下都获得了不错的效果,且随着真实端元个数的增加,精度也有提高的趋势.对RMVHU算法而言,端元提取误差与丰度反演的误差随着真实端元个数的增加,也有降低的趋势,这也说明了RMVHU算法在错误估计端元个数的情况下仍然具有较强的鲁棒性.

3.2.2 真实数据的实验

采用Cuprite子图,剔除干扰波段后,文献[8]的研究结果表明,这片区域中有14种矿物.由于有些同种类不同化学成分的矿物的光谱只存在很微小的差别,因此在解混时,将端元的个数减少至11.解混后各物质丰度图如图24所示.

图23 在错估端元个数时,各算法平均均方根误差Fig.23 The average root mean square error of diあerent algorithms when the number of endmembers is incorrectly estimated

表6 RMVHU、VCA、MVES、MVSA、SISAL算法提取端元与真实端元的光谱角距离Table 6 The spectral angle distance of real datasets via unmixing algorithms:RMVHU,VCA,MVES,MVSA and SISAL

图24 Cuprite地区估计丰度图Fig.24 Estimated abundance maps of AVIRIS Cuprite

在图24的Cuprite地区估计丰度图中,(a)为Desert Varnish的丰度图,(b)为Sphene的丰度图,(c)为Nontronite的丰度图,(d)为Kaolinite的丰度图,(e)为Dumortierite的丰度图,(f)为Chalcedony的丰度图,(g)为Pyrope的丰度图,(h)为Andradite的丰度图,(i)为Montmorillonite的丰度图,(j)为Muscovite的丰度图,(k)为Alunite的丰度图.

同时为了比较端元提取效果,将RMVHU、VCA、MVES、MVSA、SISAL算法的端元提取结果与光谱库中的标准数据对比,得到各物质的光谱角距离,如表6所示.针对所有物质,各算法得到的光谱角距离的最小值已在表中加粗.

将图24与图4相比,可以发现,RMVHU算法确实能够将高光谱图像中的物质信息有效地提取出来.

从表6中可以发现,RMVHU算法在对如下6种物质:Desert Varnish、Dumortierite、Chalcedony、Pyrope、Andrad ite以及Muscovite的端元提取中,表现最为优秀.对于其他的物质,除了Sphene,另外四种物质的端元提取效果也较为优秀,RMVHU算法的端元提取结果也非常接近由其他算法所获得的最优秀的结果.最后可以发现,RMVHU算法的平均光谱角误差是所有算法中最小的,证明了本算法在实际解混应用中性能的优越性.

4 结论

本文通过引入负数惩罚正则项,提出了自适应鲁棒最小体积高光谱解混算法(RMVHU),在求解时采用循环最小化方法拆分成可以求解的凸优化子问题,并成功应用ADMM框架对子问题进行求解.此外,文章对正则项的系数λ提出了一种自适应的调节策略,加快了收敛速度,提高了算法的稳定性.

仿真实验表明相较于其他解混算法,RMVHU算法在任意端元个数下均拥有优秀且稳定的解混表现;此外,算法对数据噪声大小,数据异常点个数以及图像像素点个数亦很不敏感,尤其在端元数估计错误的情况下也有比较稳定的表现;在不同混合度的场景下RMVHU也有很优秀的表现.

对Cuprite真实数据的解混结果表明,RMVHU算法提取出的端元及估算出的丰度图与实际地物组成基本一致.

需要注意的是,RMVHU算法的计算时间相对于诸如VCA等纯像元提取算法较长,因此在下一步的研究中将尝试采用并行化技术,以提高运算的效率.

致谢

作者感谢审稿专家在论文修改过程中提出的宝贵建议和意见.作者对哈尔滨工业大学航天学院检测实验室的研究支持表示感谢.

1 Pan Zong-Xu,Yu Jing,Xiao Chuang-Bai,Sun Wei-Dong.Spectral similarity-based super resolution for hyperspectral images.Acta Automatica Sinica,2014,40(12):2797−2807(潘宗序,禹晶,肖创柏,孙卫东.基于光谱相似性的高光谱图像超分辨率算法.自动化学报,2014,40(12):2797−2807)

2 Ni Ding,Ma Hong-Bing.Spectral-spatial classi fi cation of hyperspectral images based on neighborhood collaboration.Acta Automatica Sinica,2015,41(2):273−284(倪鼎,马洪兵.基于近邻协同的高光谱图像谱-空联合分类.自动化学报,2015,41(2):273−284)

3 Du B,Zhang L P.A discriminative metric learning based anomaly detection method.IEEE Transactions on Geoscience and Remote Sensing,2014,52(11):6844−6857

4 Lin C H,Chi C Y,Wang Y H,Chan T H.A fast hyperplanebased minimum-volume enclosing simplex algorithm for blind hyperspectralunmixing.IEEE Transactions on Signal Processing,2016,64(8):1946−1961

5 Bioucas-Dias J M,Plaza A,Dobigeon N,Parente M,Du Q,Gader P,Chanussot J.Hyperspectralunmixing overview:geometrical,statistical,and sparse regression-based approaches.IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing,2012,5(2):354−379

6 Lillesand T,Kiefer R W,Chipman J.Remote Sensing and Image Interpretation(Seventh Edition).New York:John Wiley and Sons,2014.23−61

7 Yuan Y,Fu M,Lu X Q.Substance dependence constrained sparse NMF for hyperspectralunmixing.IEEE Transactions on Geoscience and Remote Sensing,2015,53(6):2975−2986

8 Wang W H,Qian Y T,Tang Y Y.Hypergraph-regularized sparse NMF for hyperspectralunmixing.IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing,2016,9(2):681−694

9 Qu Q,Nasrabadi N M,Tran T D.Subspace vertex pursuit:a fast and robust near-separable nonnegative matrix factorization method for hyperspectralunmixing.IEEE Journal of Selected Topics in Signal Processing,2015,9(6):1142−1155

10 Li J,Bioucas-Dias J M,Plaza A,Liu L.Robust collaborative nonnegative matrix factorization for hyperspectralunmixing.IEEE Transactions on Geoscience and Remote Sensing,2016,54(10):6076−6090

11 Feng R Y,Zhong Y F,Zhang L P.An improved nonlocal sparse unmixing algorithm for hyperspectral imagery.IEEE Geoscience and Remote Sensing Letters,2015,12(4):915−919

12 Iordache M D,Bioucas-Dias J M,Plaza A.Collaborative sparse regression for hyperspectralunmixing.IEEE Transactions on Geoscience and Remote Sensing,2014,52(1):341−354

13 Meyer T R,Drumetz L,Chanussot J,Bertozzi A L,Jutten C.Hyperspectralunmixing with material variability using social sparsity.In:Proceedings of 2016 IEEE International Conference on Image Processing.Phoenix,USA:IEEE,2016.2187−2191

14 Giampouras P V,Themelis K E,Rontogiannis A A,Koutroumbas K D.Simultaneously sparse and low-rank abundance matrix estimation for hyperspectral image unmixing.IEEE Transactions on Geoscience and Remote Sensing,2016,54(8):4775−4789

15 Zheng C Y,Li H,Wang Q,Chen C L P.Reweighted sparse regression for hyperspectralunmixing.IEEE Transactions on Geoscience and Remote Sensing,2016,54(1):479−488

16 Nascimento J M P,Dias J M B.Vertex component analysis:a fast algorithm to unmixhyperspectral data.IEEE Transactions on Geoscience and Remote Sensing,2005,43(4):898−910

17 Boardman J W.Automating spectral unmixing of AVIRIS data using convex geometry concepts.In:Proceedings of Summaries of the 4th Annual JPL Airborne Geoscience Workshop.Arlington,Virginia:JPL,1993.11−14

18 Winter M E.N-FINDR:an algorithm for fast autonomous spectral end-member determination in hyperspectral data.In:Proceedings of 1999 SPIE′s International Symposium on Optical Science,Engineering,and Instrumentation.Denver,USA:SPIE,1999.266−275

19 Chan T H,Chi C Y,Huang Y M,Ma W K.A convex analysis-based minimum-volume enclosing simplex algorithm for hyperspectralunmixing.IEEE Transactions on Signal Processing,2009,57(11):4418−4432

20 Li J,Agathos A,Zaharie D,Bioucas-Dias J M,Plaza A,Li X.Minimum volume simplex analysis:a fast algorithm for linear hyperspectralunmixing.IEEE Transactions on Geoscience and Remote Sensing,2015,53(9):5067−5082

21 Bioucas-Dias J M.A variable splitting augmented Lagrangian approach to linear spectral unmixing.In:Proceedings of the 1st Workshop on Hyperspectral Image and Signal Processing:Evolution in Remote Sensing.Grenoble,France:IEEE,2009.1−4

22 Boyd S,Parikh N,Chu E,Peleato B,Eckstein J.Distributed optimization and statistical learning via the alternating direction method of multipliers.Foundations and Trends in Machine Learning,2010,3(1):1−122

23 Chan T H,Ma W K,Chi C Y,Wang Y.A convex analysis framework for blind separation of non-negative sources.IEEE Transactions on Signal Processing,2008,56(10):5120−5134

24 StrangG.Linear Algebra and Its Applications.San Diego,CA:Thomson,2005.

25 Chen S S,Donoho D L,Saunders M A.Atomic decomposition by basis pursuit.SIAM Review,2001,43(1):129−159

26 Shi Z W,An Z Y,Tan X Y,Zhu Z X,Jiang Z G.Hyperspectralunmixing using non-negative matrix factorization with automatically estimating regularization parameters.In:Proceedings of the 7th International Conference on Natural Computation.Shanghai,China:IEEE,2011.1836−1840

27 Chen Bao-Lin.Optimization Theories and Algorithms(Second Edition).Beijing:Tsinghua University Press,2005.180−193(陈宝林.最优化理论与算法.第2版.北京:清华大学出版社,2005.180−193)

28 Clark R N,Swayze G A,Wise R,Livo E,Hoefen T,Kokaly R,Sutley S J.USGS digital spectral library[Online],available:http://speclab.cr.usgs.gov/spectral-lib.html,September 13,2016.

29 AVIRIS.AVIRIS data-ordering free AVIRIS standard data products[Online],available:http://aviris.jpl.nasa.gov/html/aviris.freedata.html,September 13,2016.

30 Clark R N,Swayze G A,Livo K E,Kokaly R F,Sutley S J,Dalton J B,McDougal R R,Gent C A.Imaging spectroscopy: earth and planetary remote sensing with the usgstetracorder and expert systems[Online],available:http://speclab.cr.usgs.gov/PAPERS/tetracorder,September 13,2016.

A Robust Minimum Volume Based Algorithm with Automatically Estimating Regularization Parameters for Hyperspectral Unmixing

WANG Tian-Cheng1LIU Xiang-Zhen1DONG Ze-Zheng1WANG Hai-Bo1

Hyperspectral unmixing aims at fi nding hidden endmembers and their corresponding abundances from hyperspectral images with low spatial resolution.Based on the well-known minimum volume(MV)rule in geometrical based approaches,a robust minimum volume based algorithm with automatically estimating regularization parameters for hyperspectral unmixing(RMVHU)is proposed.In this algorithm,the ANC constraint is replaced with a negative number punished regularizer which may lead to a more robust result to outliers and noise.A cyclic minimization algorithm is used to split the nonconvex RMVHU problem into convex subproblems,and ADMM is referred to sovle the large scale optimization problem with the increasing number of pixels in the image.To improve the convergence of the algorithm,a strategy to estimate the regularization parameters of the regularizer automatically is proposed.Compared with some existing geometrical based methods,experimental results show the superiority of the RMVHU algorithm on both synthetic datasets and real datasets.

Hyperspectral unmixing,alternating direction method of multipliers(ADMM),convex optimization,minimum volume,automatically estimating regularization parameters

Wang Tian-Cheng,Liu Xiang-Zhen,Dong Ze-Zheng,Wang Hai-Bo.A robust minimum volume based algorithm with automatically estimating regularization parameters for hyperspectral unmixing.Acta Automatica Sinica,2017,43(12):2141−2159

2016-09-13 录用日期2016-12-10

September 13,2016;accepted December 10,2016

本文责任编委胡清华

Recommended by Associate Editor HU Qing-Hua

1.上海卫星工程研究所上海201109

1.Shanghai Institute of Satellite Engineering,Shanghai 201109

王天成,刘相振,董泽政,王海波.一种自适应鲁棒最小体积高光谱解混算法.自动化学报,2017,43(12):2141−2159

DOI10.16383/j.aas.2017.c160653

王天成 上海卫星工程研究所助理工程师,2016年获哈尔滨工业大学工学硕士学位.主要研究方向为高光谱解混、系统建模与仿真.本文通信作者.

E-mail:wangtcsa@163.com

(WANG Tian-Cheng Assistant engineer at Shanghai Institute of Satellite Engineering.He received his master degree from Harbin Institute of Technology in 2016.His research interest covers hyperspectral unmixing,system modeling and simulating.Corresponding author of this paper.)

刘相振 上海卫星工程研究所高级工程师,2006年获上海航天技术研究院工学硕士学位.上海卫星工程研究所信息与仿真中心主任.主要研究方向为复杂系统的建模与仿真.

E-mail:sirruslxz@163.com

(LIUXiang-Zhen Senior engineer at Shanghai Institute of Satellite Engineering.He received his master degree from Shanghai Academy of Space fl ight Technology in 2006.Director of the Information and Simulation Center for Satellite Engineering,Shanghai Institute of Satellite Engineering.His research interest covers modeling and simulating of complex system.)

董泽政 上海卫星工程研究所工程师,2010年获南京航空航天大学工学硕士学位.主要研究方向为系统建模与仿真.

E-mail:dzzh520@hotmail.com

(DONG Ze-ZhengEngineer at Shanghai Institute of Satellite Engineering.He received his master degree from Nanjing University of Aeronautics and Astronautics in 2010.His research interest covers system modeling and simulating.)

王海波 上海卫星工程研究工程师,2012年获哈尔滨工业大学工学博士学位.主要研究方向为系统建模与仿真.

E-mail:13766898363@163.com

(WANG Hai-BoEngineer at Shanghai Institute of Satellite Engineering.He received his Ph.D.degree from Harbin Institute of Technology in 2012.His research interest covers system modeling and simulating.)