图像传感器噪声的数学模型及标定

目录

0. 写在前面

这篇文章应该是博客里到目前为止最枯燥的一篇了…… 写这篇文章的原因有二:一是在我的学位论文中有一部分内容需要涉及图像噪声模型,正好这几天读完 Radiometric CCD Camera Calibration and Noise Estimation,权当记录备忘;二是因为搜了一圈发现目前网上基本没有对相机传感器噪声的数学模型进行详尽介绍的中文资料,所以也算是做一点微小的工作,填补这一块的空白吧 🙄

实际上,对于 CV 领域内大部分应用来说,$f=\mu+\epsilon$ 这种经典的模型已经足以对图像的像素值进行表征了。但是对于一些更偏工程、需要与相机硬件直接打交道的应用,往往还需要对成像过程中各阶段引入的噪声进行更加具体的建模和分析,从而将其最大程度抑制。正因如此,我默认阅读这篇文章的读者都具备了辐射度学、相机成像原理、图像处理以及统计学相关的基础知识,文中不再对一些基本概念进行赘述。

此外,图像噪声这一领域并不是我的研究方向,文中难免存在一些瑕疵或错误,欢迎指正。

1. 相机响应值模型(Camera raw response formation model)

不管是 CMOS 也好 CCD 也好,图像传感器的本质都是将光信号转换为电信号。对于图像传感器上的单个感光单元(collection site/photosensory cell/sensor element,随便什么叫法,总之仅指单个「像素」)来说,其由光电效应释放出的电子个数 $I$ 为
\begin{equation}
I = T\int_\lambda\int_x\int_y\,E(x,y,\lambda)\,S_r(x,y)\,q(\lambda)\,\textrm{d}x\,\textrm{d}y\,\textrm{d}\lambda\,,
\label{eq:electrons_num}
\end{equation}

其中 $x$、$y$ 表示该单元的二维位置坐标,$\lambda$ 是波长,$T$ 为积分时间(假设各参数在快门帘开启的这段时间内均不随时间变化而变化),$E(x,y,\lambda)$ 为该单元内光敏面处的光谱辐照度(incident spectral irradiance),$S_r(x,y)$ 为感光单元的空间响应,$q(\lambda)$ 为光电转换函数。

从量纲角度分析会更好理解一些:积分时间 $T$ 的单位为 $s$,光谱辐照度 $E$ 的单位为 $W/(m^2\cdot{}nm)$,$S_r$ 无量纲,光电转换函数 $q$ 的单位为 $e^-/(J\cdot{}nm)$(在指定波长上每吸收1焦耳能量能释放出的电子个数),因此有 $s\cdot{}(J\cdot{}s^{-1}\cdot{}m^{-2})\cdot{}(e^-\cdot{}J^{-1})\cdot{}nm^{-1}\cdot{}m^2\cdot{}nm = e^-$。

我们注意到,式 \eqref{eq:electrons_num} 中的 $x$、$y$ 均被限制在了单个感光单元内,在这么小的面积内,可以认为各个参数均不随位置变化而变化,因此位置坐标 $x$、$y$ 可以被省略。式 \eqref{eq:electrons_num} 被简化为
\begin{equation}
I = TS_rA\int_\lambda\,E(\lambda)\,q(\lambda)\,\textrm{d}\lambda\,,
\label{eq:electrons_num_simple}
\end{equation}

其中 $A$ 为感光单元内的有效感光面积。

此外,被摄物体表面的光谱辐亮度(spectral radiance)和感光单元内光谱辐照度之间还满足如下关系:
\begin{equation}
E(\lambda) = UL(\lambda)t(\lambda)\,,
\label{eq:radiance2irradiance}
\end{equation}

其中 $L(\lambda)$ 为物方平面上被摄表面的光谱辐亮度,$t(\lambda)$ 为整个系统的光谱透过率函数(包括镜头、低通、Bayer filters 等等),$U$ 是一个由 vignetting、shading 等因素导致的空间非均匀性调制函数。

因此,联立式 \eqref{eq:electrons_num_simple} 与 \eqref{eq:radiance2irradiance} 有:
\begin{equation}
I = TS_rAU\int_\lambda\,L(\lambda)\,t(\lambda)\,q(\lambda)\,\textrm{d}\lambda\,,
\label{eq:formation}
\end{equation}

式 \eqref{eq:formation} 就是数字成像系统将光信号变为电信号的理想物理模型。

式 \eqref{eq:radiance2irradiance} 里的 $U$ 显然是一个偷懒的简化。

实际上,很多文献里都会提到,$L(\lambda)$ 与 $E(\lambda)$ 之间满足那个著名的 cos 四次方定理:
\begin{equation}
E(\lambda)\approx\frac{\pi}{4F^2}cos^4(\alpha)L(\lambda)t(\lambda)\,,
\label{eq:cos4}
\end{equation}

其中 $F$ 是光圈数(焦距/光阑直径),$\alpha$ 是某个像素对应的主光线相对于光轴的偏离角。

