CN114037625A

CN114037625A - 一种基于物理散射模型的图像恢复方法

Info

Publication number: CN114037625A
Application number: CN202111261649.5A
Authority: CN
Inventors: 田建东; 周士筠; 张箴
Original assignee: Shenyang Institute of Automation of CAS
Current assignee: Shenyang Institute of Automation of CAS
Priority date: 2021-10-28
Filing date: 2021-10-28
Publication date: 2022-02-11
Also published as: WO2023070958A1; US20240046421A1

Abstract

本发明涉及一种基于真实物理参数散射物理模型的图像恢复方法。建立了一种不同于通用大气光散射模型的新的散射模型，能对每一个散射环境给出相对应的物理解释，阐述了光源位置和散射情况的关系，解决了非均匀光照下的散射去除问题。首先使用积分球来去除相机暗角，然后通过改进的类暗通道先验方法找出纯散射部分，最后根据散射介质的光衰减和后向散射分布计算天顶角、方位角、消光系数和其他相关物理量来消除后向散射和光的衰减。理论上，该模型可以用于任何均匀散射介质，如水下和雾环境，并且在高密度介质中工作良好。最终的实验结果表明，该模型解决了强散射介质下的后向散射去除问题，是正确的、有效的，并且具有实际的物理意义。

Description

一种基于物理散射模型的图像恢复方法

技术领域

本发明涉及一种散射模型，具体说是一种基于物理参数的图像恢复方法，解释了光源位置和后向散射的关系。

背景技术

近年来，科技的进步允许人们在强散射环境下采集图像和视频。然而在实际的工程环境中，去除后向散射一直是一个持续存在的问题。后向散射导致图像颜色失真和模糊等情况。这严重影响了标准计算机视觉算法的性能，因为一般来讲，这些算法是为在晴朗的空气中工作而开发的。因此，为了使用这些算法获得合理的结果，必须首先将图像消除散射。而著名的大气散射模型基于Beer-Lambert法则：

I(x)＝J(x)t(x)+A(1-t(x)) (1)

t(x)＝e^-σz (2)

其中，大气光A取决于深度图，t(x)与衰减系数σ和深度z有关，输入图像为I(x)，清晰图像为J(x)。该模型较简单，只需使用几个参数，即可去除空气和水中的散射。这也带来一些缺点，即不能涵盖全部复杂的情况。

当午间光从正南照射过来时，天顶角接近一天中的最小值，此时图像整体亮度相较于其他时间点更亮。天顶角是光线入射方向与天顶方向的夹角，其变化周期可以为一整年。向光源拍摄时，照片上的图像亮度是关于中间轴两面对称的，同理而言，拍摄角度与光照角度为90°时，图像亮度是关于图像的长轴单调的。一张图像中，离光源越远的部分亮度越低，同理，镜头离光源越远，图像整体亮度越低。在水下时，这体现为水面离镜头的距离。相机的设置，比如镜头光圈大小，焦距，快门速度等也会影响图像整体的亮度以及像素分布情况。这些由光引起的现象都会影响后向散射的去除。所以将这些因素全部考虑是必要的，且物理的。

由上可得，光源位置的变化会引起图像像素值的变化，以及后向散射分布的变化，而大气散射模型不能考虑光源的信息。考虑光源位置和散射之间的关系是本发明提出的方法的核心思想。

发明内容

针对上述技术不足，本发明提供了一种基于物理参数散射模型的图像恢复方法。1)提出了一种基于物理参数的新模型，用于消除散射介质中拍摄图像的后向散射和光的衰减带来的图像不清晰问题。2)发现了著名的大气光衰减模型的局限性，并通过提出一个与其完全不同的模型方法来解决非均匀光照造成的图像不清晰的问题。解释了大气散射环境和光源位置的关系，理论上，该模型可用于任何均匀散射介质(例如水下和烟雾环境)中拍摄的图像。3)该模型在解决带有高密度散射问题的图像时尤其有效。

