具体实施方式

在不同的附图中，相同的特征由相同的附图标记表示。特别地，各种实施例之间共有的结构和/或功能特征可以具有相同的附图标记，并且可以设置相同的结构、尺寸和材料特性。

除非另有说明，当对连接在一起的两个元件进行引用时，这表示除了导体之外没有任何中间元件的直接连接，而当对链接或耦接在一起的两个元件进行引用时，这表示这两个元件可以连接或者它们可以经由一个或多个其他元件链接或耦接。

在以下公开内容中，除非另有说明，当提及绝对位置限定符(诸如术语“前部”、“后部”、“顶部”、“底部”、“左部”、“右部”等)或提及相对位置限定符(诸如术语“上方”、“下方”、“更高”、“更低”等)、或提及取向限定符(诸如“水平”、“竖直”等)，参考图中示出的取向或者参考在正常使用期间所取向的微测辐射热计。

除非另有说明，否则表述“大约”、“近似”、“基本上”和“在……的量级”表示在10％以内，优选地在5％以内。

图1示出了包括对IR光敏感的像素阵列102的IR成像设备100。IR成像设备100例如形成IR相机的一部分。在一些实施例中，像素阵列102对长波IR光敏感，诸如具有在7至16μm或更高范围内的波长的光。

在图1的示例中，像素阵列102包括以12行12列布置的144个微测辐射热计像素104。在替代性实施例中，像素阵列102可以包括任意数量的像素行和列。通常，像素阵列102例如包括640×480或1024×768微测辐射热计像素。

在图1的示例中，阵列102的每一列像素与相对应的参考结构106相关联。尽管在功能上不是图片元素，但是通过与成像(或有源)微测辐射热计像素104的结构类比，这个结构在本文中将被称为“参考像素”。另外，输出块(OUTPUT)108耦接到像素阵列102的每一列和耦接到参考像素106中的每一个，并提供原始图像RAW，该原始图像RAW包括由像素阵列102结合参考像素106捕获的信号或读数。

例如，在转让给本申请人的美国专利US 7,700,919中更详细地讨论了测辐射热计类型的像素阵列的示例，该专利的内容在法律允许的范围内通过引用结合于此。

控制电路(CTRL)110例如向像素阵列102、向参考像素106和向输出块108提供控制信号。

原始图像RAW例如被提供给图像处理电路(IMAGE PROCESSING)112，该图像处理电路例如对图像的像素应用两点图像校正以及还应用余辉校正，以产生经校正的图像S。

在像素阵列102的读取操作期间，像素行例如一次读出一行。

图2更详细地示出了根据示例实施例的图1的图像处理电路112。

图像处理电路112的功能例如以软件实施，并且图像处理电路112包括处理设备(PROCESSING DEVICE)202，该处理设备具有处于存储在指令存储器(INSTR MEMORY)204中的指令的控制下的一个或多个处理器。在替代性实施例中，图像处理电路112的功能可以至少部分地由专用硬件实施。在这种情况下，处理设备202例如包括ASIC(applicationspecific integrated circuit，专用集成电路)或FPGA(field programmable gatearray，现场可编程门阵列)，并且可以省略指令存储器204。

处理设备202接收原始输入图像RAW，并生成经校正的图像S，该经校正的图像例如被提供给成像设备的显示器(未示出)。

处理设备202还耦接到例如存储用于实施2点图像校正的偏移值(OFFSET)208和增益值(GAIN)210的数据存储器(MEMORY)206。

在一些实施例中，执行无快门图像校正。例如，偏移值呈表示结构列扩散的向量V_COL以及表示由像素阵列102引入的2D非列结构弥散的矩阵OFF_DISP的形式。例如，列扩散主要是由于在每一列中使用参考像素106而导致，而列参考像素的行通常不是完全均匀的。2D非列弥散例如主要是由像素阵列的有源测辐射热计之间的局部物理和/或结构差异引起的，该差异例如是由工艺过程弥散引起的。