然而,对于一个真实存在的光学系统来说,照度的非均匀性绝对不仅仅是 cos 四次方衰减那么简单,还需要考虑各种复杂的原因,比如 vignetting、镜头设计层面的照度补偿等等等等。所以为了不失一般性(也为了表述方便一点),我直接简单粗暴地把这些非均匀性统一放到这个 $U$ 里去了。

之后,感光单元内产生的电信号将经由放大电路(调 ISO 其实就是在调节放大电路的增益系数)进行放大,并经 ADC 转换为数字信号缓存在 RAM 中,待后续存储或调用:
\begin{equation}
D = round\left(\frac{g^\prime{}I + V_{offset}}{\eta}\right)\,,
\label{eq:amplify}
\end{equation}

其中 $g^\prime$ 是模拟电路的增益倍数,$V_{offset}$ 是人为加入的一个偏置电压,通常可以通过这种方式来避免小于零的输出信号被 ADC 截零从而改变噪声分布的对称性(一些相机的暗电流就是这么来的),$\eta$ 是 ADC 中的量化步长(quantization step),单位为 $volts/DN$(DN 就是数字信号的数值,digital number,理论上是无量纲的,但是为了表述方便一般会人为加上一个量纲),$round(\cdot)$ 表示 ADC 的取整操作。

最后这个供输出的数字信号 $D$,我习惯称之为相机的原始响应值(camera raw response),有一些资料里也会把它称作 DN (digital number) 或 ADU (analog-to-digital units)。

感光单元释放出的电荷如何进入 ADC 这一细节有待确定。老一些的相机,尤其是 CCD 相机,似乎会为传感器阵列的每一行配备一个寄存器,每个寄存器接收一整行感光单元传来的电荷后再统一传给放大器,各行间并行处理;但是对于新一些的 CMOS,甚至会为每一个像素单独配置一个寄存器和放大器,所有信号的处理以像素为单位并行,从而大幅缩短两次采集之间的间隔。相机硬件的设计我并不在行,以上内容仅仅是道听途说,不过并不影响接下来将要进行的噪声的分析。

此外,现在的一些相机还会为传感器配置两组具有不同电容的感光单元电荷输出电路,从而在电信号进入放大电路之前就有两组电压值可供切换,即所谓的 Dual Native ISO

最后,联立式 \eqref{eq:formation} 和 \eqref{eq:amplify},就得到了最终的光信号到数字信号的完整模型
\begin{equation}
\boxed{\begin{split}
{}\\[-.2em] \qquad{}D = round\left(gTS_rAU\int_\lambda\,L(\lambda)\,t(\lambda)\,q(\lambda)\,\textrm{d}\lambda + D_{offset}\right)\,,\qquad\\[-.2em] {}
\end{split}}
\label{eq:response_model}
\end{equation}

显然,这里的 $g=g^\prime/\eta$,$D_{offset}=V_{offset}/\eta$。

2. 图像噪声来源(Noise sources)

理论上,根据式 \eqref{eq:electrons_num_simple}、\eqref{eq:amplify},在完全相同的条件下反复按下快门,对于 raw 图像上的同一个像素,在不同帧上应该具有完全相同的原始响应值。但是实际情况中,在以上光信号 $\rightarrow$ 电信号 $\rightarrow$ 数字信号的过程中必然伴随着噪声的引入。

在对图像传感器的噪声进行建模时需要再次考虑各感光单元的位置信息,因此,将式 \eqref{eq:response_model} 扩展到整个传感器平面上:
\begin{equation}
D(i,j) = round\left[Tg(i,j)S_r(i,j)A(i,j)U(i,j)\int_\lambda\,L(\lambda)\,t(i,j,\lambda)\,q(i,j,\lambda)\,\textrm{d}\lambda + D_{offset}(i,j)\right]\,.
\label{eq:response_model_spatial}
\end{equation}

具体来说,光信号到数字信号的过程中将伴随着以下几类噪声:

i. 像素响应非均匀性(Pixel response non-uniformity)

由于制造加工过程中不确定因素的存在,传感器各个感光单元之间在有效感光面积、光电转换效率、非均匀性、滤色片光谱透过率等属性上必然存在一定的差异。这种像素之间差异通常被称为像素响应非均匀性(pixel response non-uniformity, PRNU)。稍微做一些简化,如果不考虑 $q(\lambda)$ 在各个波长上的偏差,仅认为其存在幅值上的误差的话,对于某一个感光单元,可以使用非均匀性分布函数
\begin{equation}
K(i,j) = \frac{g(i,j)}{g_0}\frac{S_r(i,j)}{S_{r,0}}\frac{A(i,j)}{A_0}\frac{t(i,j)}{t_0}\frac{q(i,j)}{q_0}
\label{eq:prnu}
\end{equation}

来表征各分项误差组合后的整体误差,其中 $g_0$ 表示一个无任何制造误差的感光单元所对
应的 $g$ 值(综合增益),$S_{r,0}$、$A_0$、$t_0$、$q_0$ 同理。基于现有的半导体制造工艺, $K$ 的数值仅在单位强度附近很小的一段区间内波动,因此有
\begin{equation}
\bar{E}(K) = 1\,,\qquad\bar{var}(K)\ll 1\,.
\label{eq:prnu_stat}
\end{equation}

