光电工程  2018, Vol. 45 Issue (4): 170617      DOI: 10.12086/oee.2018.170617     
大气湍流波前压缩感知测量重建研究
李灿1 , 蔡冬梅2 , 贾鹏2 , 刘建霞1 , 李娟娟2     
1. 太原理工大学信息工程学院,山西 太原 030024;
2. 太原理工大学物理与光电工程学院,山西 太原 030024
摘要:压缩感知技术用于大气湍流波前斜率测量能在很大程度上提高波前信号的测量速度,同时降低波前测量系统的硬件压力。与现有波前斜率测量方法不同,压缩感知波前测量方法增加了从波前斜率的稀疏测量值到波前斜率信号的重建过程,因此将压缩感知技术用于波前测量,需要快速、高精度的波前斜率重建算法。Smoothed L0 Norm (SL0)算法是一种近似L0范数估计的优化迭代重建算法,与其它算法相比,不需要事先知道信号的稀疏度,计算量低且估计精度高。本文以SL0算法为基础,对波前斜率信号分区域测量,再结合并行运算,通过理论分析和仿真实验实现了一种能够快速、高精度重建信号的分区域并行算法—Block-Smoothed L0 Norm (B-SL0)。实验结果表明,B-SL0在计算时间和精度都明显优于现有的其它重建算法,对压缩感知技术用于大气湍流波前测量的可行性进行了初步探索。
关键词压缩感知    波前斜率    SL0算法    分区域    并行运算    
Research on reconstruction of atmospheric turbulence wavefront compressed sensing measurement
Li Can1, Cai Dongmei2, Jia Peng2, Liu Jianxia1, Li Juanjuan2     
1. Institute of Information Engineering, Taiyuan University of Technology, Taiyuan, Shanxi 030024, China;
2. Institute of Physics and Optoelectronics, Taiyuan University of Technology, Taiyuan, Shanxi 030024, China
Abstract: Compressed sensing technology for atmospheric turbulence wavefront slope measurement can greatly improve the wavefront signal measurement speed, while reducing the pressure of wavefront measurement system hardware. Different from the existing wavefront slope measurement method, the compressed sensing wavefront measurement increase a process which from sparse measurement of wavefront slope value to the reconstruction of the wavefront slope signal. Therefore, a fast and accurate wavefront slope reconstruction algorithm is needed if the compressed sensing technology is used for wavefront measurement. Smoothed L0 Norm (SL0) algorithm is an optimized iterative reconstruction algorithm with approximate L0 norm estimation, and compared with other algorithms, it is not necessary to know the sparsity of the signal in advance, and the calculation is low and the estimation accuracy is high. Based on the SL0 algorithm, this paper implements a subregion parallel algorithm-Block-Smoothed L0 Norm (B-SL0) which can quickly and accurately reconstruct the signal by measuring the wavefront slope signal in subarea and parallel operations through theoretical analysis and experiments. The experimental results show that B-SL0 is significantly better than other existing reconstruction algorithms in the calculation time and accuracy, and explore the feasibility of compressed sensing technology for measurement of atmospheric turbulence wavefront preliminarily.
Keywords: compressed sensing    wavefront slope    SL0 algorithm    subregion    parallel operation    

1 引言

为了克服大气湍流干扰,大型地基天文望远镜在天文成像观测中都陆续配备了自适应光学系统[1-4]。波前探测器能够实时测量大气湍流波前动态误差,是自适应光学系统的关键器件之一。Shack-Hartmann传感器(Shack-Hartmann wavefront sensor, SHWFS)[5]是一种常用的波前测量器,它通过规则排列的微透镜阵列实时测量大气湍流相位ϕ(x, y)的斜率ϕx, ϕy来实现大气湍流波前相位的测量。随着天文光学望远镜的口径的增大,成像分辨率的提高,对大气湍流波前测量系统提出了更高的要求。Rostami等人[6-7]提出将压缩感知(compressed sensing, CS)技术[8-12]应用在大气湍流波前测量上,即波前斜率压缩感知(derivative compressed sensing, DCS)技术。DCS通过对波前斜率的稀疏测量,再通过非线性重建算法精确重建波前斜率分布,最后恢复出大气湍流波前相位。压缩感知技术的应用,突破了Nyquist采样定律,可以大大减少波前斜率的测量数目,提高测量速度,减轻波前测量系统的硬件压力和数据传输与存储的压力。与现有波前斜率测量方法不同,大气湍流波前斜率的压缩感知测量,增加了一个从波前斜率的稀疏测量值到波前斜率的重建过程,这会导致波前数据处理的时间加长。由于大气湍流波前测量属于实时测量,所以需要在保证波前斜率信号重建精度的同时,加快斜率的重建速度。