在转让给本申请人的于2015年4月24日提交的美国专利申请US 14/695539中更详细地描述了无快门图像校正技术的示例，特别是矢量V_COL和矩阵OFF_DISP的生成，以及基于这个矢量和矩阵的像素值的校正的示例，该申请的内容在法律允许的范围内通过引用结合于此。

在替代性实施例中，可以实施不同类型的2点图像校正，包括涉及使用快门的解决方案。

存储器206还例如存储为一个或多个先前捕获的图像的像素计算的一个或多个余辉估计(REMANENCE ESTIMATE(S))212，如将在下面更详细地描述的那样。在一些实施例中，存储器206还存储表示余辉衰减的模型W(MODEL W)，如将在下面更详细地描述的那样。

图3表示其中存在余辉伪影302的图像300。例如，这个图像是由微测辐射热计阵列捕获的，并且在阵列是有效的同时，由阵列接收眩光以及特别是相对强的热通量。微测辐射热计阵列的每个像素的热层的电阻应该作为温度的函数而变化，温度越高，电阻越低。在图像300中，为了可读性，灰度级已经被反转：高温由较暗的像素表示。相对强热的通量导致相关像素的电阻方面的显著、持续的降低，从而导致图像300中对应于穿过微测辐射热计阵列的通量的路径的暗迹302。当由于太阳的热量导致时，这种伪影通常被称为晒斑。

图4是表示根据示例实施例的作为时间的函数的微测辐射热计阵列(诸如图1的阵列102)的像素x的余辉ω_x的示例的曲线图。余辉ω_x具有可以通过以下等式建模的指数衰减：

[数学式2]

其中：t₀是在眩光击中像素的瞬间之后，像素回到其正常温度的时间；x＝(row,column)是像素的坐标向量；α、β和τ是参数；以及：

[数学式3]

是特征函数。在图4的示例中，τ＝10、α＝1和β＝0.5。

余辉是加性扰动。因此，受余辉影响的图像的像素值f_x(t)可以建模为：

[数学式4]

其中S_x(t)是没有眩光的图像，G_ωx是由于眩光引起的增益变化矩阵，以及f_x(t)例如对应于如上所述的2点图像校正。在法国专利FR2987210中更详细地描述了一种用于基于增益变化矩阵执行校正的方法，该专利的内容在法律允许的范围内通过引用结合于此。然而，由项G_ωx(t-t₀)S_x(t)表示的像素的灵敏度方面的变化可以被认为相对于ω_x很小，并且在下文中，通过将项G_ω视为等于1，其贡献将被忽略。

当IR成像设备(诸如图1的设备100)被用于捕获IR图像时，存在可以被区分的几种不同的用例。

可以考虑动态状况(以下称为情况1)，其中像素阵列的视场中的图像可以被认为是变化的，因为IR相机正在移动和/或因为场景正在改变。对于这种类型的动态状况，以上模型是合适的，因为余辉ω_x(t)和图像S_x(t)两者随时间变化：

[数学式5]

f_x(t)＝ω_x(t)+S_x(t)

另一状况(以下称为情况2)对应于其中像素阵列的视场中的图像基本恒定的情况。因此，只有余辉随时间变化。在这种情况下，受余辉影响的图像的像素值f_x(t)可以建模为：

[数学式6]

f_x(t)＝ω_x(t)+S_x

其中可以假设没有眩光的图像S_x是恒定的。

在情况1和情况2两者中，图像流例如被捕获和存储以形成视频文件和/或近乎实时地被中继到IR相机的图像显示器。

第三种状况(以下称为情况3)对应于这样的情况，在这种情况下，IR相机的视场中的图像和余辉都变化得足够慢，以至于它们两者可以被认为是随时间恒定的，换句话说，时间衰减非常慢，以至于不能在评估的时间范围内对其进行测量。在这种情况下，图像的像素值f_x可以被建模为：

[数学式7]

f_x＝ω_x+S_x

现在将更详细地描述用于至少部分去除由情况1和情况2的余辉引起的伪影的方法。在本公开中将不会涉及情况3，该情况等同于没有关于丢失像素信息的附加时间信息的单个图像的图像修饰。

