CN107796788A - 基于变分贝叶斯期望最大算法的传感矩阵测量方法 - Google Patents

Abstract

本发明公开了一种基于变分贝叶斯期望最大算法的传感矩阵测量方法。使用数字微透镜阵列(DMD)作光波前调制，进行标定测量y_p＝|Dx_p|，p∈{1，2，....P}；处理图像数据得到B＝|ΦA|，通过计算a_m∈A替代计算传感矩阵D；引入均场假设，执行变分贝叶斯期望最大(VBEM)算法求解一个最大后验估计问题得到a_m；组合a_m得到矩阵A，进而得到传感矩阵D。本发明使用DMD调制光波前，设备速度快、像素规模大，算法复杂度合理，能够较好得完成传感矩阵的测量。

Description

基于变分贝叶斯期望最大算法的传感矩阵测量方法

技术领域

本发明属于并行压缩感知领域里的传感矩阵测量技术，特别是一种基于变分贝叶斯期望最大算法的传感矩阵测量方法。

背景技术

压缩感知(cs)是近年来在信号采集和处理领域中一种新兴的理论方法。不同于传统采样理论的先采样再压缩，它的核心思想是直接获取压缩过的信号，即把信号压缩过程结合在成像系统中。这种方法所取得的采样数目远远小于基于香农定理的传统方法，可以显著节省信号处理传输所需要的存储空间。

Cs理论主要包括信号的稀疏表示、传感矩阵的设计和重建算法三个方面。其中，作为cs理论的核心，传感矩阵的设计得到了广泛的关注。目前，美国莱斯(rice)大学研制的单像素相机已经成为了压缩感知实现和传感矩阵构建的经典案例，随之而来的并行压缩感知技术也已经成为了研究热点。

目前，散射介质成像已经成为了并行压缩感知技术的研究前沿。当光穿过散射介质(如生物组织、ZNO薄片、白漆等)时，弹道光很快衰减并伴随着随机散射光斑模式的产生。这种光斑模式被证明可以通过光波前调制的方式人为控制。例如，可以使用SLM(空间光相位调制器)调制光波前来控制探测器端的光斑模式直至其聚焦。传感矩阵可以完全描述光从调制器件经散射介质直到探测器传输过程中的所有变化。传感矩阵的使用对于散射介质应用于并行压缩感知成像和聚焦领域具有重大的现实意义。

现有传感矩阵测量方法是使用SLM调制光波前，使用CCD测量成像系统出射光，利用四步相移技术恢复计算传感矩阵。这种方法的缺点是SLM过于昂贵且规模较小，运行速度相比于DMD(数字微透镜阵列)比较缓慢。

发明内容

本发明的目的在于提供一种基于变分贝叶斯期望最大(VBEM)算法的传感矩阵测量方法。采集和处理图像数据，在均场假设前提下把传感矩阵的计算问题转化为一个贝叶斯框架下的最大后验估计问题，进而用VBEM算法求解传感矩阵。

实现本发明目的的技术解决方案为：一种基于变分贝叶斯期望最大算法的传感矩阵测量方法，包括以下步骤：

第一步用数字微透镜阵列(DMD)分P次随机调制光波前x_p，成像并采集图像数据y_p。此过程可描述为：y_p＝|Dx_p|，p∈{1，2，....P}，D为待测传感矩阵；

第二步组合上述P个等式，得到矩阵方程Y^M×P＝|D^M×NX^N×P|。转置共轭矩阵方程可得：Y^H＝|X^HD^H|。记B＝Y^H，Φ＝X^H，A＝D^H，原矩阵方程转化为：B^P×M＝|Φ^P×NA^N×M|。由此，可通过计算矩阵A替代计算传感矩阵D。；

第三步把上步中的矩阵方程B＝|ΦA|拆分为b_m＝|Φa_m|，m∈{1，2，....M}其中矩阵A可通过分M次计算向量a_m并组合得到；

第四步为求解向量a_m，在给定的均场假设约束F下，迭代VBEM算法至Kullback-Leibler散度(KL散度)值最小，使a_m的估计值逼近实际值；