参考文献[6-7]对压缩感知技术应用于大气湍流波前测量的可行性进行了研究,但并未涉及波前斜率稀疏信号的重建。波前斜率能否快速、高精度重建对压缩感知技术在波前测量上的应用至关重要。本文针对上述问题,对应用在大气湍流波前斜率重建上的压缩感知重建算法进行研究。随着压缩感知技术的发展和应用,目前有多种重建算法,如匹配追踪(matching pursuit,MP)算法[13]、正交匹配追踪(orthogonal matching pursuit,OMP)算法[14]、基追踪(basic pursuit,BP)算法[15]、梯度投影法[16]、Smoothed L0 Norm (SL0)算法[17-18]等。其中Smoothed L0 Norm (SL0)算法是一种近似L0范数估计的优化迭代重建算法,与其它类型算法相比,它不需要提前知道信号的稀疏度,计算量低且估计精度高。本文在DCS[6]基础上,对波前斜率的稀疏测量信号进行分区域处理,结合SL0算法与并行运算,实现一种分区域并行波前斜率重建算法(B-SL0)。通过仿真实验对B-SL0算法进行了详细的理论推导和性能分析,实验结果表明,B-SL0算法不仅大大缩短了重建运行时间,还提高了相位ϕ(x, y)的重建精度。

2 波前斜率压缩感知测量原理

定义ϕ(x, y)为大气湍流波前相位分布,其在xy方向的斜率信号分别为ϕx, ϕy,将ϕx, ϕy在变换域即稀疏字典ΨRM×N下进行稀疏表示为

$ {\phi _x} = \mathit{\boldsymbol{ \boldsymbol{\varPsi} }}{\mathit{\boldsymbol{\theta }}_x},\;\;{\phi _y} = \mathit{\boldsymbol{ \boldsymbol{\varPsi} }}{\mathit{\boldsymbol{\theta }}_y}\;\;, $ (1)

式中:θxRM×NθyRM×N为稀疏矩阵,即矩阵的非零元素数目满足Kx << N2Ky << N2。设计一个大小为M×N且满足RIP准则[19-20]的观测矩阵Ф(其中M << N)对波前斜率进行稀疏测量,得到斜率测量值fxfy,如式(2)所示:

$ {f_x} = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\phi _x}, \;{f_y} = \mathit{\boldsymbol{ \boldsymbol{\varPhi} }}{\phi _y}\;\;\;, $ (2)

式中:fxRM×NfyRM×N

结合式(1)与式(2)得到式(3):

$ {f_x} = \mathit{\boldsymbol{ \boldsymbol{\varPhi} \boldsymbol{\varPsi} }}{\mathit{\boldsymbol{\theta }}_x} = \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{\theta }}_x}, \;{f_y} = \mathit{\boldsymbol{ \boldsymbol{\varPhi} \boldsymbol{\varPsi} }}{\mathit{\boldsymbol{\theta }}_y} = \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{\theta }}_y}\;\;, $ (3)

式中A=ФΨRM×N是感知矩阵。

波前斜率的压缩感知测量过程可由图 1表示。

图 1 波前斜率压缩感知测量过程 Fig. 1 The measurement process of wavefront slope compressed sensing

图 1所示,SHWFS是通过一个N×N规则排布的微透镜阵列测量波前斜率值ϕx, ϕy,而压缩感知测量是先通过观测矩阵Φ得到稀疏随机测量值fxfy,再根据fxfy利用重建算法重建波前斜率,最后恢复出大气湍流波前相位分布ϕ(x, y)。

压缩感知波前测量方法可以有效地降低波前测量值数目,提高测量速度,同时降低探测器的硬件以及数据的传输与存储要求。然而与SHWFS的测量过程相比,增加了从稀疏测量值fxfyϕx, ϕy的重建过程。根据大气湍流波前测量实时、高精度的应用需求,从fxfyϕx, ϕy的重建过程要求计算速度快,精度高,这样压缩感知技术在大气湍流波前测量的应用才有实际的工程意义。