本发明解决其技术问题所采用的技术方案是：基于物理散射模型的图像恢复方法，建立关于光源物理位置与散射之间的模型映射关系，将衰减图像分为后向散射分量和直接辐射分量计算；将图像采样取点，通过多次拟合逼近最真实的物理参数；将计算出的物理参数反代入模型中，以计算整幅图像的衰减值和后向散射，从而恢复清晰图像。

该方法包括以下步骤：

步骤1.对衰减图像利用积分球校准相机的渐晕；

步骤2.将图像像素坐标转换为世界坐标系下的坐标；

步骤3.建立去散射数学模型，将衰减后的图像分为直接辐射和后向散射两部分，分别推导模型映射关系；

步骤4.对去除渐晕的校正图像进行采样，OSTU二值化处理，获取纯散射的部分像素区域，用于拟合和参数计算；

步骤5.使用计算得到的参数来去除整幅图像的后向散射和衰减。

所述使用积分球校准相机的渐晕，包括：

建立渐晕表；

使用积分球拍摄原始衰减图像，获取每个像素位置上的校正系数，并利用该系数对原始衰减图像逐像素使用查找表方法来消除渐晕，获取校正图像。

所述直接辐射项可以用下式表示：

所述后向散射项可以用下式表示：

所述清晰图像用下式表示：

其中，I^c(x)为输入的待恢复图像，L_b ^c(x,y)为后向散射项，

为直接辐射项。

所述获取纯散射的部分像素区域为提取只包含散射分量且具有相似深度的像素用于估计物理参数：

a.将输入图像由RGB空间转换至HSV空间，提取其饱和度通道和强度通道；

b.利用OSTU最大类间方差法计算最佳全局阈值，获取高强度和饱和度的图像区域；

c.进一步在图像的x和y方向上间隔采样，生成点集V；

d.在图像点集V内应用改进的类DCP法则：

即得到J^c(x)≈0的图像区域。

所述物理参数包括太阳入射的天顶角、方向角、相机焦距、光照度、消光系数、物体到相机的深度、各向异性常数、水中物体到水面的距离。

在拟合和参数计算时，由于不同的参数对函数的影响有区别，所以根据参数的拟合优先级次序依次进行拟合求解。

该方法根据光源位置和散射分布的关系，对非均匀光照下的图像进行恢复。

本发明具有以下有益效果及优点：

1.本发明方法建立了一种新的图像去散射模型方法，与当下流行的图像增强方法不同，该模型根据图像上的后向散射情况，反推拍摄场景的物理参数，对去除后向散射的过程做出了物理解释，通过物理方法的逆过程求解来解决问题，使得图像恢复效果更自然。

2.本发明方法创新性地解释了光源位置和后向散射的关系，解决了非均匀光照下的图像去散射问题。

3.本发明方法在对去散射方面要优于现有的算法，且表现稳定，尤其在图像出现强散射时效果更为突出。

附图说明

图1是物理参数散射物理模型的图像恢复方法流程示意图；

图2是纯散射无物体情况下拍摄的图像；

图3是非均匀光照的原因；

图4是光源的天顶角和方位角；

图5是后向散射的相关计算参数；

图6是实验环境；

图7是关于θ_A变化的像素值变化情况；

图8是关于θ_Z变化的像素值变化情况；

图9是关于g,σ变化的像素值变化情况；

图10是关于f变化的像素值变化情况；

图11是关于Z_t,M,Y_u变化的像素值变化情况；

图12是关于表2对应的高散射图像的恢复图像；

图13是本发明与其他方法的对比实验结果。

具体实施方式

下面结合实施例对本发明做进一步的详细说明。结合附图对方法步骤进行说明。

一种基于物理参数的散射模型，包括以下步骤：

1.使用均匀光源(如积分球)校准相机的渐晕及去除图像噪声；

2.建立去散射数学模型，将衰减后的图像分为直接辐射和后向散射两部分分别进行计算；

3.对采集的图片进行OSTU二值化处理，对纯散射的部分像素拟合和参数计算。

4.使用计算得到的参数来去除整幅图像的后向散射和衰减。

所述校准渐晕包括：