第五步重复步骤四使m历遍{1,2,…..M},得到所有a_m∈A；

第六步组合向量a_m得到矩阵A。执行转置共轭操作，传感矩阵D＝A^H。

第一步中，x_p∈R^N代表第p次经DMD调试的入射光波前，向量x_p有N个元素表示DMD有N个调制元；y_p∈R^N代表探测器第p次测得的图像数据，y_p代表强度测量即光波的相位信息在成像过程中遗失了；D∈R^M×N代表成像系统的待测传感矩阵。

第三步中，已知b_m＝|Φa_m|，通过矩阵Φ和模值向量b_m求解复向量a_m是一个经典的相位恢复问题，为求解此问题进行均场假设如下：θ的分布其中p(θ_m)＝1/2π，θ_m代表第p次测量时探测器第m个像元进行模值探测失去的相位信息。即在一次测量中，探测器上M个像素点失去的相位信息取任意值的可能性均等。

第三步中由上述均场假设，模值探测中失去的相位信息得到考虑。向量a_m的计算问题在贝叶斯框架下转化为一个最大后验估计问题：即在模值向量b_m发生的基础上，使向量a_m发生的后验概率最大，此时a_m的值即为所求值。

第四步中在上述假设和问题转化的基础上，引入约束条件F如下：

至此VBEM迭代算法可用于向量a_m的求解。每执行一次算法迭代，得到一个a_m和θ的分布q(a_m，θ_m)。以此分布作为后一次算法迭代的先验条件即后验联合分布p(a_m，θ|b_m)。当迭代算法的结果q(a_m，θ_m)与其先验p(a_m，θ|b_m)之间的KL散度值小于给定值，认为算法收敛，算法结果q(a_m，θ_m)逼近实际值。至此得到了向量a_m。

本发明与现有技术相比，其显著优点：(1)使用DMD代替SLM作为光波前调制设备，系统运行速度快，规模大，设备便宜。(2)DMD相比于SLM稳定性好，对噪声的抵抗能力强。(3)在标定测量达到一定数目时，本方法可以较好得还原出传感矩阵。

附图说明

图1是本发明的方法流程图。

图2是本发明采集图像数据所用的成像系统光路图。

图3是使用本发明求得的传感矩阵与实际传感矩阵的相关性随标定测量增加而变化的趋势。图3(a)由50次模拟仿真求均值得到；图3(b)由任1次模拟仿真得到。

具体实施方式

下面结合附图对本发明作进一步的详细描述。

本发明提出了一种基于贝叶斯框架及均场假设的传感矩阵测量方法，步骤如下：

首先，简述一下成像系统的光路设置和其具体特性。由光路图2，由激光器产生的激光经由扩束系统产生平行光打在数字微透镜阵列DMD上；DMD对入射光进行随机的二元振幅调制(0-1调制)，DMD的每一个模式状态可被认为一个光波前输入；DMD调制的波前光经透镜聚焦在散射介质上，光在介质中输运时发生随机散射，如图2右端所示光在穿过介质后形成一个输出模。此过程可以数学化得表示为y＝Dx，其中x代表经DMD调制后的光波前输入，传输矩阵D代表散射介质和透镜系统对光波的影响，y代表光波输出。

由于探测器探测到的数值是某点光强在探测器曝光时间内的积分，在此过程中，光的相位信息可以认为被丢弃了，所以成像系统进行一次图像数据采集的过程可以看作y＝|Dx|。其中x∈R^N表示DMD上有N个调制元即一个调制光波前有N个元素；y∈R^M表示探测器M个像素的强度测量值；D∈R^M×N表示传感矩阵。

当DMD处于任意给定模式x下时，探测器端有对应的特定输出模y。控制DMD，使各像元随机变化，产生一个的光波前向量x₁，探测器端对应的输出为向量y₁，且y₁＝|Dx₁|。对DMD进行P次模式控制，相应得有P个光波前向量{x₁，x₂，...x_p}和对应的模值探测向量{y₁，y₂，...y_p}，且y₁＝|Dx₁|,y₂＝|Dx₂|,…y_P＝|Dx_P|。

