CN111352126A - 一种基于大气散射介质调制的单像素成像方法 - Google Patents

Abstract

本发明公开了一种基于大气散射介质调制的单像素成像方法，包括：激光束进入大气散射介质，以大气散射介质作为空间调制器，计算当前条件下的大气散射介质的单次调制样式；在该调制样式下，大气反射方向的光经过透镜到达目标，桶探测器收集所有来自目标的光信号，获得探测值；随机改变激光束进入大气的传输方向，获得大气对光场调制的一系列不同的调制样式；以单次调制样式为行向量，将一系列不同的调制样式对应的行向量组成大气散射介质测量矩阵；并获得一系列不同的调制样式下对应的一系列探测值组成的测量向量；利用图像重构算法恢复目标信号。本发明的方法利用大气散射介质作为单像素成像的空间调制器，实现被散射介质包围的目标的清晰成像。

Description

一种基于大气散射介质调制的单像素成像方法

技术领域

本发明涉及光学成像领域，具体涉及一种基于大气散射介质调制的单像素成像方法。

背景技术

光学成像在对地观测、自动驾驶等领域发挥了重要作用，散射介质给光学成像带来了挑战。散射介质改变了光波的传输方向，造成强度衰减的同时，改变了光场的空间分布，导致无法通过散射介质对目标清晰成像。近年来光学相机的技术指标不断提高，具有更高的空间分辨率、更优的调制传递函数、更小的重量/尺寸，然而散射介质的存在，严重降低了相机获取的图像质量，限制了光学相机的探测能力。

近年来，针对散射条件下光学成像的新理论、新技术和新方法不断涌现。当前的研究方法大致可以分为两大类：抑制散射介质的影响和利用散射介质的影响。抑制散射介质影响的核心思想是找到并分离出未受散射介质干扰的光子，即所谓“弹道光子”，采用的方法包括时间分辨、关联成像、压缩感知等。由于弹道光子的数量随传播距离呈指数下降，因此这些方法均不适用于长距离条件下的成像。利用散射介质的影响以“记忆效应”和“透过矩阵”为代表。记忆效应是指光波从散射介质的不同方向穿过，形成的散斑形状不变，仅位置发生偏移。基于记忆效应，出现了波前整形、散斑相关和频谱测量等方法。记忆效应与散射介质的厚度成反比，因此这些方法只适用于薄的散射介质。基于透过矩阵的方法把散射介质看作光波的多路复用器，通过空间光调制器和全场干涉测量获取透过矩阵，这种方法要求散射介质不随时间变化，只适用于静态散射介质。综上所述，目前针对时变散射介质的长距离光学成像问题，尚无有效的解决方法。

发明内容

本发明针对散射条件下的光学成像难题，提出一种基于大气散射介质调制的单像素成像方法，以准直激光束穿过大气散射介质，利用大气散射介质的自然随机特性和时间变化特性，形成一系列随机的瞬时调制样式照明观测目标，通过单像素探测器收集来自目标的光强，采用合适的算法重构目标图像，为强散射条件下的光学成像提供一种新的技术途径。

为实现上述目的，本发明提供了一种基于大气散射介质调制的单像素成像方法，所述方法包括：

激光束进入大气散射介质，以大气散射介质作为空间调制器，计算当前条件下的大气散射介质的单次调制样式；在该调制样式下，大气反射方向的光经过透镜到达目标，桶探测器收集所有来自目标的光信号，获得探测值的一个元素值；

随机改变激光束进入大气的传输方向，获得大气对光场调制的一系列不同的调制样式；以单次调制样式为行向量，将一系列不同的调制样式对应的行向量组成大气散射介质测量矩阵；并获得一系列不同的调制样式下对应的一系列探测值组成的测量向量；

基于测量向量和大气散射介质测量矩阵，利用图像重构算法恢复目标信号。

作为上述方法的一种改进，所述以大气散射介质作为空间调制器，计算当前条件下的大气散射介质的单次调制样式，具体为：

基于激光波长、激光传输方向、大气水平气象视距、大气分子和气溶胶光学厚度廓线信息，通过蒙特卡罗模拟方法计算大气散射介质的单次调制样式。