建立包括相机所有焦距等设置相对应的去渐晕表；对待处理的图像使用查表法，使用对应的去渐晕表去除渐晕。渐晕是一种位置相关的照明衰减现象，对于所有相机镜头来说都是不可避免的。通常，当从光轴上的对象点穿过光圈的光束大于离轴对象点的成像光束时，会导致渐晕。就渐晕消除而言，理想的第四余弦定律在实践中无法验证，因此有许多类似曲线的近似工作，以及对第四余弦定律的改进工作。在本文中，我们使用查找表(LUT)方法来消除渐晕。此操作应在暗室中进行，设置具有低镜面反射和已知输入强度的均匀光源。我们使用直径为1000mm的积分球，内置12盏灯，出光口直径为100mm，确保光源均匀。每个像素位置上的校正系数为：

其中，I_ref(i,j)表示位置(i,j)处的像素值，I_ref,max表示I_ref(i,j)的最大值，I_LUT(i,j)表示坐标(i,j)处的校正因子。计算所有位置的校正因子，并将其组合到去渐晕表中。此后，使用相同相机设置拍摄的图像可以使用查表法进行校正：

I′(i,j)＝I(i,j)I_LUT(i,j) (4)

其中I(i,j)表示像素(i,j)处的输入强度值，I′(i,j)表示该像素渐晕校正后的输出强度值。

所述后向散射为光不经过物体反射，直接通过散射介质中的悬浮粒子发生随机散射到达相机镜头引起的光学效应。所述直接辐射分量为光经过物体表面反射，再通过散射介质中的悬浮粒子发生散射到达相机镜头引起的光学效应。

所述去散射数学模型包括：

1.非均匀光照概述

图2为两幅在水下拍摄的纯散射无物体图像(a)、(e)，以大气光散射而论，它被标记出的两条线的曲线形状应该趋向为常值。然而如图2所示，被标出的像素值的曲线形状与光源位置有很大关系。对于如(b)、(f)中标记的垂直线(d)、(h)，其像素值随着纵坐标的减小而减小，这是因为光源距离纵坐标较大的位置较远。对于如(b)、(f)中标记的水平线(c)、(g)，条件更加复杂。我们发现，当地球平面上的拍摄方向和光源方向之间的角度为kπ(k＝1,2,3…)时，其曲线形状保持关于中轴对称。当地球平面上的拍摄方向和光源方向之间的角度为kπ+π/2(k＝1,2,3…)时，其曲线形状保持单调递增或单调递减。当拍摄角度介于上述两种情况之间时，其曲线的像素值将有极大值或极小值。上述问题可以简单地概括为离光源越近，照度越高，像素值越大。此外，我们还发现，由于光在密度更大的介质中衰减更快，随着介质浊度的增加，这种现象变得更严重。这些问题是不可忽视的，在稠密介质中进行深度估计和图像恢复时会出现问题。

如图3所示，光传播导致照明和强度不均匀。图3(a)表明了摄像机接收到的光由两部分组成：物体的直接衰减部分和光源的散射部分。然而，广泛使用的大气光散射模型没有考虑从光源到物体的光传输。如图3(b)所示，光D₂比光D₁传播更长的距离，这意味着光线D₂衰减得更厉害。大气光散射模型还假设散射光是平行的，并且传播到散射点的距离相等。然而，如图3(c)所示，在不同的摄像机视角(O1和O2)上，从光源到散射点(L1和L2)的散射光线也不同。图3(b)和图3(c)所体现的两个现实因素导致了照明其实并不均匀。

2.模型建立

将摄像机的电荷耦合器件(CCD)作为世界坐标系的原点。令图像平面上的x轴和y轴平行于世界的X轴和Y轴，Z轴与摄像机的光轴对齐。对于图像平面上的点x＝(x₀,y₀)，令其对应物体表面上的点X＝(X₀,Y₀,Z₀)，他们之间的关系为

其中f表示相机的焦距。

那么对于一个图像坐标系上的点(x,y),其对应在世界坐标系上的坐标可以表示为：

其中x₀和y₀是图像的中心坐标，w_x和w_y是根据CCD宽度计算的像素的物理尺寸，f为焦距大小。

采用单散射模型，忽略前向散射，因此到达相机的辐射量可以被分为直接辐射和后向散射：