3 B-SL0算法原理及实现

下面以波前斜率ϕx的重建过程为例对重建算法进行讨论。对波前斜率ϕx的稀疏测量值fx重建求解的核心问题可以表示为

$ \mathop {{\rm{min}}}\limits_{{\mathit{\boldsymbol{\theta }}_x}} {\left\| {{\mathit{\boldsymbol{\theta }}_x}} \right\|_0}\;\;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;{f_x} = \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{\theta }}_x}\;\;。$ (4)

对于式(4),定义函数l(θx)为

$ l({\mathit{\boldsymbol{\theta }}_x}) = \left\{ {\begin{array}{*{20}{c}} 1&{{\mathit{\boldsymbol{\theta }}_x} = 1}\\ 0&{{\mathit{\boldsymbol{\theta }}_x} = 0} \end{array}} \right.\;\;。$ (5)

因此,L0范数可以等价为式(6):

$ {\left\| {{\mathit{\boldsymbol{\theta }}_x}} \right\|_0} = L({\mathit{\boldsymbol{\theta }}_x}) = \sum\limits_{s = 1}^N {l({\mathit{\boldsymbol{\theta }}_{{x_s}}})}\;\;。$ (6)

因为L0范数求解是一个NP-hard问题,无法直接求解,所以出现了很多重建算法。根据大气湍流波前测量实时性的要求,通过对多个重建算法的分析,SL0算法无需提前知道信号的稀疏度,计算量小、运算速度快,适合用于波前斜率的重建。

3.1 波前斜率信号分区域并行处理

大气湍流波前斜率信号属于二维分布。对于二维信号,现有的重建算法大部分都是将其转化为一维信号或按列依次重建,通过串行计算进行重建,这样不但使重建算法计算时间大大增加,还会破坏波前斜率信号列与列之间的联系,降低波前斜率重建精度,影响波前相位的恢复。针对这些问题,本文考虑分区域重建波前斜率信号,以x方向上的斜率ϕx为例,由于SHWFS是通过一个N×N规则排布的微透镜阵列测量波前斜率值ϕx, ϕy,我们将N×N规则排布的微透镜阵列分割成nB×B大小的区域,其中n=N2/B2,然后得到nB×B大小的ϕx的子集ϕxi,由于SL0算法针对一维信号重建,故将ϕxi每列首尾相连,得到nB2×1大小的列矩阵yfi,将nyfi向量合并在一起得到处理后的波前斜率信号yf,最后通过随机采样Фf得到yf的稀疏随机测量值Yf。虽然本文也采用将二维信号转化为一维信号处理的方法,但是并没有完全破坏信号列与列之间的联系,在一定程度上较SL0算法减轻了对信号内部关系的破坏,提高了波前相位重建精度。分区域处理后的波前斜率信号yf图 2所示。

图 2 分区域处理后的波前斜率信号yf Fig. 2 Wavefront slope signal yf after subarea processing

所以式(4)转化为式(7):

$ \mathop {{\rm{min}}}\limits_{{\mathit{\boldsymbol{\theta }}_f}} {\left\| {{\mathit{\boldsymbol{\theta }}_f}} \right\|_0}\;\;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;{\mathit{\boldsymbol{Y}}_f} = {\mathit{\boldsymbol{A}}_f}{\mathit{\boldsymbol{\theta }}_f}\;\;, $ (7)

式中Af=ФfΨ

对于每一块yfi而言,式(4)转化为式(8):

$ \mathop {{\rm{min}}}\limits_{{\mathit{\boldsymbol{\theta }}_{fi}}} {\left\| {{\mathit{\boldsymbol{\theta }}_{fi}}} \right\|_0}\;\;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;{\mathit{\boldsymbol{Y}}_{fi}} = {\mathit{\boldsymbol{A}}_f}{\mathit{\boldsymbol{\theta }}_{fi}}\;\;, $ (8)

式中i=1, 2, …, n

对于每个区域而言,彼此之间是互不相关的,所以得到的波前斜率ϕxi也是互不依赖的,即处理后的斜率信号yf的列与列之间互不依赖,能同时对yf每列利用SL0算法并行重建再合并为处理后的斜率重建值${\mathit{\boldsymbol{\hat y}}_f}$,再通过逆变换得到最终的波前斜率重建信号${\hat \phi _x}$,这样不但提高信号的重建速度,还在一定程度上保护了斜率信号列与列之间的关联,提高了重建精度。

3.2 B-SL0算法原理

B-SL0算法的基础是SL0算法,SL0算法属于一种近似L0范数估计方法,基本思想是利用一个连续代理函数来逼近L0范数。定义一个标准高斯函数${f_\sigma }({\mathit{\boldsymbol{\theta }}_{{f_i}}})$作为处理后的波前斜率yf的某一列yfi的L0范数的近似估计函数,如式(9)所示:

$ {f_\sigma }({\mathit{\boldsymbol{\theta }}_{{f_i}}}) = \exp (-\frac{{\mathit{\boldsymbol{\theta }}_{{f_i}}^2}}{{2{\sigma ^2}}})\;\;, $ (9)

式中σ是一个常数参变量。

从式(9)可以得出:

$ \mathop {\lim }\limits_{\sigma \to 0} {f_\sigma }({\mathit{\boldsymbol{\theta }}_{{f_i}}}) = \left\{ {\begin{array}{*{20}{c}} 0&{{\mathit{\boldsymbol{\theta }}_{{f_i}}} \ne 0}\\ 1&{{\mathit{\boldsymbol{\theta }}_{{f_i}}} = 0} \end{array}} \right.\;\;, $ (10)

${F_\sigma }({\mathit{\boldsymbol{\theta }}_{{f_i}}}) = \sum\limits_{s = 1}^{{B^2}} {{f_\sigma }({\mathit{\boldsymbol{\theta }}_{{f_i}_s}})} $

进而得到:

$ {\left\| {{\mathit{\boldsymbol{\theta }}_{{f_i}}}} \right\|_0} = L({\mathit{\boldsymbol{\theta }}_{{f_i}}}) = N-\mathop {\lim }\limits_{\sigma \to 0} {F_\sigma }({\mathit{\boldsymbol{\theta }}_{{f_i}}})。$ (11)

由此可以看出,式(4)的求解就是求${F_\sigma }({\mathit{\boldsymbol{\theta }}_{{f_i}}})$函数的最大值。显然σ值越小,式(10)越逼近L0范数,但同时会造成函数的不光滑,从而可能陷入局部最大值;反之,σ值越大,式(9)越光滑,包含的局部最大值也就越少,式(10)近似L0范数的程度也就越差。SL0算法采用逐步降低σ值的方法来解决上述问题,即首先设置σ为一个足够大的值σ0来确保此时的Fσ平滑且包含很少的局部最大值,然后赋予σ值一个逐渐降低的序列,多次循环直到达到设置的截止条件为止。SL0算法重建处理后的波前斜率信号yf的某一列yfi的主要步骤为

输入:观测向量Yfi(i=1, 2, …, n),Af

运算:

1) 设置θfi的初始值${{\mathit{\boldsymbol{\hat \theta }}}_{{f_{{i_0}}}}}$${\mathit{\boldsymbol{Y}}_{{f_i}}} = {\mathit{\boldsymbol{A}}_f}{\mathit{\boldsymbol{\theta }}_{{f_i}}}$的最小二乘解,即${{\mathit{\boldsymbol{\hat \theta }}}_{{f_{{i_0}}}}} = {({\mathit{\boldsymbol{A}}_f}^{\rm{T}}{\mathit{\boldsymbol{A}}_f})^{-1}}{\mathit{\boldsymbol{A}}_f}^{\rm{T}}{\mathit{\boldsymbol{Y}}_{{f_i}}}$

2) 设置逐渐减小的参数σ的序列{σ1, σ2σJ};