下面，通过表示t＝n和t+δt＝n+1对时间进行采样。因此，在时间t和像素x处矩阵M的值将被表示为M_n,x，而不是M(t,x)。

总体方案

图5是示出根据本公开的示例实施例的校正图像序列中的余辉的方法500中的操作的一般流程图。该方法例如由图1和图2的图像处理电路112实施。替代性地，该方法可以由另一类型的处理设备来实施，该处理设备可以形成或不形成IR相机的一部分。

图像捕获

在操作501中，该方法在捕获图像f_n时开始。这个图像例如对应于已经对其应用了图像校正(诸如上述的2点图像校正)的原始图像RAW_n。

变化检测

在操作501之后，例如实施操作502，其中基于图像f_n和先前捕获的图像f_n-1检测变化场景的存在。操作502之后是余辉估计方法503，包括例如估计当前图像n的每个像素x处的当前余辉平均值μ_ω,n,x。操作502中的变化场景的检测用于在执行余辉估计操作503的替代性方法504和505之间进行选择。特别地，如果在操作502中检测到变化的场景，则在随后的余辉估计方法503中，对应于图5中的方法504，应用情况1的方法。相反，如果检测到静态场景，则在随后的余辉估计方法503中，应用对应于图5中的方法505的情况2的方法。

例如，在操作502中，可以通过将时间变化与阈值进行比较来检测场景的动态，该阈值被设置为高于时间噪声水平，例如是时间噪声水平的两倍或三倍。可以对图像全局地执行这样的评估，诸如基于平均时间变化，在这种情况下，例如根据相同的方法(情况1的方法504或情况2的方法505)处理图像的所有像素。替代性地，可以基于图像的各个像素或各个子区域处的时间变化来执行局部评估，并且可以逐个像素或逐个子区域地在应用情况1的方法504或情况2的方法505之间进行选择。

方法504在下文标题为情况1的小节中更详细地描述，以及方法505在下文标题为情况2的小节中更详细地描述。

虽然图5示出了包括方法504和505两者的执行的方法，这些方法之间的选择是基于场景变化检测针对每个帧进行的，但是在替代性实施例中，可以仅实施方法504或者仅实施方法505。

余辉校正

在操作503之后，例如实施操作506，包括基于所估计的余辉来校正图像。例如，通过对图像f_n的每个像素x应用以下等式来校正图像，以生成图像S_n：

[数学式8]

S_n,x＝f_n,x-μ_ω,n,x

在操作506之后，n递增，并且操作501至506例如在捕获后续图像n之后重复。

情况1

余辉区域检测

在情况1中，方法503包括应用图5的方法504，包括操作510至513。

在操作510中，更新一个或多个余辉区域。

例如，为每个图像f_n生成指示余辉区域的掩模。这些区域可能大于其中像素被余辉影响的实际区域。然而，在一些实施例中，余辉区域对应于图像的至多75％的像素。例如，可以通过检测其中像素值明显高于周围区中的区域、通过跟踪例如太阳的轨迹来生成掩模，尽管也可以使用其他技术，包括基于边缘检测的技术。

修复

然后，在操作511中，对图像中的那些余辉区域应用修复。修复是本领域已知的图像重构技术，并且包括基于来自目标区域外部的图像信息来校正目标区域内的像素值。例如，在C.Guillemot和O.Le Meur的题为“Image inpainting:Overview and recentadvances(图像修复：概述和最新进展)”，IEEE signal processing magazine,31(1),127-144的2014年出版物(其内容在法律允许的范围内通过引用结合于此)中更详细地描述了修复的方法。

余辉测量

在操作512中，基于经修复的图像，针对每个像素处理图像f_n中的像素x处的余辉的噪声度量ω*_n,x。例如，通过逐像素执行以下减法来计算这个噪声度量：

[数学式9]

其中S*_n,x是像素x处的第n个经修复的图像。

余辉估计