其中，I_b为后向散射，I_d为直接辐射分量。因此，精确计算每个像素的后向散射强度和衰减，即可得到清晰的图像。

根据比尔-朗伯定律，光被认为随着从光出射点到透镜的距离而衰减，光的传播与消光系数有关。其中，消光系数基于光波长而变化。在水中，可见光在最长的波长处衰减，最终在人眼中呈现蓝色，而红色往往会比绿色和蓝色衰减得更快，所以不同通道的消光系数不应该是相同的。此外，当太阳光从远处照射来时，向每个方向的散射的部分也应当列入考虑，这由相位函数和散射系数所决定。

3.直接辐射分量

对于图像平面上的点x＝(x₀,y₀)，假设其对应于世界坐标系内物体上的坐标X＝(X₀,Y₀,Z₀)。那么对于该点来说，直接辐射分量为：

其中,‖D(X)‖表示光传播到物体的距离，I₀表示进入散射介质前光源的辐射强度，c表示颜色通道。距离‖D(X)‖包括水中的距离‖τD(X)‖(0≤τ≤1)和空气中的距离‖1-τD(X)‖。此后，光经过物体反射，反射率为ρ^c(X)，并经过距离‖X‖到达相机镜头。经过所有过程后，到达相机镜头的直接辐射度可以写为：

I_d(x)＝E^ce^{-σ(‖τD(X)‖+‖X‖)} (8)

令Y_u为相机到水平面的距离。如图5所示，考虑到图4所示的天顶角θ_Z，光在水中被物体反射之前的传播距离等于：

令r(x,y)等于物体到相机镜头的距离‖X‖,那么得到：

最终可用物理参数表示直接辐射分量：

其中，J^c(x)表示清晰图像，r表示r(x,y)像素点到原点的距离。

当光源在散射介质中时，Y_u指的是物体到光源的距离；否则Y_u指的是物体到水平面的距离或者到大气层顶部的距离。

4.后向散射分量

后向散射是在光不经过物体表面反射，只经过介质散射到达相机镜头，由水中悬浮颗粒的小角度随机散射引起的光学响应。使用相位函数来描述光在特定方向上的散射程度，并使用散射系数β来表示光散射的程度。在相函数方面，最经典的是Heyney-Greenstein相函数，它可以用分子各向异性g和散射角α两个参数来表示,通过g∈(-1,1)来描述前向和后向散射的相对数量:

其中，散射角α取值范围为[0,π],为物点方向与光源方向的夹角：

u＝[sin(θ_Z) sin(θ_A) cos(θ_Z) sin(θ_Z) cos(θ_A)] (16)

其中，θ_Z表示天顶角，θ_A表示方位角。

根据原理u·X＝‖u‖‖X‖cosα，散射角可以表示为现实的物理量：

由式(10)得，

后向散射是通过对从相机到散射点深度的积分得到的：

其中，Z_s表示深度，β^c表示消光系数，σ^c表示散射系数。

化简式子，令

M^c＝E^cβ^c,得到后向散射分量如下，由焦距、散射角、深度、消光系数等组成。

所述像素拟合和参数计算包括：

根据式(6)(12)(20)，可得到清晰图像的公式为：

其中，I^c(x)为输入的待恢复图像，L_b ^c(x,y)为后向散射项，

为直接辐射项。

在去除图像中的渐晕后，将图像平面坐标转换到真实坐标系中进行去散射的计算。对图像进行步长为50像素的采样，对采样点采用类似于暗通道先验的方法来寻找直接辐射度接近0的部分，设为集合V。对集合V,采用优化拟合方法计算中公式中未知参数M^c,σ,Y_u,θ_z,f,θ_z,θ_A,Z_s。最后通过得到的参数，计算整幅图像上的后向散射部分，并减去后向散射项，除以衰减项可得到最终结果。

1.实验准备

为了验证模型的准确性，可对实验环境布置如图6所示。我们使用在充满散射介质(稀释牛奶)的水缸中进行实验。光源为太阳，使用一个尺寸为90×90×45厘米的水箱，用清水注满。往里添加牛奶，多次增加散射环境的浊度。我们设计了实验来测试和验证我们的图像恢复模型。由于方程式与实际物理值相关，可做如下设置，避免客观因素对物理参数计算的干扰：