这里的 $\bar{E}(\cdot)$ 和 $\bar{var}(\cdot)$ 分别表示在空间维度上计算均值(mean)与样本方差(sample variance),下同。

ii. 热噪声(Thermal noise)

光电转换器件的热效应将产生一定数量的服从泊松分布的自由电子,这些自由电子与光电效应所激发的光生电子组合后共同从感光单元中释放出来,从而使信号出现一定程度的波动,因此此类噪声也被称为热噪声(thermal noise)$N_{Th}$ 。热噪声的强度仅仅取决于曝光时间 $T$ 以及环境温度,而与光电效应产生的电子数 $I$ 无关。

热噪声是一随机变量,根据泊松分布的性质有
\begin{equation}
E(N_{Th}) = var(N_{Th}) = \mu_{Th}\,,
\label{eq:thermal_noise_stat}
\end{equation}

其中 $\mu_{Th}$ 是一个与曝光时间及环境温度有关的变量,$E(\cdot)$、$var(\cdot)$ 分别表示计算数学期望(expectation)与方差(variance),下同。

iii. 散粒噪声(Shot noise)

根据量子理论,从感光单元中释放出的电子在数量上存在一个随机的涨落,该涨落即为散粒噪声(photon shot noise)$N_S$。散粒噪声是由大量单个事件的统计不确定性引起的,其中主要包括输入光子散粒噪声、光生电流散粒噪声与热效应散粒噪声。散粒噪声的存在使得感光单元最终释放出的有效电子数服从泊松分布:$(I^\ast+N_S)\sim{}Poisson(I^\ast)$,其中 $I^\ast=KUI_0+N_{Th}$ 表示在不考虑散粒噪声的情况下该感光单元应释放出的电子个数,且
\begin{equation}
I_0 = TS_{r,0}A_0\int_\lambda{}L(\lambda)t_0(\lambda)q_0(\lambda)\mathrm{d}\lambda
\label{eq:I0}
\end{equation}

表示一个“完美”的不含任何噪声的感光单元理论上释放出的电子数

根据泊松分布的性质,有 $E(I^\ast+N_S)=var(I^\ast+N_S)=I^\ast$,因此,对于散粒噪声 $N_S$ 自身而言,有
\begin{equation}
E(N_S) = 0\,,\qquad{}var(N_S) = I^\ast = KUI_0+N_{Th}\,.
\label{eq:shot_noise_stat}
\end{equation}

iv. 固定模式噪声(Fixed pattern noise)

对于 CMOS 来说,其各个感光单元所对应的模拟电路中有可能存在一定的噪声,此类噪声通常被称为固定模式噪声(fixed pattern noise, FPN)$N_{FP}$。固定模式噪声与像素响应非均匀性类似,均属于系统误差,两者的区别在于,固定模式噪声是由电路的制造误差引起的,与光电效应无关,且通常呈现出一定的空间分布规律。

v. 读出噪声(Read-out noise)

信号放大单元把模拟电压信号进行放大的过程中会引入一定的读出噪声(read-out noise)$N_R$。读出噪声的期望为零,其方差与放大单元的增益系数 $g$ 成线性正相关。

vi. 量化误差(Quantization error )

由于模数转换单元输出的数字信号需要以整数的形式存储与寄存器中,因此这一连续信号转化为离散信号的过程中必然会引入一定的量化误差(quantization error)$N_Q$。一些文献中量化误差也被称为取整误差。若进入模数转换单元的模拟电压信号等概率地分布于实数轴上,则量化误差将服从区间 $[-\tfrac{1}{2}, \tfrac{1}{2}]$ 内的均匀分布:$N_Q\sim{}\mathcal{U}(-\tfrac{1}{2}, \tfrac{1}{2})$,因此有
\begin{equation}
E(N_Q) = 0\,,\qquad{}var(N_Q) = \tfrac{1}{12}\,.
\label{eq:quantization_error_stat}
\end{equation}

综合考虑以上六种类型的噪声的话,式 \eqref{eq:response_model_spatial} 变为
\begin{equation}
\boxed{\begin{split}
{}\\[-.2em] \quad{}D(i,j) & = g_0\Big[K(i,j)U(i,j)I_0 + N_{Th}(i,j) + N_S(i,j) + N_{FP}(i,j) + N_R(i,j)\Big]\quad\\[.6em] & + D_{offset}(i,j) + N_Q(i,j)\,,\\[-.2em] {}
\end{split}}
\label{eq:noise_model}
\end{equation}

其中 $I_0$ 的含义同式 \eqref{eq:I0},表示一个完美的感光单元理论上释放出的电子数。

式 \eqref{eq:noise_model} 即为带有噪声项的数字相机原始响应值构成模型

显然,由于热噪声、散粒噪声、读出噪声、量化误差均为随机变量,因此式 \eqref{eq:noise_model} 中的原始响应值 $D$ 也为一随机变量。