在操作513中，通过应用余辉估计操作530来估计余辉。例如，余辉估计包括对每个像素x的瞬时平均余辉μ_ω,n,x的估计。例如，对于每个像素x，这包括生成余辉的当前估计ω_n,x，作为当前余辉度量ω*_n,x和余辉的先前估计ω_n-1,x之间的加权和。在一些实施例中，操作530包括生成估计ω_n,x简单地作为在前m个图像中所测量的余辉{ω*_k,x}_k＝n-M:n的平均值。在替代性实施例中，操作530包括基于考虑余辉的时间衰减模型的Kalman滤波器生成估计ω_n,x。例如，可以使用指数模型W。使用Kalman滤波器，ω_n,x的估计量例如是它的均值μ_ω,n,x和完全表征了它的其方差σ²ω_,n,x。

方法504允许基于本发明人进行的观察来去除余辉，即当像素阵列的视场中的图像正在改变时，像素值的时间变化可以被认为是具有以μ_ω,n,x为中心的正态分布的噪声。

图6是示出在以例如每秒25帧的图像捕获速率捕获的300幅图像的序列n上受到余辉影响的像素x的模拟像素值(PIXEL VALUE)的示例的曲线图。图6中的点表示300个像素值，以及连续曲线602表示随时间的变化的所估计的余辉平均值。图6中的虚线604表示真实余辉ω_x，其具有指数衰减。通过获得给定图像的这个真实余辉ω_x的合理精确估计，可以将其从像素值中减去，以便提供表示图像场景的基础像素信息。

实际上，考虑到像素阵列的视场中的图像随时间变化，场景可以被认为是如下定义的稳态噪声：

[数学式10]

其中“～”表示“根据……”分布，并且N(μ,σ²)是具有平均值μ和方差σ²的正态分布，平均值和方差在此在时间和空间上是恒定的。

然而，余辉是非稳态噪声，因为它具有随时间衰减的平均值：

[数学式11]

因此，在图像f_n的像素x处观察到的图像是两个独立随机变量的和：

[数学式12]

图7是示出了微测辐射热计阵列上从0到P的像素(PIXEL x)的模拟像素值(PIXELVALUE)的曲线图，其中假设区域K中的像素在晒斑的路径中。

图7以连续曲线702表示观察图像f_n,x的示例，其是由虚线704表示的图像场景μ_S,n的平均值和由连续曲线706表示的平均余辉μ_ω,n,x的总和。因此，通过估计图像场景μ_S,n,x的平均值，可以通过计算以下获得μ_ω,n,x：

[数学式13]

μ_ω,n,x＝f_n,x-μ_S,n,x

在图5的方法504的操作513中，本发明人提出使用修复来获得图像场景μ_S,n的平均值的估计，现在将参考图8更详细地对这进行描述。

图8是示出了微测辐射热计阵列上从0到P的像素(PIXEL x)的像素值(PIXELVALUE)的曲线图，其中类似于图7的示例，假设区域K中的像素在晒斑的路径中。

图8中的连续曲线802表示观察到的图像f_n的线，而虚线804表示被修复的这条线，表示为S*_n，其可以被表达如下：

[数学式14]

其中K_bar是晒斑区域K之外的区域。换句话说，来自晒斑区域之外的区域的像素信息用于修复区域K。

图8中的曲线806代表S的连续样本。它们证明修复曲线804没有偏差：

[数学式15]

E[S^*]＝μ_S

其中E[]是期望值。

虚线808表示余辉ω_n，并且连续曲线810表示相减f_n-S*_n的结果。

可以假设图像的修复部分提供了图像场景的平均值的无偏差估计μ_S：

[数学式16]

以及因此提供μ_ω,n,x的无偏差度量：

[数学式17]

图9是表示余辉测量ω*_n,x与平均余辉μ_ω,n,x之间的差的幅值(MAGNITUDE)的分布(其对应于正态分布)的曲线图。实际上，在每个时间瞬间，对应于以下的ω*_n的度量：

[数学式18]

提供了余辉的期望值E[ω*_n,x]＝μ_ω,n,x的无偏差噪声测量。