除了放置镜头那一面外，使用黑色遮光布挡住水缸的里外所有面，包括水缸底面，防止阳光在玻璃间反射造成计算干扰；

使用指南针校准水缸方向，使得水缸的南北向与地轴的南北极平行；在水平的屋顶放置水缸，防止天顶角计算不准确；

放入牛奶后搅拌均匀，静置等水面静止后再进行拍摄；

关闭相机的自动补光、自动调校、周边光量矫正等功能，使用最原始的raw数据进行处理；。

拍摄照片时，尽量去除或避开吸附在拍摄这一侧的玻璃表面的水滴；

拍摄照片时将相机水平握持，尽量避免倾斜。否则可能导致方向角等参数计算不准确。

2.实验过程

2.1软件环境配置

本次实验的计算机平台的配置为CPU处理器Intel(R)Core(TM)i9-10900K CPU@3.70GHz，缓存32G，显卡Nvidia GeForce RTX2080 SUPER。软件环境为Windows 10专业版。实验的软件开发环境为Matlab2020b。

2.2参数计算

为了计算后向散射分量I_b，需要找到一些区域，其像素值与后向散射值相近，而直接辐射分量等于0。

所以需要提取一些只包含散射分量且具有相似深度的像素来估计物理参数。为此，首先将输入图像由RGB空间转换至HSV空间，提取其饱和度通道和强度通道。然后，使用Matlab中的graythresh函数，即OSTU最大类间方差法，计算出最佳全局阈值，提取出一些强度和饱和度都很高的区域。在这些区域中，进一步在图像的x和y方向上以50像素的间隔采样，以生成点集V。对于属于V的每个点，取其五个像素相邻区域的最小值如下。

公式(24)受He等人提出的暗通道先验法的启发。不同之处在于，本发明没有提取不同颜色通道的最小值，因为这样计算可以适配更多光源类型。此外，为了确保用于计算的点只包含后向散射项，即J^c(x)≈0，我们在点集V上应用了此改进的类DCP法则，而不是在整个图像上应用。

2.3曲线拟合

使用提取出的只包含后向散射的区域，来计算对应的计算物理参数。对于我们在公式(20)中提出的后向散射项，它包含八个参数,M^c,σ^c,Y_u,g,f,θ_z,θ_A,Z_t。理论上，对于八个参数的方程，取大于等于十个点就可以计算出所有参数的解析解。但是逐一求解比较困难耗时，而且取点受噪声影响较大，拟合是比较省时的方法。由于水下扰动的因素太多导致大量的噪音，还有参数太多的原因，方程的解很容易在拟合过程中发生偏移，对参数的优化是必要的。在拟合计算中，每个参数有三个对应的值，分别为初始值，拟合上限，拟合下限。选对初始点尤为重要。

在实验中，我们使用牛奶作为散射介质。根据以前工作中牛奶各向异性参数的计算，我们可以大致将g的值范围从[-1,1]减少到[0,1]。由于滴下的牛奶量远小于水箱的体积，因此g的初始值可以设置为较小的值0.1。如图所示，g主要影响三维坐标中沿纵轴的坡度，这体现在亮度差上，随水深的增加而逐渐变化。

表1

使用控制变量法，可知以表1所示的参数带入后向散射函数得到的数学关系。θ_A表示从光源到相机的方位角，它主要控制像素值在x轴方向上的倾斜。如图7(a)(b)(c)所示，当两个函数的θ_A关于180°对称时，他们的像素值z关于x＝x_s/2对称，倾斜方向在[0,π]和[0,π]之间也各不同。其中，x_s为图像坐标系内的最大值。数学推导如下。

定义函数

则

f(-x)＝f(2π-θ_A) (26)

而相位函数

显然

P(-x)＝P(2π-θ_A) (28)

而在式(20)中，相位函数是唯一包含x和θ_A的部分，因此关于后向散射的函数

L_b(-x)＝L_b(2π-θ_A) (29)

因为公式(20)是以图像中心为坐标原点，而计算时同常以图像中的第一个像素点为坐标原点，所以上式可以更改为