作为上述方法的一种改进，所述基于激光波长、激光传输方向、大气水平气象视距、大气分子和气溶胶光学厚度廓线信息，通过蒙特卡罗模拟方法计算大气散射介质的单次调制样式，具体包括：

设激光波长为λ，随机进入大气散射介质的方向相对于目标的天顶角和方位角为θ₀和

水平气象视距为v；当10⁶个光子被发射，光子距目标的垂直高度为z，对应的大气分子和气溶胶光学厚度分别为τ_m(z)和τ_a(z)，基于蒙特卡罗模拟方法，以4个介于[0,1]之间的随机数完整模拟光子在大气散射介质中的传输过程：

确定碰撞性质：

已知某一高度z上的分子光学厚度与气溶胶光学厚度为τ_m(z)和τ_a(z)，借助一个在[0,1]区间内均匀取值的随机数r₀，当0≤r₀≤τ_m(z)/[τ_m(z)+τ_a(z)]时，确认光子与分子发生碰撞；当τ_m(z)/[τ_m(z)+τ_a(z)]<r₀≤1时，确认光子与气溶胶发生碰撞；

确定光子自由路径：

对于一个随机碰撞过程，非碰撞概率r₁为一个在[0,1]间均匀取值的随机数，由公式l＝-ln(1-r₁)计算得到与r₁对应的光学距离值l；

确定碰撞后的方向：

分子的散射相函数采用Rayleigh散射相函数P(θ)：

其中，θ为散射角；

气溶胶粒子的散射相函数采用改进的Henyey-Greenstein相函数：

其中，g为不对称因子；

定义归一化的散射相函数p(θ,g)为：

则p(θ,g)的取值位于[0,1]之间；

利用两个在[0,1]间均匀取值的随机数r₂和r₃，确定碰撞后光子传输的散射角和方位角

r₂＝p(θ,g)

碰撞与传输过程重复进行，直至光子到达目标表面；以目标像元为中心，一定范围内均为分为若干间隔，记录每个间隔内接收的散射光子的数目，即获得当前传输条件下大气散射介质的调制样式。

作为上述方法的一种改进，所述基于测量向量和大气散射介质测量矩阵，利用图像重构算法恢复目标信号，具体包括：

当大气散射介质测量矩阵Φ满足RIP条件时，目标信号x通过下面的优化目标函数得到：

x＝Ψx'

其中，Ψ为正交基矩阵，x'是信号x在正交基矩阵上的投影系数向量；τ为常数因子，‖·‖_p代表lp范数，y为m维矢量的测量向量，由m个探测值组成，m也为调制次数，即大气散射介质测量矩阵Φ的行数为m。

作为上述方法的一种改进，所述图像重构算法采用基于l₁范数最小的BP算法、基于l₁范数最小的GPSR算法或者基于l₀范数最小的贪婪追踪系列算法。

本发明的优势在于：

1、云雾等大气强散射条件下的光学成像是亟待解决的科学和技术难题，本发明利用大气散射介质的自然随机特性和时间变化特性，作为单像素成像的空间调制器，实现被散射介质包围的目标的清晰成像。

2、在本发明的技术方案中，大气散射介质不再是观测目标与探测器之间的“干扰”，而是成像系统本身的理想的组成部分，从而将云雾等大气强散射介质对光学成像的不利影响转换成有利条件，提升强散射介质条件下的光学成像能力。

附图说明

图1为现有技术的单像素成像典型光路示意图；

图2为现有技术的压缩感知数学模型；

图3为本发明的基于散射介质调制的单像素成像示意图；

图4为本发明的大气散射介质测量矩阵生成示意图；

图5为目标图像(32×32)；

图6为目标图像的重构图像。

具体实施方式

下面结合附图给出本发明的实施流程，以详细说明本发明的技术方案。

1、基本概念

单像素成像的典型光路如图1所示。光源发出的光经过透镜扩束后到达空间调制器(Digital Micromirror Device，DMD)，DMD由可独立编程的微透镜阵列组成，每个微透镜可向两个方向翻转，即可做二值调制，将随机0/1测量矩阵加载于DMD上，即可实现对入射光场的二值随机调制。来自于DMD一个反射方向的调制光经透镜到达目标，经目标透射后经过透镜由桶探测器全部收集探测，即可得到探测值y的一个元素值。