对于j=1, 2, …, J

① 令σ=σj

② 使用最速下降法在可行域

$ \mathit{\Gamma } = \left\{ {{\mathit{\boldsymbol{\theta }}_{{f_i}}}\left| {{\mathit{\boldsymbol{Y}}_{{f_i}}} = {\mathit{\boldsymbol{A}}_f}{\mathit{\boldsymbol{\theta }}_{{f_i}}}} \right.} \right\}; $

来最大化Fσ,并迭代运行K次;

初始化:令${\mathit{\boldsymbol{\theta }}_{{f_i}}} = {{\mathit{\boldsymbol{\hat \theta }}}_{{f_{{i_{j-1}}}}}}$

对于k=1, 2, …, K

a) 令

$ \delta = {[{\mathit{\boldsymbol{\theta }}_f}_{{i_1}}{f_\sigma }({\mathit{\boldsymbol{\theta }}_f}_{{i_1}}), {\mathit{\boldsymbol{\theta }}_f}_{{i_2}}{f_\sigma }({\mathit{\boldsymbol{\theta }}_f}_{{i_2}}), \cdots, {\mathit{\boldsymbol{\theta }}_f}_{{i_{{B^2}}}}{f_\sigma }({\mathit{\boldsymbol{\theta }}_f}_{{i_{{B^2}}}})]^{\rm{T}}}; $