L_b(x_s-x)＝L_b(2π-θ_A) (30)

另外，如图7(d)(e)(f)所示，当θ_A＝kπ时，整体图像像素值关于x＝x_s/2对称。所以理论上来说，对同一张图像，θ_A只有一个解。在拟合时，应当将θ_A取值在0°到360°之间，设置其初始值以10°度为步长，遍历10°～350°度，下限以10度为步长，遍历0°～340°度，上限以10°度为步长，遍历20°～360°度。比较遍历过程中所有的拟合值，确定出θ_A在的区间及其初始值。

θ_Z表示光源的天顶角。如图8所示，θ_Z的变化主要控制函数振幅的变化。但是，当θ_Z增加时，式(20)中分数的分子和分母都会增加，因此与θ_Z和整体像素值没有定性关系。在拟合时，θ_Z取值在0°到90°之间，设置其初始值以5°为步长，遍历5°～85°，下限以5°为步长，遍历0°～80°，上限以5°为步长，遍历10°～90°。比较遍历过程中所有的拟合值，确定出θ_Z在的区间及其初始值。

参数σ^c，g与散射介质本身有关。Friscad等提出了计算参与介质散射特性的理论模型，Narasimhan等通过增加稀释介质的浓度测量了40种常见材料的散射参数，包括g。如图9(a)所示，g主要影响y-z坐标平面上沿z轴的斜率。当拟合时，设置遍历g的步长为0.1，遍历[0,1]寻找最合适的初始值。如图9(b)所示,σ^c表示衰减系数，其随函数整体振幅变化，几乎不影响L_b对x,y,z的梯度。根据Narasimhan及Jerlov的研究得知该值一般为10^-5/mm到10^-3/mm数量级不等，故将该值的初始值设为10^-4数量级即可，上下限分别为10^-5和1，遍历步长为0.0001。

焦距f主要控制函数x轴和y轴上的曲率度，如图10所示，其表现为对y或z的梯度大小，以及后向散射最大值和最小值之间的差值。这可以通过相机拍摄的机制来解释，当相机位置不动的时候，焦距越大，视野越窄，而亮度呈曲线变化，所以焦距变化时会造成图像内亮度差距的变化。对于相机Cannon 5D Mark2而言，它在所有拟合计算中都可以直接取上下限为[24，105]，步长为10，进行拟合。此外，可以考虑直接从照片的EXIF标签中读取。

如图11所示，参数Z_t,Y_u,M主要影响整体振幅，很少改变斜率。Z_t表示从物体到照相机的深度。Z_t越大，图像整体越浑浊，后向散射程度越深。M通常不是一个大的值，它与f的乘积与背景光的值相近。此外，需要控制上下限距离足够大，防止无法取到所有高拟合度的解。

在拟合时，由于不同的参数对函数的影响有区别，所以应当设置拟合优先级。对于参数M^c,Y_u,σ^c,f,g,θ_A,θ_Z,Z_t分别设置优先级为[5 8 1 4 7 6 2 3]，按优先级次序进行拟合，选取拟合程度最大的一个结果。

整体的算法流程在算法1中给出。

3拟合结果与对比实验

对于如图12(a)的高散射下图像，肉眼几乎看不见里面物体的内容。算出它对应的现实物理参数如表2：

表2

它的恢复图像如图12(b)，可以清楚地看出本发明的效果去除了相当浓度的散射，并且估算的物理参数也非常接近真实值。这可以验证本模型的正确性。

图13展示了本方法与其他方法的对比实验结果，对比的方法分别为peng的方法(IBLA)，huang的方法(RGHS)，UDCP方法，Song的方法(ULAP)。肉眼可见本发明在提升图像视觉方面的优越性，比其他的方法都更有效。

4小结

本发明提出了一种新的基于物理的方法来恢复密集散射图像。这种方法允许特色在于，能够优越地处理非均匀照明，这是在常用的大气光散射模型中被忽略的。本发明使用完全真实的物理参数来描述后向散射和直接辐射的分布，利用曲线拟合的思想去除后向散射。相信这项工作将为混浊介质中的去散射、水下增强和图像恢复等研究开辟一条新的途径。