组合这P个等式得到矩阵方程：Y＝|DX|。其中，Y∈R^M×P＝{y₁，y₂，...y_p}为探测器端模值探测向量的组合，X∈R^N×P＝{x₁，x₂，...x_p}为已知调制光波前的组合，D∈R^M×N为待测传感矩阵。

转置共轭矩阵方程:Y^H|X^HD^H|,为便于后续计算，使B＝Y^H，Φ＝X^H，A＝D^H。原矩阵方程可化为：B^P×M＝|Φ^P×NA^N×M|。

有上述数据采集和处理过程，可以通过计算方程B＝|ΦA|中的矩阵A替代计算方程Y＝|DX|中的传感矩阵D。

其中B∈R^P×M＝{b₁，b₂，...b_M},A∈R^N×M＝{a₁，a₂，...a_M},Φ∈R^P×N。

对于矩阵方程B＝|ΦA|，可拆分为b₁＝|Φa₁|,b₂＝|Φa₂|,…b_M＝|Φa_M|。且统一表示为：b_m＝|Φa_m|其中m＝(1，2，...M)。

由此可以通过求解b_m＝|Φa_m|中的a_m，使m历遍1,2，....M得到矩阵A。

在已知Φ和模值向量b_m的情况下求解复向量a_m是一个经典的相位恢复问题，为此，引入几个概念和模型假设：

每一个绝对值元素可以表示为：

其中，θ_p∈(0，2π)表示探测器第m个像素点第p次强度测量中失去的相位信息，n_p表示第p次测量的高斯噪声且均值为0方差为在并行压缩感知系统中，传输矩阵D代表散射介质对光波造成的影响，前人对其的研究说明传输矩阵D是一个高斯独立同分布矩阵。也就是说矩阵A中的元素a_im服从高斯独立同分布，假设其方差为

为在后续模型中考虑到相位信息θ_p，进行均场假设如下：

其中p(θ_m)＝1/2π (2)

即在一次测量中，探测器上M个像素点失去的相位信息取任意值的可能性均等。

以均场假设为基础，向量a_m的计算问题可以在贝叶斯框架下表述为一个最大后验估计问题：

即在模值向量b_m发生的基础上，使向量a_m发生的后验概率最大，此时a_m的值即为所求向量值。其中

p(a_m|b_m)＝∫_θp(a_m，θ|b_m) (4)

边缘化均场假设下的隐变量θ以避免直接计算p(a_m|b_m)(若不把θ纳入模型中，该最大后验问题将难以解决)。

由上述概念和假设，给定均场假设下的约束F，可以通过最小化KL(Kullback-Leibler)散度的方式使a_m和θ的分布q(a_m，θ_m)尽可能得逼近后验联合分布p(a_m，θ|b_m)：

其中

使用VBEM(变分贝叶斯期望最大)算法解决这个逼近问题。算法流程如下：

q(a_t)＝CN(v_t，a_t) (8)

其中

其中I₁(I₀)分别代表1阶和0阶的第一类修正贝塞尔函数。.^H代表共轭转置.^*代表标量转置。

在上述VBEM算法中，(9)、(10)、(11)、(12)，(13)皆代表迭代过程中更新的参数，(7)中的q(θ_m)代表一次迭代第m个θ值得分布，(8)中的q(a_t)代表一次迭代第i个a值的分布。迭代过程取遍M个θ组成(6)中的取遍N个a组成(6)中的由此构成约束条件F，进而在约束F下计算(5)KL散度。

每运行一次迭代算法，KL散度衰减，当衰减至一个给定的数值时，认为算法收敛。至此，

p(a_m|b_m)≈∏_tq(a_mt) (14)

即a_m以b_m为条件的条件分布p(a_m|b_m)逼近了a_m的分布q(a_m)，所有元素a_mt构成a_m。

在算法收敛时，对于所有

由此得到了向量a_m。使m历遍{1，2，...M},即得到所有a_m∈A，组合a₁，a₂，…a_M得到矩阵A。

执行转置共轭操作D＝A^H，得出传感矩阵D。