b) 设置步长μ(为一正常数),更新:

$ {\mathit{\boldsymbol{\theta }}_{{f_i}}} = {\mathit{\boldsymbol{\theta }}_{{f_i}}}-\mu \delta ; $

c) 将θfi投影到可行域Γ上,即:

$ {\mathit{\boldsymbol{\theta }}_{{f_i}}} = {\mathit{\boldsymbol{\theta }}_{{f_i}}}- {\mathit{\boldsymbol{A}}_f}^{\rm{T}}{(\mathit{\boldsymbol{A}}{}_f{\mathit{\boldsymbol{A}}_f}{\rm{T}})^{- 1}}({\mathit{\boldsymbol{A}}_f}{\mathit{\boldsymbol{\theta }}_{{f_i}}}- {\mathit{\boldsymbol{Y}}_{{f_i}}})。$

③ 设置θfij=θfi

输出:重构信号${\mathit{\boldsymbol{\hat \theta }}_{{f_{{i_j}}}}} = {\mathit{\boldsymbol{\hat \theta }}_{{f_i}}}$

从上述可知SL0算法存在两个循环,其中内循环采用最速下降法求解目标函数Fσ的极大值,并且采用固定的步长因子μ。为了提高速度,并不需要等到最速下降法到达收敛,即对每一个σ值,并不用得到Fσ的精确最大值,只需要进入Fσ最大化附近的区域,避免进入局部最大值即可,从而加快算法整体的运行速度。

B-SL0算法流程图如图 3所示。首先通过观测矩阵Фf处理后的波前斜率yf稀疏测量,得到观测数据Yf,后利用并行运算结合SL0算法重构稀疏矩阵θf估计值${{\mathit{\boldsymbol{\hat \theta }}}_f}$,再通过逆矩阵变换得到分区域处理信号yf的估计值${{\mathit{\boldsymbol{\hat y}}}_f}$。由于得到的是各分区域处理后的波前斜率的估计值,需要将${{\mathit{\boldsymbol{\hat y}}}_f}$进行逆变换,得到波前斜率信号ϕx的估计值${{\hat \phi }_x}$,同理可得到ϕy的估计值${{\hat \phi }_y}$,最后利用${{\hat \phi }_x}$${{\hat \phi }_y}$恢复出波前相位分布ϕ(x,y)的估计值$\hat \phi (\mathit{x}{\rm{, }}\mathit{y})$

图 3 B-SL0算法流程图 Fig. 3 The flow chart of B-SL0 algorithm
4 实验结果与分析

为了说明B-SL0算法对大气湍流波前斜率的重建性能,本文基于MATLAB R2016a处理平台进行了仿真实验。将细分离散余弦基(discrete cosine transform, DCT)作为波前斜率的稀疏字典Ψ,高斯随机矩阵作为测量矩阵Ф,设每个区域的大小B=16。利用功率谱反演[21-22]得到大气湍流随机相位分布ϕx, ϕy作为测量目标,通过对斜率信号的稀疏测量,利用B-SL0算法重建波前斜率值,最后重构出波前相位$\hat \phi (\mathit{x}{\rm{, }}\mathit{y})$,并与ϕ(x,y)比较,进而评价B-SL0算法的性能。