由于这个值随时间变化，本发明人提议基于当前测量ω*_n和先前估计μ_ω,n-1之间的加权和来估计每个新图像f_n的这个值，如现在将参考图10更详细描述的那样。

图10是表示沿着n＝1000个连续帧的序列的一个像素x的余辉估计的值z_n的曲线图。每个连续曲线对应于在图5的操作530中用于估计值μ_ω,n,x的不同方法。点对应于原始余辉ω*_n,x测量。线1002表示这些原始余辉度测量ω*_n,x的平均值，其可以表示为：

[数学式19]

图10中的虚线1004表示像素x的余辉的真实衰减。

图10中的连续曲线1006表示基于10个图像的窗口的移动平均值，其可以表示为：

[数学式20]

曲线1006提供了快速收敛，但是噪声相对较大。

图10中的曲线1008表示基于100个图像的窗口的移动平均值，其可以表示为：

[数学式21]

曲线1008仍然提供相对快速的收敛，并且比曲线1006噪声更少。然而，这个曲线是有偏差的，并且这种偏差有在最终经校正的图像中可见的风险。

在一些实施例中，为了在收敛速度和噪声之间获得合理的折衷，可以基于m个图像的滑动窗口来使用移动平均值，其中m在20到150的范围内，并且例如在20到80的范围内。在一些实施例中，图像的数量m可随时间变化，例如随着图像序列进展而随时间增加。在一些实施例中，如果例如结果仍然在某种程度上有偏差，并且如果需要校准滑动窗口的大小，则基于滑动窗口的这种方法可能存在缺点。

图10中的曲线1010对应于Kalman滤波器的应用，与移动平均相比，该Kalman滤波器具有相对较慢的收敛，但是以相对较高的精度收敛到由曲线1004表示的余辉的真实衰减。Kalman滤波器考虑了序列中的多个图像的余辉的测量，以及余辉的指数衰减的模型W。现在将参考图11和图12更详细地描述在操作530中Kalman滤波器的使用。

Kalman滤波

图11是示出在图5的操作513中使用Kalman滤波器估计平均余辉值μ_ω,n,x的方法的操作的流程图。该方法例如由图1和图2的图像处理电路112实施。替代性地，该方法可以由另一类型的处理设备来实施，该处理设备可以形成或不形成IR相机的一部分。

Kalman滤波器非常适合于使用噪声(观察到的)测量ω_n ^*并且在知道一些超参数的情况下在线性高斯动态时间序列的每一步骤n处估计潜在的(隐藏的)状态向量ω_n。这种时间序列是一种Markov链，并且可以用下面的简单序列来完全描述：

[数学式22]

ω_n＝Aω_n-1+v,

ω₀＝μ_ω,0+u,

其中A为转换矩阵，H为测量矩阵，并且参数v、u和η分别是协方差矩阵Γ，Σ_ω,0和R的一些零均值正态加性噪声。

图12示出了这种时间序列的示例，其中隐藏状态向量(在这个示例中具有1维)是圆，并且观察到的度量是点。来自时间(t-δt)的输入是余辉ω(t-δt)的先前校正的估计，以及将时间(t-δt)处的先前估计链接到时间t处的新估计的转换矩阵A。因为在这个时间序列中一切都是高斯的，所以状态向量完全由其均值向量μ_ω和协方差矩阵Σ_ω(它们是待确定的未知值)来表征。

在每一步骤n处，示出的是[Murphy,Machine Learning,a probabilisticperspective,2012p.638]在知道(时间常数)超参数Γ、R、U、A、H和μ₀的情况下，潜在状态ω_n可以使用以下等式估计：

[数学式23]

其中：

[数学式24]

是Kalman增益，以及I是单位矩阵。

换句话说，在每个步骤n，根据测量值ω_n ^*和先验估计(μ_ω,n,x,Σ^- _ω,n,x)来估计新的后验估计(μ_ω,n,x,Σ_ω,n,x)。

应用于余辉滤波

协方差矩阵

在本公开中，隐藏变量ω是每个像素处的余辉，并且是标量。因此，协方差矩阵Γ、R、Σ_ω,0、Σ_ω只是分别表示为γ、r、σ_ω,0 ²、σ_ω ²的方差。