下文中的内容有误,请勿阅读!我将尽快订正。

3. 噪声特性(Noise characteristics)

把上节列出的不同种类噪声整理一下。综合考虑上节中加粗的五种噪声源,对于最终输出的原始响应值信号 $D$,有(电子数理论值 $I_0$ 的下标被省略,下同)
\begin{equation}
D=(KI + N_{DC} + N_S + N_R)A + N_Q\,,
\label{eq:D}
\end{equation}

其中

  • $I$:感光单元释放的电子个数理论值,见式 \eqref{eq:electrons_num_simple};
  • $A$:放大电路 + 其他电子元件的组合增益;
  • $K$:固定模式噪声(fixed pattern noise),$E[K] = 1$,$\sigma_K^2\ll1$,分布未知,可预先标定($E[\cdot]$ 表示在 $(i,j)$ 平面计算数学期望,下同);
  • $N_{DC}$:暗电流噪声(dark current noise),$E[N_{DC}]=I_{DC}=I_{offset}+E[N_{thermal}]$,$\sigma_{DC}^2$ 未知,正态分布,在一定温度下可预先标定;
  • $N_S$:散粒噪声(shot noise),$E[N_S]=0$,$\sigma_S^2 = KI + N_{DC}$,泊松分布(平移至 $E[N_S]=0$),无法预先标定;
  • $N_R$:读出噪声(read noise),$E[N_R]=0$,$\sigma_R^2$ 未知,分布未知,无法预先标定;
  • $N_Q$:量化误差(quantization error),$E[N_Q]=0$,$\sigma_Q^2=\frac{1}{12}$,均匀分布,无法预先标定。

仅考虑无法预先标定的三种噪声源:散粒噪声、读出噪声、量化误差。三者之间不存在相关性,因此,体现在原始响应值上的最终噪声 $\sigma_N^2$ 可以表示为
\begin{equation}
\begin{split}
\sigma_N^2 & = (\sigma_S^2+ \sigma_R^2)A^2 + \sigma_Q^2\\[.5em] & = \big[KI + N_{DC} + \sigma_R^2\big]A^2 + \sigma_Q^2\,.
\end{split}
\label{eq:noise}
\end{equation}

顺带推导一下信噪比计算公式。信噪比(signal-noise-ratio)定义为信号强度的平方 $\mu^2$ 与方差 $\sigma^2$ 之比,因此有
\begin{equation}
\begin{split}
S/N & = 10\log_{10}\left(\frac{E[D^2]}{E\left[\sigma_N^2\right]}\right) = 20\log_{10}\left(\frac{E[D]}{E\left[\sigma_N\right]}\right)\\[.8em] & = 20\log_{10}\left(\frac{E\big[(KI + N_{DC} + N_S + N_R)A + N_Q\big]}{\sqrt{\big(E[KI] + E[N_{DC}] + \sigma_R^2\big)A^2 + \sigma_Q^2}}\right)\,.
\end{split}
\label{eq:snr}
\end{equation}

计算信噪比时暂时忽略小量 $\sigma_Q^2$,由此可以得到
\begin{equation}
S/N = 20\log_{10}\left(\frac{I + I_{DC}}{\sqrt{I + I_{DC} + \sigma_R^2}}\right)\,.
\label{eq:snr_simple}
\end{equation}

因此可以得出那个众所周知的结论:要提高图像的信噪比,只能增加传感器上收集的电子数。调节放大增益(ISO)或者一味地增加位深($b$)都是没有任何帮助的(cf. Noise, Dynamic Range, and Bit Depth)。而要增加电子数,要么延长积分时间 $T$,要么改进光电转换函数 $q$ 使其具有更高的转换效率(Dual Native ISO 技术可以等效视为为每个感光单元都配备了两组 $q$)。

另外,从式 \eqref{eq:noise} 也可以看出,在其他条件不变的前提下,拉高 ISO(增大 $A$)只会导致噪声的增加,因此对于处理 raw 图像的同学来说,前期拉高 ISO 和后期调节图像亮度是没有本质区别的。

4. 噪声标定(Noise calibration)

如果只是要了解图像噪声的来源和大致性质的话,看完上面一节就够了。这一节简单谈一下如何对不同种类的噪声进行标定。因为标定过程需要引入传感器平面的位置坐标,这一节里的各种公式看起来可能不会那么友好 😳

首先把式 \eqref{eq:D} 扩展到整个传感器二维平面,其位置坐标用 $i$、$j$ 进行表示,则有
\begin{equation}
D(i,j) = \big[K(i,j)I(i,j) + N_{DC}(i,j) + N_S(i,j) + N_R(i,j)\big]A + N_Q(i,j)\,.
\label{eq:D2d}
\end{equation}

下文中所有与空间位置有关的参数都将显式地用 $(i,j)$ 进行表示。若不带位置记号则表示该参数为一位置无关的常量。

我们将 $D(i,j)$ 拆解为两部分:期望不为零的部分 $\mu(i,j)$ 以及期望为零的部分 $N(i,j)$:
\begin{equation}
D(i,j) = \mu(i,j) + N(i,j)\,.
\label{eq:decomposition}
\end{equation}