测量过程可表示为：

y＝Φx+e (1)

其中信号x为n维矢量，测量值y为m维矢量，Φ为DMD调制获得的m×n维测量矩阵，e为噪声。根据压缩感知(Compressed Sensing，CS)理论，大部分的信号都是稀疏的或是可压缩的，根据信号的稀疏性，只需要远低于奈奎斯特率的少量采样，就可以将信号完美恢复。CS数学模型如图2所示，信号x通过CS算法稀疏表征：

x＝Ψx' (2)

其中x’是稀疏的。由此，测量值可表示为：

y＝ΦΨx'+e (3)

在基于CS理论的单像素计算成像中，空间光调制器DMD是核心器件，作用是生成满足单像素成像要求的测量矩阵。理论证明随机测量矩阵是较好的测量矩阵，但在实际应用中却面临硬件实现的困难，经典的DMD也只能实现部分随机测量矩阵。考虑到大气散射介质具有天然的随机特性，而且这种特性随时间变化，有可能为单像素成像提供一种纯物理的压缩感知“器件”。

基于散射介质调制的单像素成像如图3所示，成像过程与单像素成像典型光路相似，不同之处在于以大气强散射介质作为空间调制器，代替传统的DMD，激光穿越大气散射介质时，由于大气散射介质的随机特性和时变特性，介质对光场的调制也是随机且随时间变化的，一系列不同的调制样式经过透镜照明目标。激光在大气散射介质中传输的距离越长，受到的调制作用越强。以往散射成像难以克服的长距离问题，反而成为成像的有利条件。同时，为了进一步确保散射介质调制作用的随机性，不断改变光源的位置，使激光以随机的角度进入大气散射介质。针对每一次调制，桶探测器收集所有来自目标的光信号，获得一系列探测值，最后利用图像重构算法恢复目标信息。

2、测量矩阵

在基于CS的单像素成像中，测量矩阵是人为构造的。利用大气散射介质进行单像素成像，也需要获取每次测量所对应的调制样式，所有的调制样式组成最终的测量矩阵。

基于激光波长、激光传输方向、大气水平气象视距、大气分子和气溶胶光学厚度廓线等信息，通过蒙特卡罗(Monte Carlo，MC)模拟方法计算大气散射介质的单次调制样式。假设激光波长为λ，随机进入大气散射介质的方向相对于目标的天顶角和方位角为

水平气象视距为v。假设10⁶个光子被发射，光子距目标的垂直高度为z，对应的大气分子和气溶胶光学厚度分别为τ_m(z)和τ_a(z)，基于MC方法，以4个介于[0,1]之间的随机数完整模拟光子在大气散射介质中的传输过程。

(1)确定碰撞性质

大气中存在分子和气溶胶两种粒子，碰撞发生时，首先要判断光子与哪种粒子碰撞。已知某一高度z上的分子光学厚度与气溶胶光学厚度为τ_m(z)和τ_a(z)，借助一个在[0,1]区间内均匀取值的随机数r₀，当0≤r₀≤τ_m(z)/[τ_m(z)+τ_a(z)]时，确认光子与分子发生碰撞；当τ_m(z)/[τ_m(z)+τ_a(z)]<r₀≤1时，确认光子与气溶胶发生碰撞。

(2)确定光子自由路径

一般用光学距离l来表示连续两次碰撞之间光子所经历的自由路程。如取非碰撞概率密度函数p与l间的函数关系为p＝e^-l，那么光子与粒子相距光学距离l时的非碰撞概率r₁应为：

对于一个随机碰撞过程，r₁是一个在[0,1]间均匀取值的随机数，由公式l＝-ln(1-r₁)便可算得与r₁对应的光学距离值l。

(3)确定碰撞后的方向

散射相函数P(θ,g)(θ为散射角，g为不对称因子)决定光子碰撞后的去向，它由碰撞的性质决定。分子散射采用Rayleigh散射相函数：

对于气溶胶粒子，采用改进的Henyey-Greenstein相函数：

定义归一化的散射相函数：

则p(θ,g)的取值位于[0,1]之间。利用两个在[0,1]间均匀取值的随机数r₂和r₃，可以确定碰撞后光子传输的散射角和方位角。