转换矩阵A

在本公开中，隐藏状态是标量，并且转换矩阵A也是标注为a的标量，并且与余辉的模型W直接相关。如上关于图4所述：

[数学式25]

[数学式26]

在本示例中，假设参数β等于1，从而得到如下的ω(t+δt)的模型：

[数学式27]

ω(t+δt)＝aω(t),其中

转换矩阵a例如对于阵列的所有像素是公共的。替代性地，可以为每个像素存储不同的转换矩阵a_x。

测量矩阵H

如上所述，度量ω_n ^*也是实标量值，其使用度量矩阵H与隐藏状态直接相关：

[数学式28]

其中H＝1

初始化

在图11的方法530的操作1101中，通过初始化Kalman滤波器的隐藏状态来初始化Kalman参数。例如，可以针对所有像素x将初始平均值μ_ω,0和方差值σ_ω,0 ²初始化为：

[数学式29]

以及

[数学式30]

先验估计

在操作1102中，基于余辉的噪声度量以及基于Kalman参数和Kalman超参数，为平均余辉生成新的先验估计。例如，在时间n+1，为平均余辉生成新的先验估计：

[数学式31]

其中γ是模型的方差，并且例如是用户定义的。

后验估计

在操作1103中，基于Kalman参数、Kalman超参数和基于经修复的图像确定的当前测量的余辉值ω_n,x ^*，校正余辉的先前新均值和方差值：

[数学式32]

其中：

[数学式33]

是Kalman增益。

换句话说，余辉的后验估计被生成作为先验估计和度量之间的平衡，并且该平衡通过Kalman增益K进行调谐。

序列

再次参考图11，在操作1103之后，例如更新Kalman参数。如图5所示，然后可以对下一图像n+1重复方法530，包括基于经更新的Kalman参数重复操作1102和1103。实际上，参数初始化操作1101例如仅在图5的方法的第一次迭代执行。

情况2

在情况2中，图5的操作503包括方法505，该方法对应于根据本公开的另外的示例实施例的校正图像序列中的余辉的方法。该方法例如由图1和图2的图像处理电路112实施。替代性地，该方法可以由另一类型的处理设备来实施，该处理设备可以形成或不形成IR相机的一部分。

方法505包括操作520和521。

在操作520中，基于先前图像f_n-1为每个像素计算图像f_n中的余辉的度量。例如，通过执行以下减法，为图像的每个像素x计算余辉的度量ω*_n,x：

[数学式34]

ω*_n,x＝f_n,x-f_n-1,x

在操作521中，基于余辉测量和指数衰减模型W来估计每个像素的余辉ω_n,x。例如，余辉估计包括每个像素x的瞬时平均余辉的估计μ_ω,n,x。在一些实施例中，余辉平均值μ_ω,n,x及其方差σ² _ω,n,x是使用Kalman滤波器来估计的，如上文结合图11所述。然而，对于静态场景的情况，测量ω*噪声很多，并且因此，对于情况2，转换矩阵a通常应该比情况1中更精确。因此，在一些实施例中，对于每个像素x，余辉指数衰减的局部转换矩阵a_x＝1-δτ/τ_x被存储在存储器206中。这些转换矩阵例如在预处理操作期间被估计，该预处理操作例如包括回归恰好在余辉出现之后的第一帧中的两个或更多个上的指数衰减以估计τ_x。而且，在情况2中，转换矩阵是：

[数学式35]

本文描述的实施例的优点在于，对于变化场景(情况1)或静止场景(情况2)的情况，可以以相对简单和有效的方式在图像序列中校正余辉。

已经描述了各种实施例和变型。本领域技术人员将理解，可以组合这些实施例的某些特征，并且本领域技术人员将容易想到其他变型。例如，对于本领域技术人员来说显而易见的是，在为所有像素提供共同的指数衰减模型W，或者为每个像素x提供不同的模型之间的选择将取决于所需的精度和IR像素阵列的特性。