图 4为相干长度r0=0.05,相位屏尺寸D=2,分辨率为256×256的大气湍流相位屏的仿真结果,图 4(a)为大气湍流相位分布ϕ(x, y)。在压缩比r=0.3(压缩比r=M/N)下,对波前斜率稀疏测量,分别利用B-SL0和SL0重建波前斜率,进而重构波前相位$\hat \phi (\mathit{x}{\rm{, }}\mathit{y})$,见图 4(b)4(d)。比较$\hat \phi (\mathit{x}{\rm{, }}\mathit{y})$ϕ(x, y),得到重建相位分布的误差,见图 4(c)4(e)。可以明显的看出B-SL0算法重建波前斜率恢复出的相位屏和原图更为相近,误差小;SL0算法重建结果恢复出的相位屏与原图有明显的差别,误差大,所以B-SL0算法具有更好的重建性能。

图 4 r=0.3时,B-SL0和SL0重建和重建误差对比图。(a)原图;(b) B-SL0重建图;(c) B-SL0重建误差图;(d) SL0重建图;(e) SL0重建误差图 Fig. 4 The reconstruction comparison of B-SL0 and SL0 when r=0.3. (a) Original image; (b) Reconstruction of B-SL0; (c) Reconstruction error of B-SL0; (d) Reconstruction of SL0; (e) Reconstruction error of SL0

为了定量分析SL0与B-SL0算法的重建性能,实验对大气湍流波前相位ϕ(x, y)进行了50次循环重建,计算其均方误差MSE与峰值信噪比PSNR的平均值。图 5为不同压缩比下两种算法的比较结果。当压缩比r在0.2~0.5之间时,B-SL0算法的重建相位屏的PSNR高出SL0算法的1倍左右,在r=0.3时,PSNR高出了接近20 dB,MSE则低于SL0算法的22倍左右。图 6(a)6(b)为波前斜率ϕx, ϕy的重建时间比较,B-SL0算法的运行时间远低于SL0算法,且随着压缩比的上升,B-SL0算法的运行时间的变化相对比较平稳,基本都在0.1 s~0.9 s之间。B-SL0算法无论是在重建精度还是运行时间上都优于SL0算法。

图 5 合成相位ϕ峰值信噪比与均方误差分析曲线 Fig. 5 The analysis curve about curve peak signal to noise ratio and mean square error of synthetic phase ϕ

图 6 波前斜率的重建时间 Fig. 6 Reconstruction time of wavefront slope

为了证明B-SL0算法的一般性,模拟生成50幅大气湍流随机相位屏,在同一压缩比r下,分别采用B-SL0算法与SL0算法重建,并把50幅图的均方误差MSE与峰值信噪比PSNR的平均值和重建时间进行对比,见表 1

表 1 不同相位屏的SL0与B-SL0重建质量对比 Table 1 The reconstruction quality comparison of SL0 and B-SL0 algorithm under different phase screens
Algorithms PSNR/dB MSE Running time/s
ϕx ϕy
B-SL0 48.48 0.43 0.27 0.24
SL0 35.25 9.70 1.55 1.52

表 2中可以看出,B-SL0算法重建结果的PSNR比SL0算法高出13 dB左右,MSE低了接近23倍,运行时间也比SL0算法降低了5倍左右。为了更加直观体现B-SL0算法的性能,我们从这50幅不同的相位屏中任取一幅,将其原图与B-SL0算法、SL0算法分别重建的相位屏图做对比,如图 7所示,其中图 7(a)为相位屏原图,7(b)为B-SL0算法的相位屏重建图,7(c)为B-SL0算法的相位屏重建误差图 7,(d)为SL0算法的相位屏重建图,7(e)为SL0算法的相位屏重建误差图。从图中可以明显看出B-SL0比SL0算法取得了更好的重建效果。

表 2 各算法重建质量对比 Table 2 The reconstruction quality comparison of different algorithms
Algorithms PSNR/dB MSE Running time/s
ϕx ϕy
B-SL0 56.67 0.40 0.26 0.24
SL0 45.09 5.80 1.56 1.54
OMP 44.09 7.30 4.24 4.81
SP 32.73 99.88 1.59 1.59
BP 34.13 72.27 24.26 30.96

图 7 任意一幅相位屏重建对比图。(a)原图;(b) B-SL0重建图;(c) B-SL0重建误差图;(d) SL0重建图;(e) SL0重建误差图 Fig. 7 The reconstruction comparison of any one phase screen. (a) Original image; (b) Reconstruction of B-SL0; (c) Reconstruction error of B-SL0; (d) Reconstruction of SL0; (e) Reconstruction error of SL0