图3是本方案使用matlab模拟仿真本方法得到的传感矩阵与实际传感矩阵的相关性图像。横轴代表标定测量次数P与输入调制光波前自由度N(DMD调制元个数)的比值。纵轴代表求得传感矩阵与实际传感矩阵的相关性。相关性越接近1，求得传感矩阵与实际传感矩阵越接近，本方法越成功。

由图3(a)可知，当标定测量次数较小时，本方法不能完成传感矩阵的测量；随着标定测量次数的增加，本方法所得传感矩阵与实际传感矩阵的相关性升高，最终稳定于高值，说明基本可以完成传输矩阵的测量；由图3(b)可以看出，单次测量所得结果并不准确，传感矩阵的恢复效果有起伏，说明标定质量会影响本方法恢复传感矩阵的效果。

Claims (5)

1.一种基于变分贝叶斯期望最大算法的传感矩阵测量方法，其特征在于步骤如下：

第一步，用数字微透镜阵列DMD分P次随机调制光波前x_p，成像并采集图像数据y_p，此过程描述为：y_p＝|Dx_p|，p∈{1，2，...，P}，D为待测传感矩阵；

第二步，组合上述P个等式，得到矩阵方程Y^M×P＝|D^M×NX^N×P|，转置共轭矩阵方程可得：Y^H＝|X^HD^H|；记B＝Y^H，Φ＝X^H，A＝D^H，原矩阵方程转化为：B^P×M＝|Φ^P×NA^N×M|；由此，可通过计算矩阵A替代计算传感矩阵D；

第三步，把上步中的矩阵方程B＝|ΦA|拆分为b_m＝|Φa_m|，m∈{1，2，....M}，其中矩阵A通过分M次计算向量a_m并组合得到；

第四步，在给定的均场假设约束F下，迭代VBEM算法至Kullback-Leibler散度值最小，使a_m的估计值逼近实际值，求解得到向量a_m；

第五步，重复步骤四使m历遍{1,2,.....M}，得到所有a_m∈A；

第六步，组合向量a_m得到矩阵A，执行转置共轭操作，得到传感矩阵D＝A^H。

2.根据权利要求1所述的基于变分贝叶斯期望最大算法的传感矩阵测量方法，其特征在于：第一步中，x_p∈R^N代表第p次经DMD调制的入射光波前，向量x_p有N个元素，表示DMD有N个调制元；y_p∈R^N代表探测器第p次测得的图像数据，y_p代表强度测量即光波的相位信息在成像过程中遗失了；D∈R^M×N代表成像系统的待测传感矩阵。

3.根据权利要求1所述的基于变分贝叶斯期望最大算法的传感矩阵测量方法，其特征在于：在第三步中，已知b_m＝|Φa_m|，通过矩阵Φ和模值向量b_m求解复向量a_m是一个经典的相位恢复问题，为求解此问题进行均场假设如下：θ的分布其中p(θ_m)＝1/2π，θ_m代表第p次测量时探测器第m个像元进行模值探测失去的相位信息，即在一次测量中，探测器上M个像素点失去的相位信息取任意值的可能性均等。

4.根据权利要求3所述的基于变分贝叶斯期望最大算法的传感矩阵测量方法，其特征在于：在第三步中由上述均场假设，模值探测中失去的相位信息得出，向量a_m的计算问题在贝叶斯框架下表述为一个最大后验估计问题：

即在模值向量b_m发生的基础上，使向量a_m发生的后验概率最大，此时a_m的值即为所求向量值。

5.根据权利要求1-4所述的基于变分贝叶斯期望最大算法的传感矩阵测量方法，其特征在于：第四步中在上述假设和问题转化的基础上，引入约束条件F如下：至此VBEM迭代算法用于向量a_m的求解；每执行一次算法迭代，得到一个a_m和θ的分布q(a_m，θ_m)，以此分布作为后一次算法迭代的先验条件即后验联合分布p(a_m，θ|b_m)；当迭代算法的结果q(a_m，θ_m)与其先验p(a_m，θ|b_m)之间的KL散度值小于给定值，则认为算法收敛，算法结果q(a_m，θ_m)逼近实际值，至此得到了向量a_m。