其中 $\mu(i,j)$ 决定了响应值的强度,$N(i,j)$ 决定了响应值波动的剧烈程度。

根据前一节中各种噪声特性的分析结果,有
\begin{equation}
\mu(i,j) = K(i,j)I(i,j)A + N_{DC}(i,j)A\,,
\label{eq:mu}
\end{equation}

以及
\begin{equation}
N(i,j) = N_S(i,j)A + N_R(i,j)A + N_Q(i,j)\,.
\label{eq:N}
\end{equation}

4.1 噪声的分解(Noise decomposition)

为了方便后续标定过程中的计算,$N(i,j)$ 又可以进一步地被拆解为两部分:第一部分 $N_I(i,j)$ 与感光单元收集的电子数有关,第二部分 $N_C(i,j)$ 与其无关
\begin{equation}
\left\{\begin{split}
N_I(i,j) & = N_S(i,j)A\,,\\[.8em] N_C(i,j) & = N_R(i,j)A + N_Q(i,j)\,.
\end{split}\right.
\label{eq:N_decomposition}
\end{equation}

现在对噪声一项一项进行分析。首先看 $N_I(i,j)$。由于 $N_S(i,j)$ 满足泊松分布,因此 $N_I(i,j)$ 的方差 $\sigma_I^2(i,j)$ 满足
\begin{equation}
\begin{split}
\sigma_I^2(i,j) & = A^2\sigma_S^2(i,j)\\[.8em] & = A^2\big[K(i,j)I(i,j) + N_{DC}(i,j)\big]\,.
\end{split}
\label{eq:NI_var}
\end{equation}

再看 $N_C(i,j)$。$N_C(i,j)$ 的方差 $\sigma_C^2$ 也很容易计算:
\begin{equation}
\begin{split}
\sigma_C^2 & = A^2\sigma_R^2(i,j) + \frac{1}{12}\\[.2em] &= A^2\sigma_R^2 + \frac{1}{12}\,,
\end{split}
\label{eq:NC_var}
\end{equation}

其中,读出噪声 $N_R$ 是一个全局噪声,并不依赖于某个感光单元所在的具体位置,因此其方差 $\sigma_R^2$ 的表达式中可以省去位置坐标。

至此,可以得到式 \eqref{eq:noise} 的改进版本,即,带有位置坐标的最终噪声表达式为
\begin{equation}
\begin{split}
\sigma_N^2(i,j) & = \sigma_I^2(i,j) + \sigma_C^2\\[.8em] & = A^2\big[K(i,j)I(i,j) + N_{DC}(i,j)\big] + \sigma_C^2
\end{split}
\label{eq:noise_position}
\end{equation}

4.2 传感器整体噪声估计(Estimation of sensor noise)

💡 在第 4.5 小节中会提到,由于光学系统的限制,即使对于完全均匀的场景,传感器不同区域也会有不同的照度 $E(i,j,\lambda)$,从而对应不同的释放电子数 $I(i,j)$。对于位于传感器中心与边缘的感光单元,两者有可能相差数倍以上。因此从这一小节开始,对于实际噪声的标定,若无特别说明的话,均指对于传感器某一特定区域进行计算。若要对整个传感器进行标定,可采用分块计算的方法($5\times5$ 一般就够了),最后再对各块的参数进行插值从而得到全局的噪声参数估计。

对一台相机在原始响应值上的最终噪声进行估计通常有两种方法,一是采集两幅 raw 图像 $D_1(i,j)$ 及 $D_2(i,j)$,并在差值图像 $D_\Delta(i,j) = D_1(i,j) - D_2(i,j)$ 的基础上进行噪声方差估计;二是采集 $S$ 幅图像,然后在均值图像 $D_{avg}(i,j) = \frac{1}{S}\sum\limits_{k=1}^SD_k(i,j)$ 的基础上进行估计。第一种方法所需工作量略小一些,并且当目标区域含有的像素数足够多时,其一致性估计的结果(consistent estimates)可以认为非常接近真值,因此这里仅对第一种方法进行介绍。

对于两幅在完全一致条件下(相同的拍摄环境、相同的光源、相同的被摄物体)捕获的 raw 图像,根据式 \eqref{eq:decomposition},有
\begin{equation}
\left\{\begin{split}
D_1(i,j) & = \mu(i,j) + N_1(i,j)\,,\\[.6em] D_2(i,j) & = \mu(i,j) + N_2(i,j)\,.
\end{split}\right.
\label{eq:D1D2}
\end{equation}

对这两幅 raw 图像逐像素计算差值,则有
\begin{equation}
\begin{split}
D_\Delta(i,j) & = D_1(i,j) - D_2(i,j)\\[.6em] & = N_1(i,j) - N_2(i,j)\,.
\end{split}
\label{eq:delta}
\end{equation}

显然,$E\big[D_\Delta(i,j)\big]=0$,且由于 $N_1(i,j)$、$N_2(i,j)$ 不具备相关性,$D_\Delta(i,j)$ 的方差为 $2\sigma_N^2(i,j)$。

首先对图像中的响应值强度进行估计。因为 $N(i,j)$ 的期望为零,因此可以利用 $D_1(i,j)$ 与 $D_2(i,j)$ 计算得到传感器某一区域内 $\mu$ 的一致性估计:
\begin{equation}
\hat{\mu} = \frac{1}{2mn}\left[\,\sum_{j=1}^n\sum_{i=1}^m\,D_1(i,j) + \sum_{j=1}^n\sum_{i=1}^m\,D_2(i,j)\,\right]\,,
\label{eq:mu_estimate}
\end{equation}

其中 $m$、$n$ 分别为传感器在该区域内的横、纵方向上的感光单元个数。

同理,也可以对传感器噪声进行估计:
\begin{equation}
\hat{\sigma}_N^{\,2} = \frac{1}{2(mn-1)}\,\sum_{j=1}^n\sum_{i=1}^m\,\left[D_\Delta(i,j)-\overline{D}_\Delta\right]^2\,,
\label{eq:N_estimate}
\end{equation}

其中 $\overline{D}_\Delta$ 是 $D_\Delta(i,j)$ 在该区域内的空间平均值。

当光信号较强或积分时间较长时,感光单元释放出足够的电子数,此时散粒噪声 $N_S$ 的泊松分布可用正态分布进行近似(cf. Normal approximation to Poisson distribution)。再次忽略量化误差 $N_Q$,则总体噪声仍然为正态分布。此时,$\sigma_N^2$ 的估计值 $\hat{\sigma}_N^2$ 是一个关于真值的高阶小量:
\begin{equation}
var({\hat{\sigma}_N^{\,2}}) = \frac{2(\sigma_N^2)^2}{mn-1} \approx \frac{2(\hat{\sigma}_N^{\,2})^2}{mn-1}\,,
\label{eq:N_estimate_var}
\end{equation}

因此这种估计方法是相当可靠的。

4.3 增益系数 $A$ 及电子数无关噪声 $\sigma_C^2$ 的估计(Estimation of amplification factor and electrons-independent noise)

这部分的内容会稍微 tricky 一些 😛

对式 \eqref{eq:mu} 与式 \eqref{eq:noise_position} 同时计算空间平均值,有
\begin{equation}
\left\{\begin{split}
\overline{\mu} & = \overline{I}A + I_{DC}A\,,\\[.8em] \overline{\sigma}_N^{\,2} & = A^2(\overline{I} + I_{DC}) + \sigma_C^2\,.
\end{split}\right.
\label{eq:mu_sigma_means}
\end{equation}

显然,$\overline{\mu}$ 与 $\overline{\sigma}_N^{\,2}$ 满足线性关系:
\begin{equation}
\overline{\sigma}_N^{\,2} = A\overline{\mu} + \sigma_C^2\,.
\label{eq:linear}
\end{equation}

因此,只要拍摄 $P$ 组具有不同信号强度的 $D_1(i,j)$、$D_2(i,j)$(每一组对应不同的 $\overline{I}$,可通过改变光源亮度或者在相机前放置中性滤光片来实现),并使用上一小结中的方法对每一组的信号强度和噪声进行估计,用 $\hat{\mu}$、$\hat{\sigma}_N^{\,2}$ 替代 $\overline{\mu}$、$\overline{\sigma}_N^{\,2}$,即可在 $(\mu,\sigma_N^2)$ 平面内对 $P$ 个样本点回归出一条斜率为 $A$、截距为 $\sigma_C^2$ 的直线。当不同组对应的 $\hat{\sigma}_{N_q}^{\,2}$ 来自同一个满足正态分布的总体时,采用最大似然估计,有
\begin{equation}
\hat{A},\,\hat{\sigma}_C^{\,2} = \operatorname*{argmin}_{A,\,\sigma_C^2}\sum\limits_{q=1}^P\frac{\left[\hat{\sigma}_{N_q}^{\,2}-\left(A\hat{\mu}_q+\sigma_C^2\right)\right]^2}{var(\hat{\sigma}_{N_q}^{\,2})}\,
\label{eq:mle}
\end{equation}

其中 $(\cdot)_q$ 表示该数据来自第 $q$ 组图像($0\lt{}q\le{}P$)。最大似然估计的详细计算过程见 Radiometric CCD Camera Calibration and Noise Estimation 附录。

4.4 暗电流噪声 $N_{DC}$ 的估计(Estimation of dark current noise)

💡 注意:这一小节里关于暗电流噪声的估计是指整个传感器平面上的全局估计,不必再考虑照度 $E(i,j,\lambda)$ 的影响——因为在这里的实验设置下根本就没有照度了嘛。

暗电流噪声的标定应该是各种噪声中最简单的一部分了。

假设没有任何光信号进入相机系统(例如在暗室中拍摄,或者直接盖上镜头盖和 EVF),则根据式 \eqref{eq:D2d},有
\begin{equation}
D_{dark}(i,j) = \big[N_{DC}(i,j) + N_S(i,j) + N_R(i,j)\big]A + N_Q(i,j)\,.
\label{eq:D2d_dark}
\end{equation}

易知,此时 $D_{dark}(i,j)$ 在时域上的期望为 $N_{DC}(i,j)A$,方差与式 \eqref{eq:noise_position} 相同,仍为 $\sigma_N^2(i,j)$。

由于每次按下快门捕获的 $D_{dark}(i,j)$ 之间并不具备相关性,因此,只要拍摄 $S$ 幅这样的全黑图像并进行平均,就可以得到关于 $N_{DC}(i,j)$ 的一致性估计:
\begin{equation}
\hat{N}_{\!DC}(i,j) = \frac{\sum\limits_{k=1}^SD_{dark,k}(i,j)}{S\hat{A}}\,,
\label{eq:I_DC_estimate}
\end{equation}

且 $var\big[\hat{N}_{\!DC}(i,j)\big] \approx \frac{\hat{\sigma}_N^{\,2}}{S\hat{A}^2}$。

在标定暗电流噪声时偏置电流 $I_{offset}$ 的重要性就体现出来了。如果没有偏置电流,则有可能对于某一个像素,$N_S + N_{thermal} + N_R$(散粒噪声 + 热噪声 + 读出噪声)的理论值小于零,但是经过 ADC 后该信号必然被置为零。此时 $D_{dark}(i,j)$ 在时域上的期望必然会略大于 $N_{DC}(i,j)A$,如果再使用式 \eqref{eq:I_DC_estimate} 进行估计,得到的 $\hat{N}_{\!DC}(i,j)$ 自然也就不再准确了。

4.5 固定模式噪声 $K(i,j)$ 的估计(Estimation of fixed pattern noise)

结合式 \eqref{eq:D2d} 和 \eqref{eq:N},有
\begin{equation}
D(i,j) = \big[K(i,j)I(i,j) + N_{DC}(i,j)\big]A + N(i,j)\,.
\label{eq:D2d_decomposition}
\end{equation}

因为 $N(i,j)$ 的期望为零,且 $N_{DC}(i,j)$ 的估计可以由式 \eqref{eq:I_DC_estimate} 得到,因此,只要拍摄足够多的 raw 图像,就能够得到关于 $K(i,j)I(i,j)A$ 的估计:
\begin{equation}
\begin{split}
e(i,j) & = K(i,j)I(i,j)A \\[.6em] & \approx D_{avg} - \hat{N}_{\!DC}(i,j)\hat{A}\\[.4em] & = \frac{1}{S}\sum\limits_{k=1}^SD_k(i,j) - \hat{N}_{\!DC}(i,j)\hat{A}\,,
\end{split}
\label{eq:KI_estimate}
\end{equation}

在理想的情况下,如果能保证传感器各个感光单元都具有严格相等的的 $I(i,j)$,则可以轻松地估计出 $K(i,j)$:
\begin{equation}
\hat{K}(i,j) = norm\big[e(i,j)\big]\,,
\label{eq:K_estimate_simple}
\end{equation}

其中 $norm[\cdot]$ 表示空间维度上的归一化(因为 $E[K(i,j)] = 1$)。

但是!对于一个实际的成像系统,哪怕视场中的光源和被摄物体表面绝对均匀(比如用这种积分球产生一个面光源),其传感器上不同位置的感光单元硅晶元面上的照度,也会受到一个与位置相关的因子 $\kappa(\rho,\lambda)$ 的调制,其中 $\rho$ 表示该感光单元距离传感器中心的距离(假设镜组和传感器的装配过程中不存在误差),$\lambda$ 表示波长,且有 $\kappa(\rho,\lambda)\le1$。以最简单的单透镜镜组为例,对于轴外的感光单元,其照度与夹角 $\alpha$ 的余弦值的四次方成正比,其中 $\alpha$ 表示该感光单元与被摄点的连线与光轴的夹角(cf. 式 \eqref{eq:rad2irra} 以及 Derivation of the "Cosine Fourth" law for falloff of illuminance across a camera image)。显然,距离 $\rho$ 越大,夹角 $\alpha$ 也越大。同时,由于实际镜组不可能做到无穷大,必然还有渐晕效应的存在,进一步加剧了轴外像素的照度衰减。当然,现在的镜头经过精心的设计,在一定程度上缓解了这种衰减,但是在实际相机捕获的图像中,这种「中心亮四周暗」的情况仍然或多或少存在着,比如下图。


传感器平面的照度衰减示意图。工业界习惯把这种 falloff 称为 lens shading。由于色差(chromatic aberration 那个色差,不是 color difference 哦)的存在,轴外感光单元的照度衰减 $\kappa(\rho,\lambda)$ 是波长相关的函数,从而导致 raw 图像的响应值不单是幅值下降,还有可能出现偏色,俗称 color lens shading。这种 color lens shading 受到光源光谱功率分布、物体光谱反射比的共同影响,理论上是无法进行标定的(因为你永远无法预测相机将要拍摄的内容)。

这种光学系统层面产生的非均匀性往往需要在 ISP 中单独进行校正,因此不在本文中讨论范畴之内。我们只要知道,对于实际的成像系统,基本不可能保证传感器各处的 $I(i,j)$ 是均匀的。

既然 $I(i,j)$ 无法被视作一常数,我们必须寻找其他方法对 $I(i,j)$ 进行估计。

类似第 4.3 小节里的那样,还是采集具有不同信号强度的 $P$ 组图像,每组图像使用完全相同的条件拍摄 $S$ 张并求平均,以消除 $N(i,j)$ 的影响。此时,对于第 $q$ 组强度对应的图像($1\le{}q\le{}P$),根据式 \eqref{eq:KI_estimate},可以得到关于 $K(i,j)I_q(i,j)A$ 的估计 $e_q(i,j)$:
\begin{equation}
\begin{split}
e_q(i,j) & \approx D_{avg,q} -\hat{N}_{\!DC}(i,j)\hat{A}\\[.6em] & = \frac{1}{S}\sum\limits_{k=1}^SD_{k,q}(i,j) - \hat{N}_{\!DC}(i,j)\hat{A}\,,
\end{split}
\label{eq:KI_estimate_q_configuration}
\end{equation}

其中 $D_{avg,q}$ 表示第 $q$ 组取平均后的图像,$D_{k,q}(i,j)$ 表示在第 $q$ 组信号强度下拍摄的第 $k$ 幅图像。

对于图像中的每一个像素 $(i,j)$,如果我们可以找到一个以其为中心且大小合适的区域,使得该区域内 $E[K(i,j)] \approx 1$,那么就可以得到 $e(i,j)$ 在该区域内的平均值 $\tilde{e}(i,j) \approx I(i,j)A$,进而利用
\begin{equation}
K(i,j) = \frac{K(i,j)I(i,j)A}{I(i,j)A} \approx \frac{e(i,j)}{\tilde{e}(i,j)}
\label{eq:K_estimate}
\end{equation}

估计出 $K(i,j)$。需要注意的是,这个区域既不能选得太小,否则无法保证 $E[K(i,j)] \approx 1$;也不能选得太大,否则 $I(i,j)$ 将受到 $\kappa(\rho,\lambda)$ 的影响出现波动,从而导致 $e(i,j)$ 在该区域内的平均值不再等于 $I(i,j)A$。对于实际的相机,这个区域取 $9\times9$ 至 $15\times15$ 都没有问题。

根据这种方法,对于 $D_{avg,q}$ 中的每一个像素,都可以得到一组 $e_q(i,j)$ 和 $\tilde{e}_q(i,j)$ 数据,对应了 $(e_q, \tilde{e}_q)$ 平面内一个样本点。那么对于 $P$ 个样本点,必然能够拟合出一条经过原点的直线,该直线的斜率便对应了该像素的 $K(i,j)$。利用最小二乘法,有
\begin{equation}
\hat{K}(i,j) = \frac{1}{P}\sum_{q=1}^P\frac{e_q(i,j)}{\tilde{e}_q(i,j)}\,.
\label{eq:K_estimate_lsq}
\end{equation}

在实际标定实验的操作过程中,被摄物体表面有可能存在脏点,因此建议在 $P$ 组图像的采集过程中,每拍完一组都重新移动一次物体,这样即使在 $P$ 个样本中存在一两个异常点,也可以很方便地进行 outlier removal,不会影响到最终 $K(i,j)$ 的拟合精度。

References

11 Comments

  1. Anonymous 2019-05-05 Reply

    关于科技论文公式的写法有点疑问,例如博主的公式(1),dxdy中的x,y貌似用了斜体?\lambda用了正体?
    正好在写论文,想问问看正规的写法应该是哪个?

  2. Anonymous 2019-03-04 Reply

    文章写的真好,讲明了很多道理,祝你在以后有更大的发展

  3. laixintao 2019-02-28 Reply

    我现在发现 Chrome 也好了。首页的排版也没有问题。有可能是网络环境的问题吧…现在不好重现了。 :)

  4. nick 2019-02-28 Reply

    现在是在虹软上班了吗?

    • Author
      Captain 2019-02-28 Reply

      没有,后来决定不去了

  5. nick 2019-02-28 Reply

    看了乐高入坑指南,发现博主很久没有更新乐高了,是退坑了吗

    • Author
      Captain 2019-02-28 Reply

      退坑是不可能了,只是抽不出时间玩了

  6. laixintao 2019-02-26 Reply

    Chrome 有问题,显示不出来,Safari 是可以的

    • Author
      Captain 2019-02-27 Reply

      是不是 Chrome 把 JavaScript 给禁止了?另外我今天发现用 2K 显示器时,博客首页在 Chrome 下排版会出错,你那里有这种情况么?

  7. laixintao 2019-02-25 Reply

    博主,你的公式渲染是否有问题,还是我浏览器的问题?

    • Author
      Captain 2019-02-25 Reply

      我这里看着没有问题啊,你用什么浏览器?

Leave a reply to Captain Click here to cancel the reply

Your email address will not be published.