为了进一步分析B-SL0算法在重建精度和运行时间上的性能,我们将B-SL0与一些经典重建算法进行对比。表 2是在压缩比r=0.3的情况下,各种算法重建性能比较。表 2表明,无论是峰值信噪比PSNR,还是均方误差MSE,亦或是运行时间,B-SL0的重建相位屏的效果都要优于其他几种算法。

最后研究了B-SL0算法的鲁棒性,在压缩比r=0.3的条件下,比较了各算法在不同的信噪比SNR下的重建性能。图 8是不同算法在不同信噪比下的峰值信噪比PSNR与均方误差MSE的比较曲线。

图 8 不同信噪比下的PSNR与MSE对比图 Fig. 8 Comparison of PSNR and MSE under different signal to noise ratio

从图中可以明显看出相比较其它几种重建算法,B-SL0算法在不同的噪声条件下,重建效果依然比较平稳,体现了B-SL0算法的鲁棒性。

5 总结

将压缩感知技术应用在大气湍流波前斜率重建上,虽然能够很大程度上降低对波前测量硬件系统的压力,减少数据传输与存储的压力,但是也会增加波前数据处理时间,这对压缩感知重建算法提出了更高的要求。本文在DCS的基础上,针对其重构精度和运行速度两方面的不足,将SL0算法与分区域并行运算相结合,提出了分区域并行算法—B-SL0算法。由于SL0算法是针对一维信号重建,对于波前斜率这类二维信号采用逐列串行重建方法,一方面这属于串行运算,增加了算法重建时间;另一方面会破坏波前斜率信号列与列之间的关联,降低了波前斜率重建精度,B-SL0算法则采用分区域并行运算,不但减少了重建算法的运行时间,还降低了对斜率信号列与列之间的联系的破坏,提高了波前相位的重建精度。仿真实验结果表明,B-SL0算法不但在波前斜率的运行时间上优于SL0算法,而且其重建的波前斜率恢复出的波前相位精度也要好于SL0算法,此外,与一些经典算法相比(如OMP、SP和BP等),在相同的条件下该算法不但在运行时间上有了很大的改进,而且其重建出的斜率还能够更好地恢复大气湍流波前相位ϕ(x, y),体现了该算法在重建相位屏ϕ(x, y)上的良好的性能。从算法本身来讲,本文的并行运算是利用MATLAB软件工具箱自带的并行运算函数,未对其并行运算部分采用GPU加速进行优化。本文的工作主要是对将压缩感知用于大气湍流波前测量的可行性进行了初步探索,对于实际的工程应用还有相当的距离,需要进一步改进算法本身的性能,提高计算速度和重建精度。