(4)获得特定传输条件下的调制样式

碰撞与传输过程重复进行，直至光子到达目标表面。以目标像元为中心，一定范围内均为分为若干间隔，记录每个间隔内接收的散射光子的数目，即可以获得当前传输条件下大气散射介质的调制样式。

(5)测量矩阵生成

随机改变激光束进入大气的传输方向，获得大气对光场调制的一系列调制样式，以单次调制样式为行向量，组成大气散射介质测量矩阵，具体过程如图4所示。图中，调制次数即为测量次数，决定测量矩阵的行数。

3、重构算法

当测量矩阵Φ满足RIP(restricted isometry property)条件时，信号x可通过下面的优化目标函数得到：

其中τ为一常数因子，‖·‖_p代表lp范数。

常用的方法包括基于l₁范数最小的BP(Basis Pursuit)、GPSR(GradientProjection for Sparse Reconstruction)算法，基于l₀范数最小的贪婪追踪系列算法，如OMP(Orthogonal Matching Pursuit)、CoSaMP(Compressiv Sampling Matching Pursuit)等算法。

通过仿真与实验相结合的方法展示本申请的技术效果。

1、大气散射介质瞬时调制样式模拟

首先通过数值方法模拟激光在大气中传输时大气散射介质的瞬时调制样式。大气环境的构造采用MODTRAN辐射传输模型。该模型将0-100km的大气分为34层，对对每一层的物理参数，从气压、温度到分子吸收系数、消光系数等，都进行了定义。设定激光波长为532nm，采用1976US Standard Atmosphere模式和RURAL气溶胶消光模式。由于本技术利用大气的时变散射特性进行成像，要求大气处于浑浊状态，且随时间变化，因此设置水平气象视距在[0.5km,2km]范围内随机取值。为验证本技术适用于长距离成像，设置激光穿越整层大气，因此传输距离为100km。为确保大气调制的随机性，我们设计激光进入大气的方向是随机的，因此设置激光传输方向的天顶角在[150°,180°]范围内随机取值，方位角在[0°,360°]范围内随机取值。通过MC方法，获得一系列大气散射介质瞬时调制样式。

2、成像实验

搭建了基于大气散射介质调制的单像素成像实验光路。实验中成像目标如图5所示，图像分辨率为32×32。根据实验光路图，半导体激光器发出的波长532nm的激光经过10倍扩束器扩展后投射到DMD上。该DMD分辨率为1920×1080。注入DMD的图样为模拟生成的大气散射介质瞬时调制样式，像素为32×32。DMD上每32×32个微透镜组成调制样式的一个像素。大气调制后的结构光照明目标，经目标反射后由透镜组成像在光电倍增管(photomultiplier tube，PMT)的光敏面上，由数据采集卡进行模数转换和量化。最后由定制程序对DMD和采集卡进行控制。根据目标图像尺寸，为保证全采样，注入DMD的大气散射介质瞬时调制样式共32×32＝1024个，PMT对应获得1024个测量值(测量次数也可以选择为32×32的一半)。

3、图像重构

基于1024个测量值，采用CoSaMP算法进行图像重构。CoSaMP算法的输入、输出及流程如下：

输入：测量矩阵A，观测值y，迭代次数为t＝1。

输出：重构信号估计X。

算法流程：

①初始化残差r₀＝y，索引集

②计算相关系数|<r_t-1,a_j>|，取其绝对值后降序排列，将其前2K项构成集合

对应位置构成集合J₀；

③更新索引集和新的测量矩阵：Λ_t＝Λ_t-1∪J₀,

④求解最小二乘解：X_t＝argmin_x||y-A_tx||₂；

⑤从X_t中取绝对值最大的前K项，记为X_tK,更新集合J_tK；

⑥更新残差：

达到迭代次数停止，否则返回②。

图像重构结果如图6所示。左侧为重构的二维图像，右侧为重构结果的三维展示。通过大气散射介质对光波的瞬时调制，利用单像素探测器进行信号收集，通过图像重构方法成功实现目标成像，验证了本方法的可行性。

最后所应说明的是，以上实施例仅用以说明本发明的技术方案而非限制。尽管参照实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，对本发明的技术方案进行修改或者等同替换，都不脱离本发明技术方案的精神和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (5)