参考文献
[1]
Lin X D, Xue C, Liu X Y, et al. Current status and research development of wavefront correctors for adaptive optics[J]. Chinese Optics, 2012, 5(4): 337-351.
林旭东, 薛陈, 刘欣悦, 等. 自适应光学波前校正器技术发展现状[J]. 中国光学, 2012, 5(4): 337-351.
[2]
Xian H. Design and optimization of wavefront sensor for adaptive optics system[D]. Chengdu: University of Electronic Science and Technology of China, 2008.
鲜浩. 自适应光学系统波前传感器设计与优化[D]. 成都: 电子科技大学, 2008.
[3]
Niu C J, Yu S J, Han X E. Analysis about effect of wavefront sensorless adaptive optics on optical communication[J]. Laser & Optoelectronics Progress, 2015, 52(8): 080102.
牛超君, 于诗杰, 韩香娥. 无波前探测自适应光学对光通信性能影响分析[J]. 激光与光电子学进展, 2015, 52(8): 080102.
[4]
Zhang Q, Jiang W H, Xu B. Study of zonal wavefront reconstruction adapting for Hartmann-Shack wavefront sensor[J]. High Power Laser and Particle Beams, 1998, 10(2): 229-233.
张强, 姜文汉, 许冰. 用于Hartmann-Shack波前探测器的区域法算法研究[J]. 强激光与粒子束, 1998, 10(2): 229-233.
[5]
Yazdani R, Fallah H. Wavefront sensing for a Shack-Hartmann sensor using phase retrieval based on a sequence of intensity patterns[J]. Applied Optics, 2017, 56(5): 1358-1364. DOI:10.1364/AO.56.001358
[6]
Rostami M, Michailovich O, Wang Z. Image deblurring using derivative compressed sensing for optical imaging application[J]. IEEE Transactions on Image Processing, 2012, 21(7): 3139-3149. DOI:10.1109/TIP.2012.2190610
[7]
Polans J, Mcnabb R P, Izatt J A, et al. Compressed wavefront sensing[J]. Optics Letters, 2014, 39(5): 1189-1192. DOI:10.1364/OL.39.001189
[8]
Donoho D L. Compressed sensing[J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289-1306. DOI:10.1109/TIT.2006.871582
[9]
Ren Y M, Zhang Y N, Li Y. Advances and perspective on compressed sensing and application on image processing[J]. Acta Automatica Sinia, 2014, 40(8): 1563-1575.
任越美, 张艳宁, 李映. 压缩感知及其图像处理应用研究进展与展望[J]. 自动化学报, 2014, 40(8): 1563-1575.
[10]
Tsaig Y, Donoho D L. Extensions of compressed sensing[J]. Signal Processing, 2006, 86(3): 549-571. DOI:10.1016/j.sigpro.2005.05.029
[11]
Candès E J, Romberg J, Tao T. Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information[J]. IEEE Transactions on Information Theory, 2006, 52(2): 489-509. DOI:10.1109/TIT.2005.862083
[12]
Candes E J, Tao T. Near-optimal signal recovery from random projections: universal encoding strategies[J]. IEEE Transactions on Information Theory, 2006, 52(12): 5406-5425. DOI:10.1109/TIT.2006.885507
[13]
Mallat S G, Zhang Z F. Matching pursuits with time-frequency dictionaries[J]. IEEE Transactions on Signal Processing, 1993, 41(12): 3397-3415. DOI:10.1109/78.258082
[14]
Tropp J A, Gilbert A C. Signal recovery from random measurements via orthogonal matching pursuit[J]. IEEE Transactions on Information Theory, 2007, 53(12): 4655-4666. DOI:10.1109/TIT.2007.909108
[15]
Yin W, Morgan S, Yang J F, et al. Practical compressive sensing with Toeplitz and circulant matrices[J]. Proceedings of SPIE, 2010, 7744: 77440K.
[16]
Applebaum L, Howard S D, Searle S, et al. Chirp sensing codes: Deterministic compressed sensing measurements for fast recovery[J]. Applied & Computational Harmonic Analysis, 2009, 26(2): 283-290.
[17]
Mohimani H, Babaie-zadeh M, Jutten C. A fast approach for overcomplete sparse decomposition based on smoothed 10 norm[J]. IEEE Transactions on Signal Processing, 2009, 57(1): 289-301. DOI:10.1109/TSP.2008.2007606
[18]
Mohimani G H, Babaie-zadeh M, Jutten C. Fast sparse representation based on smoothed l0 norm[C]//Proceedings of the 7th International Conference on Independent Component Analysis and Signal Separation. Springer-Verlag, 2007: 389–396.
[19]
Candes, Emmanuel J. The restricted isometry property and its implications for compressed sensing[J]. Comptes rendus-Mathematique, 2008, 346(9): 589-592.
[20]
Cai T T, Wang L, Xu G W. New bounds for restricted isometry constants[J]. IEEE Transactions on Information Theory, 2010, 56(9): 4388-4394. DOI:10.1109/TIT.2010.2054730
[21]
Cai D M, Wang K, Jia P, et al. Sampling methods of power spectral density method simulating atmospheric turbulence phase screen[J]. Journal of physics, 2014, 63(10): 104217.
蔡冬梅, 王昆, 贾鹏, 等. 功率谱反演大气湍流随机相位屏采样方法的研究[J]. 物理学报, 2014, 63(10): 104217. DOI:10.7498/aps.63.104217
[22]
Zhang Z L. Research on the simulation system of indoor atmospheric turbulence[D]. Taiyuan: Taiyuan University of Technology, 2017.
张智露. 室内大气湍流模拟系统的研究[D]. 太原: 太原理工大学, 2017.
[23]
Li Y J, Zhu W Y, Rao R Z. Simulation of random phase screen of non-Kolmogorov atmospheric turbulence[J]. Infrared and Laser Engineering, 2016, 45(12): 1211001.
李玉杰, 朱文越, 饶瑞中. 非Kolmogorov大气湍流随机相位屏模拟[J]. 红外与激光工程, 2016, 45(12): 1211001.