1.一种基于大气散射介质调制的单像素成像方法，所述方法包括：

激光束进入大气散射介质，以大气散射介质作为空间调制器，计算当前条件下的大气散射介质的单次调制样式；在该调制样式下，大气反射方向的光经过透镜到达目标，桶探测器收集所有来自目标的光信号，获得探测值的一个元素值；

随机改变激光束进入大气的传输方向，获得大气对光场调制的一系列不同的调制样式；以单次调制样式为行向量，将一系列不同的调制样式对应的行向量组成大气散射介质测量矩阵；并获得一系列不同的调制样式下对应的一系列探测值组成的测量向量；

基于测量向量和大气散射介质测量矩阵，利用图像重构算法恢复目标信号。

2.根据权利要求1所述的基于大气散射介质调制的单像素成像方法，其特征在于，所述以大气散射介质作为空间调制器，计算当前条件下的大气散射介质的单次调制样式，具体为：

基于激光波长、激光传输方向、大气水平气象视距、大气分子和气溶胶光学厚度廓线信息，通过蒙特卡罗模拟方法计算大气散射介质的单次调制样式。

3.根据权利要求2所述的基于大气散射介质调制的单像素成像方法，其特征在于，所述基于激光波长、激光传输方向、大气水平气象视距、大气分子和气溶胶光学厚度廓线信息，通过蒙特卡罗模拟方法计算大气散射介质的单次调制样式，具体包括：

设激光波长为λ，随机进入大气散射介质的方向相对于目标的天顶角和方位角为θ₀和φ₀，水平气象视距为v；当10⁶个光子被发射，光子距目标的垂直高度为z，对应的大气分子和气溶胶光学厚度分别为τ_m(z)和τ_a(z)，基于蒙特卡罗模拟方法，以4个介于[0,1]之间的随机数完整模拟光子在大气散射介质中的传输过程：

确定碰撞性质：

已知某一高度z上的分子光学厚度与气溶胶光学厚度为τ_m(z)和τ_a(z)，借助一个在[0,1]区间内均匀取值的随机数r₀，当0≤r₀≤τ_m(z)/[τ_m(z)+τ_a(z)]时，确认光子与分子发生碰撞；当τ_m(z)/[τ_m(z)+τ_a(z)]<r₀≤1时，确认光子与气溶胶发生碰撞；

确定光子自由路径：

对于一个随机碰撞过程，非碰撞概率r₁为一个在[0,1]间均匀取值的随机数，由公式l＝-ln(1-r₁)计算得到与r₁对应的光学距离值l；

确定碰撞后的方向：

分子的散射相函数采用Rayleigh散射相函数P(θ)：

其中，θ为散射角；

气溶胶粒子的散射相函数采用改进的Henyey-Greenstein相函数：

其中，g为不对称因子；

定义归一化的散射相函数p(θ,g)为：

则p(θ,g)的取值位于[0,1]之间；

利用两个在[0,1]间均匀取值的随机数r₂和r₃，确定碰撞后光子传输的散射角和方位角

r₂＝p(θ,g)

碰撞与传输过程重复进行，直至光子到达目标表面；以目标像元为中心，一定范围内均为分为若干间隔，记录每个间隔内接收的散射光子的数目，即获得当前传输条件下大气散射介质的调制样式。

4.根据权利要求3所述的基于大气散射介质调制的单像素成像方法，其特征在于，所述基于测量向量和大气散射介质测量矩阵，利用图像重构算法恢复目标信号，具体包括：

当大气散射介质测量矩阵Φ满足RIP条件时，目标信号x通过下面的优化目标函数得到：

x＝Ψx'

其中，Ψ为正交基矩阵，x'是信号x在正交基矩阵上的投影系数向量；τ为常数因子，‖·‖_p代表lp范数，y为m维矢量的测量向量，由m个探测值组成，m也为调制次数，即大气散射介质测量矩阵Φ的行数为m。

5.根据权利要求4所述的基于大气散射介质调制的单像素成像方法，其特征在于，所述图像重构算法采用基于l₁范数最小的BP算法、基于l₁范数最小的GPSR算法或者基于l₀范数最小的贪婪追踪系列算法。

Images