这篇文章里有超多公式，建议在横向分辨率大于 1920px 的屏幕上阅读。
目录
0. 写在前面
这篇文章应该是博客里到目前为止最枯燥的一篇了……
写这篇文章的原因有二：一是在我的学位论文中有一部分内容需要涉及图像噪声模型，正好这几天读完 Radiometric CCD Camera Calibration and Noise Estimation ，权当记录备忘；二是因为搜了一圈发现目前网上基本没有对相机传感器噪声的数学模型进行详尽介绍的中文资料，所以也算是做一点微小的工作，填补这一块的空白吧 🙄
实际上，对于 CV 领域内大部分应用来说，$f=\mu+\epsilon$ 这种经典的模型已经足以对图像的像素值进行表征了。但是对于一些更偏工程、需要与相机硬件直接打交道的应用，往往还需要对成像过程中各阶段引入的噪声进行更加具体的建模和分析，从而将其最大程度抑制。正因如此，我默认阅读这篇文章的读者都具备辐射度学、相机成像原理、图像处理以及统计学的相关基础知识。
此外，图像噪声这一领域并不是我的研究方向，文中难免存在一些瑕疵或错误，欢迎指正。
1. 相机响应值构成模型
在分析信号中的噪声之前，我们当然首先需要搞明白信号究竟是如何构成的。
不管是 CMOS 也好 CCD 也好，图像传感器的本质都是将光信号转换为电信号。对于图像传感器上的单个感光单元（有叫作 collection site 的，或者 photosensory cell，或者 sensor element，总之指的是「单个像素」）来说，其感光后由光电效应释放出的电子个数 $I$ 可以表示为
\begin{equation}
I = T\int_\lambda\int_x\int_y\,E(x,y,\lambda)\,S(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(x,y)$ 表示感光单元内 的灵敏度空间波动，$q(\lambda)$ 为相机的光电转换函数。
从量纲角度分析会更好理解一些：积分时间 $T$ 的单位为 $s$，光谱辐照度 $E$ 的单位为 $W/(m^2\cdot{}nm)$，$S$ 无量纲，光电转换函数 $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 = T\overline{S}A\int_\lambda\,E(\lambda)\,q(\lambda)\,\textrm{d}\lambda
\label{eq:electrons_num_simple}
\end{equation}
其中 $\overline{S}$ 表示 $S(x,y)$ 在单个感光单元内的期望值（均值），$A$ 表示单个感光单元的有效感光面积。
此外，感光单元内光谱辐照度和被摄物体表面的光谱辐亮度（spectral radiance）之间还满足如下关系：
\begin{equation}
E(\lambda) = \frac{\pi}{4F^2}UL(\lambda)t(\lambda)
\label{eq:radiance2irradiance}
\end{equation}
其中 $F$ 是相机的相对孔径（也就是摄影领域中的那个光圈数 F），$L(\lambda)$ 为物方平面上被摄物体表面的光谱辐亮度，$t(\lambda)$ 为整个相机系统的光谱透过率函数（包括镜头、低通滤镜、Bayer filters 等等），$U$ 是一个由 vignetting、shading 等因素导致的空间非均匀性调制函数，无量纲。
式 \eqref{eq:electrons_num_simple} 的原理不展开说了，感兴趣的话可以简单参考一下下图。
单透镜理想光学系统成像模型 所以，联立式 \eqref{eq:electrons_num_simple} 与 \eqref{eq:radiance2irradiance}，有：
\begin{equation}
I = \frac{\pi}{4F^2}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$ 是某个像素对应的主光线相对于光轴的偏离角。这个照度的四次方衰减定律一般也称为自然渐晕（natural vignetting）。
然而，对于一个真实的相机系统来说，CMOS 平面上的照度非均匀性可不仅仅是 cos 四次方衰减那么简单，必然还包含了各种复杂的原因，比如光学渐晕（optical vignetting）、镜头设计时人为加入的照度调制等等（参考下图）。所以为了不失一般性（好吧说实话，其实是因为这种复杂的非均匀性基本不可能给出解析表达式），我直接简单粗暴地把 $cos^4$ 以及其他的非均匀调制统一丢到这个 $U$ 里去了。
Luminance shading 在这之后，感光单元内产生的电荷信号将经由放大电路（调 ISO 其实就是在调节放大电路的增益系数）进行放大，并经 ADC 转换为数字信号缓存在 RAM 中，待后续存储或调用：
\begin{equation}
D = round\left(\frac{g^\prime{}I + V_\mathrm{offset}}{\eta}\right)
\label{eq:amplify}
\end{equation}
其中 $g^\prime$ 表示模拟电路的增益倍数，$V_\mathrm{offset}$ 是人为加入的一个偏置电压，通常可以通过这种方式来避免小于零的输出信号被 ADC 截零从而改变噪声分布的对称性（一些相机的暗电流就是这么来的），$\eta$ 是 ADC 中的量化步长（quantization step），单位为 $volts/DN$（DN 就是 digital number，表示数字信号的数值，理论上是无量纲的，但是为了表述方便通常会人为给它加上一个量纲），与相机所使用的位深（bit-depth，平时我们说一幅图像是8位的或者12位的或者14位的，指的就是这个位深）有关，$round(\cdot)$ 则表示 ADC 的取整过程。
最后这个输出的数字信号 $D$，我习惯称之为相机的原始响应值 （raw response），也就是相机 raw 图像中可以直接读出的那个数字大小。一些资料里也习惯把它称作 DN (digital number) 或 ADU (analog-to-digital units)。
感光单元释放出的电荷如何进入 ADC 这一细节有待确定。老一些的相机，尤其是 CCD 相机，似乎会为传感器阵列的每一行配备一个寄存器，每个寄存器接收一整行感光单元传来的电荷后再统一传给放大器，各行间并行处理；但是对于新一些的 CMOS，甚至会为每一个像素单独配置一个寄存器和放大器，所有信号的处理以像素为单位并行，从而大幅缩短两次采集之间的间隔。此外，现在的一些相机还会为传感器配置两组具有不同电容的感光单元电荷输出电路，从而在电信号进入放大电路之前就有两组电压值可供切换，即所谓的
Dual Native ISO 技术。
最后，联立式 \eqref{eq:formation} 和 \eqref{eq:amplify}，我们就得到了完整的光信号到数字信号的理想转换模型，或者说是理想的数字信号构成模型 （response formation model）：
\begin{equation}
D = round\left[g\,TS_rAU\int_\lambda\,L(\lambda)\,t(\lambda)\,q(\lambda)\,\textrm{d}\lambda + D_\mathrm{offset}\right]
\label{eq:response_model}
\end{equation}
显然，这里有 $g=g^\prime/\eta$，$D_\mathrm{offset}=V_\mathrm{offset}/\eta$，我把它们分别叫作综合增益 与数字偏置量 （$D_\mathrm{offset}$ 基本就是相机的暗电流数值啦）。
2. 图像噪声来源
第一节中讲的一大堆内容都是理想情况下的成像过程，而在实际情况中，以上光信号 $\rightarrow$ 电信号 $\rightarrow$ 数字信号的过程中必然存在着噪声。
由于图像噪声中既存在空间无关噪声也存在空间相关噪声，所以在对图像传感器的噪声进行建模之前，我们需要把各感光单元所处的位置坐标也考虑进来。将式 \eqref{eq:response_model} 扩展到整个传感器平面上，可以写作〔其实就是在式 \eqref{eq:response_model} 的每个符号后面加上位置坐标 $(i,j)$ 而已~〕：
\begin{equation}
\begin{split}
D(i,j) = round\bigg[&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\bigg.\\[.5em]
\bigg.& + D_\mathrm{offset}(i,j)\bigg]
\end{split}
\label{eq:spatial_response_model}
\end{equation}
接下来，我们来看一下相机成像过程中存在的几类主要噪声：
● 像素响应非均匀性
由于制造加工过程中的一些原因，传感器各个感光单元之间在有效感光面积、光电转换效率、非均匀性、滤色片光谱透过率等指标上往往存在一定的差异。这种像素之间 的差异通常被称为像素响应非均匀性（pixel response non-uniformity, PRNU）。严格意义上 PRNU 应该属于「系统误差」而不是「噪声」，但是出于表述方便的需求我仍然把它划分到传感器噪声的范畴里啦。
稍微做一些简化，如果不考虑波长 $\lambda$ 上的偏差，那么对于某一个感光单元，可以使用非均匀性分布函数
\begin{equation}
K(i,j) = \frac{g(i,j)}{g_0}\frac{\overline{S}(i,j)}{\overline{S}_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$ 表示一个理想的、无任何制造误差的感光单元所对应的综合增益，$\overline{S}_0$、$A_0$、$t_0$、$q_0$ 同理。在现有的半导体制造工艺下， $K$ 的数值仅在单位强度附近很小的一段区间内波动，因此有
\begin{equation}
\overline{E}(K) = 1\,,\qquad\overline{\mathrm{var}}(K)\ll 1
\label{eq:prnu_prop}
\end{equation}
这里的 $\overline{E}(\cdot)$ 和 $\overline{\mathrm{var}}(\cdot)$ 分别表示在空间维度上计算均值 （mean）与样本方差 （sample variance），下同。
💡 各个像素的 $K(i,j)$ 其实在相机出厂时就已经确定了，并不是一个随机变量，所以只能计算均值与样本方差，不存在数学期望或者方差之说；而后文中的一些噪声属于随机变量，因此可以计算数学期望与方差。期望与均值、方差与样本方差之间的区别需要特别注意一下。
把式 \eqref{eq:prnu} 代入式 \eqref{eq:spatial_response_model} 后，有
\begin{equation}
D(i,j)=round\left[g_0K(i,j)U(i,j)I_0+D_\mathrm{offset}(i,j)\right]
\label{eq:spatial_response_model_with_prnu}
\end{equation}
其中
\begin{equation}
I_0 = \frac{\pi{}T}{4F^2}\overline{S}_0A_0\int_\lambda\,L(\lambda)t_0(\lambda)q_0(\lambda)\,\mathrm{d}\lambda
\label{eq:I0}
\end{equation}
表示一个理想的、无任何制造误差的感光单元理论上释放出的电子数 。这个 $I_0$ 在后文中会多次出现。
● 热噪声
光电转换器件的热效应将产生一定数量的服从泊松分布的自由电子，从而使信号出现一定程度的波动，因此此类噪声也被称为热噪声（thermal noise）$N_{Th}$ 。热噪声的强度仅取决于曝光时间 $T$ 以及环境温度，而与光电效应产生的电子数 $I$ 无关。
热噪声是一个随机变量，根据泊松分布的性质有
\begin{equation}
E(N_{Th}) = \mathrm{var}(N_{Th}) = \mu_{Th}
\label{eq:thermal_noise_prop}
\end{equation}
其中 $\mu_{Th}$ 是一个与曝光时间及环境温度有关的变量，$E(\cdot)$、$\mathrm{var}(\cdot)$ 分别表示计算数学期望 （expectation）与方差 （variance），下同。
把热噪声考虑进来后，式 \eqref{eq:spatial_response_model_with_prnu} 就变成了
\begin{equation}
D(i,j)=round\big\{g_0\left[K(i,j)U(i,j)I_0+N_{Th}(i,j)\right]+D_\mathrm{offset}(i,j)\big\}
\label{eq:spatial_response_model_with_thermal_noise}
\end{equation}
● 散粒噪声
根据量子理论，从感光单元中释放出的电子在数量上存在一个随机的涨落，该涨落即为散粒噪声（photon shot noise）$N_S$。散粒噪声是由大量单个事件的统计不确定性引起的，其中主要包括输入光子散粒噪声、光生电流散粒噪声与热效应散粒噪声。散粒噪声的存在使得感光单元最终释放出的有效电子数服从泊松分布：$(I^\ast+N_S)\sim{}Poisson(I^\ast)$，其中 $I^\ast=KUI_0+N_{Th}$ 表示在不考虑散粒噪声的情况下该感光单元应释放出的电子个数。
由于泊松分布有 $E(x)=\mathrm{var}(x)=\lambda$ 的性质，那么对于加入散粒噪声后的信号 $(I^\ast+N_S)$ 而言，有
\begin{equation}
E(I^\ast+N_S) = \mathrm{var}(I^\ast+N_S) = I^\ast = KUI_0+N_{Th}
\label{eq:signal_with_shot_noise_charact}
\end{equation}
因此对于散粒噪声自身而言，有
\begin{equation}
E(N_S) = 0\,,\qquad{}\mathrm{var}(N_S) = I^\ast = KUI_0+N_{Th}
\label{eq:shot_noise_charact}
\end{equation}
把散粒噪声考虑进来后，式 \eqref{eq:spatial_response_model_with_thermal_noise} 就变成了
\begin{equation}
\begin{split}
D(i,j)=round\Big\{&g_0\big[K(i,j)U(i,j)I_0+N_{Th}(i,j)+N_S(i,j)\big]\Big.\\[.6em]
\Big.& + D_\mathrm{offset}(i,j)\Big\}
\end{split}
\label{eq:spatial_response_model_with_shot_noise}
\end{equation}
● 固定模式噪声
对于 CMOS 来说，其各个感光单元所对应的模拟电路中有可能存在一定的噪声，此类噪声通常被称为固定模式噪声（fixed pattern noise, FPN）$N_{FP}$。固定模式噪声与像素响应非均匀性类似，均属于系统误差，两者的区别在于，固定模式噪声是由电路的制造误差引起的，与光电效应无关，且通常呈现出比较明显的空间分布规律。
把固定模式噪声考虑进来后，式 \eqref{eq:spatial_response_model_with_shot_noise} 就变成了
\begin{equation}
\begin{split}
D(i,j)=round\Big\{&g_0\big[K(i,j)U(i,j)I_0+N_{Th}(i,j)+N_S(i,j)+N_{FP}(i,j)\big]\Big.\\[.6em]
\Big.& + D_\mathrm{offset}(i,j)\Big\}
\end{split}
\label{eq:spatial_response_model_with_fpn}
\end{equation}
● 读出噪声
信号放大单元把模拟电压信号进行放大的过程中会引入一定的读出噪声（read-out noise）$N_R$。读出噪声的期望为零，其方差与放大单元的增益系数 $g$ 成线性正相关。
把读出噪声考虑进来后，式 \eqref{eq:spatial_response_model_with_fpn} 就变成了
\begin{equation}
\begin{split}
D(i,j)=round\bigg\{&g_0\Big[K(i,j)U(i,j)I_0+N_{Th}(i,j)+N_S(i,j)+N_{FP}(i,j)\Big.\bigg.\\[.2em]
\bigg.\Big.& + N_R(i,j)\Big] + D_\mathrm{offset}(i,j)\bigg\}
\end{split}
\label{eq:spatial_response_model_with_readout_noise}
\end{equation}
● 量化误差
由于 ADC 输出的数字信号需要以整数的形式存储在寄存器中，因此这一连续信号转化为离散信号的过程中必然会引入一定的量化误差（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{}\mathrm{var}(N_Q) = \tfrac{1}{12}
\label{eq:quantization_error_charact}
\end{equation}
把量化误差考虑进来后，式 \eqref{eq:spatial_response_model_with_readout_noise} 就变成了
\begin{equation}
\begin{split}
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]\\[.8em]
& + D_\mathrm{offset}(i,j) + N_Q(i,j)
\end{split}
\label{eq:noise_model}
\end{equation}
这个就是带有噪声项的 数字相机原始响应值构成模型。再次注明一下，$I_0$ 代表的是一个完美的感光单元理论上释放出的电子数〔见式 \eqref{eq:I0}〕。
显然，由于热噪声、散粒噪声、读出噪声、量化误差均为随机变量，因此式 \eqref{eq:noise_model} 中的原始响应值 $D$ 也是一个随机变量。
把以上六种噪声或者误差的统计特性进行了一个简单的整理，如下表所示。
噪声来源 符号 分布类型 期望/样本均值 方差/样本方差 像素响应非均匀性 $K$ – $\overline{\mu}_K=1$ $\overline{\sigma}_K^2\ll{}1$ 固定模式噪声 $N_{FP}$ – $\overline{\mu}_{FP}$ 未知 $\overline{\sigma}_{FP}^2$未知 热噪声 $N_{Th}$ 泊松分布 $\mu_{Th}$ 未知 $\sigma_{Th}^2$ 未知 散粒噪声 $N_S$ 泊松分布（有效信号+散粒噪声） $\mu_S=0$ $\sigma_S^2=KUI_0+N_{Th}$ 读出噪声 $N_R$ 未知 $\mu_R=0$ $\sigma_R^2$ 未知 量化误差 $N_Q$ 均匀分布 $\mu_Q=0$ $\sigma_Q^2=\tfrac{1}{12}$
3. 噪声估计与误差修正
如果只是要了解图像噪声的来源和大致特性的话，看完上面一节就够了。这一节谈一下如何对传感器的噪声进行参数估计，以及如何利用估计出的噪声参数对像素响应非均匀性（PRNU）以及固定模型噪声（FPN）这两类系统误差进行修正。
为什么我们需要大费周章地对图像传感器的噪声和误差进行分析呢？首先，在进行对成像精度有很高要求的图像处理时，我们希望能够把相机当作严格的测量仪器 而不是一个简单的光信号传感器。既然是测量仪器，那么系统的可靠性和测量结果的可复现性自然是一个重要的性能评价指标，因此系统误差的标定和修正就成为了一个必不可少的步骤。其次，如果不通过建模的手段对噪声和误差的参数进行估计，那么当相机的拍摄参数改变时，或者当场景发生变化时（因为噪声与电子数有关，而电子数与场景自身的辐射特性有关），已有的标定信息将不可避免地失效。因此这一「建模+参数估计」的过程实际上是将相机的标定流程从离线切换至了在线，从而允许我们更加灵活地对系统进行修正。
不过呢，因为推导过程需要引入传感器平面的位置坐标，所以这一节里的各种公式可能看起来并不是那么友好 😳
这一节中所有与空间位置有关的参数都将显式地用 $(i,j)$ 进行表示。若不带该空间坐标则表示该参数为一位置无关的常量。
3.1 一点准备工作：噪声的分解
为了后续推导的需要，我们需要把式 \eqref{eq:noise_model} 中的原始响应值 $D(i,j)$ 拆解为两部分：第一部分 $\mu(i,j)$ 表示其中的直流项（也就是信号的期望），第二部分 $N(i,j)$ 表示由随机变量构成的噪声项：
\begin{equation}
D(i,j) = \mu(i,j) + N(i,j)
\label{eq:decomposition0}
\end{equation}
\begin{equation}
\left\{
\begin{split}
\mu(i,j) & = g_0\Big[K(i,j)U(i,j)I_0 + N_{FP}(i,j) + \mu_{Th}\Big] + D_\mathrm{offset}(i,j)\\[.5em]
N(i,j) & = g_0\Big[N_{Th}(i,j) – \mu_{Th} + N_S(i,j) + N_R(i,j)\Big] + N_Q(i,j)
\end{split}
\right.
\label{eq:decomposition}
\end{equation}
显然，$N(i,j)$ 是一个期望为零 的随机变量，且根据前文，它的方差可以表示为
\begin{equation}
\sigma_N^2(i,j) = g_0^2\left[\sigma_S^2(i,j)+\sigma_R^2\right] + \sigma_Q^2
\label{eq:noise_variance}
\end{equation}
💡 式 \eqref{eq:noise_variance} 中第二个等式里，热噪声 $N_{Th}$ 的方差正好等于它的期望（因为是泊松分布），因此两项相减恰好为零；读出噪声的方差 $\sigma_R^2$ 是一个全局常数，与感光单元所处的位置无关，因此可以省略位置坐标 $(i,j)$。
另外，既然都写到这里了，那就顺带推导一下信噪比的计算公式吧。
信噪比（signal-noise-ratio）定义为信号强度的平方 $\mu^2$ 与方差 $\sigma^2$ 之比 ，因此有
\begin{equation*}
\begin{split}
S/N & = 10\log_{10}\left[\frac{E\left(D^2\right)}{E\left(\sigma_N^2\right)}\right] = 20\log_{10}\left[\frac{\mu}{E(\sigma_N)}\right]\\[1em]
& = 20\log_{10}\Bigg\{\frac{g_0\left(KUI_0 + N_{FP} + \mu_{Th}\right)}{E\left[\sqrt{g_0^2\left(\sigma_S^2+\sigma_R^2\right)+\sigma_Q^2}\,\right]}\Bigg\}
\end{split}
\label{eq:snr}
\end{equation*}
计算信噪比的时候不妨忽略空间非均匀性调制函数 $U$ 和小量 $\sigma_Q^2$，也不考虑 PRNU 和 FPN 两种系统误差，并假设暗电流偏置为零，此时，上式可以简化为
\begin{equation*}
\begin{split}
S/N & = 20\log_{10}\bigg(\frac{I_0 + \mu_{Th}}{\sqrt{I_0 + \mu_{Th} +\sigma_R^2}}\bigg)
\end{split}
\label{eq:snr0}
\end{equation*}
由此我们可以得出那个众所周知的结论：要提高图像的信噪比，只能增加传感器上收集的电子数。调节增益 $g_0$ 或者一味地增加位深都是没有任何帮助的（cf. Noise, Dynamic Range, and Bit Depth ）。而要增加电子数，要么延长积分时间 $T$，要么改进光电转换函数 $q$ 使其具有更高的转换效率。
另外，从信噪比公式其实也可以得出结论：对于处理 raw 图像的同学来说，前期拉高 ISO（增大 $g_0$）和后期调节图像亮度其实并没有本质区别。
进一步地，根据噪声与电子数 $I$ 之间的相关性，式 \eqref{eq:decomposition} 中的噪声项 $N(i,j)$ 可以再被拆解为与电子数相关噪声 $N_I(i,j)$ 以及与电子数无关噪声 $N_C(i,j)$：
\begin{equation}
\left\{
\begin{split}
N_I(i,j) & = g_0N_S(i,j)\\[.8em]
N_C(i,j) & = g_0\left[N_{Th}(i,j) + N_R(i,j) – \mu_{Th}\right] + N_Q(i,j)
\end{split}
\right.
\label{eq:noise_decomposition}
\end{equation}
根据散粒噪声的特性，有：
\begin{equation}
\begin{split}
\sigma_I^2(i,j) & = g_0^2\sigma_S^2(i,j)\\[.8em]
& = g_0^2\left[K(i,j)U(i,j)I_0(i,j) + N_{Th}(i,j)\right]
\end{split}
\label{eq:photon_dependent_noise}
\end{equation}
同时，根据前文还有
\begin{equation}
\sigma_C^2(i,j) = g_0^2\left(\mu_{Th} + \sigma_R^2\right) + \frac{1}{12}
\label{eq:photon_independent_noise}
\end{equation}
正如上面提到的，热噪声的期望与读出噪声的方差均与感光单元所处的位置无关，因此 $\mu_{Th}$ 与 $\sigma_R^2$ 可省略位置坐标。此外，需要注意一下，$N_C(i,j)$ 的表达式 \eqref{eq:noise_decomposition} 中减号右边的那个 $\mu_{Th}$ 是一个常数项，所以在计算方差时得到的其实就是零。
3.2 整体噪声估计
在对图像传感器的系统误差进行估计之前，我们需要先对整体噪声水平 $N(i,j)$ 进行参数估计，而 $N(i,j)$ 的参数当然指的就是它的方差 $\sigma_N^2(i,j)$ 啦。
正如渐晕图片 里显示的那样，即使对于完全均匀的场景，传感器不同区域也会对应于不同的照度，从而对应不同的释放电子数。对于位于传感器的中心与边缘的两个感光单元来说，两者有可能相差数倍以上。因此从这一小节开始，对于实际噪声的参数估计，若无特别说明，均指对于传感器某一特定区域 进行计算。若要对整个传感器进行标定，可先分区域处理（$5\times5$ 或者 $7\times7$ 一般就够了），最后再对各区域的参数进行插值从而得到全局的噪声参数估计。
对一台相机在原始响应值上的最终噪声进行估计通常有两种方法，一是采集两幅 raw 图像 $D_1(i,j)$ 及 $D_2(i,j)$，并在其差值图像的基础上进行噪声方差估计；二是采集多幅图像，然后在均值图像的基础上进行估计。第一种方法涉及的计算稍微麻烦一些，所以我就只介绍第二种方法啦。
在对整体噪声方差 $\sigma_N^2(i,j)$ 进行估计之前，我们需要先计算一下期望值 $\mu(i,j)$。显然，一种最简单的方法就是在时域 上对带噪声的信号进行平滑。
假如我们对着同一固定场景使用完全相同的参数拍摄 $n$ 张图像，那么有
\begin{equation}
\left\{\begin{split}
D_1(i,j) & = \mu_1(i,j) + N_1(i,j)\\[.5em]
& \vdots\\[.5em]
D_{n}(i,j) & = \mu_{n}(i,j) + N_{n}(i,j)
\end{split}\right.
\label{eq:multiple_images}
\end{equation}
由于噪声项 $N(i,j)$ 具有 $E\left[N(i,j)\right] = 0$ 的特性，而当图像数量 $n$ 足够大时，我们可以使用样本均值作为数学期望的一致性近似 。因此，$\mu(i,j)$ 的一致性估计 结果也就可以通过在时域上计算均值的方法得到：
\begin{equation}
\begin{split}
\widehat{\mu}(i,j) & \approx \widetilde{E\,}\left[\mu(i,j)\right]\\[.5em]
& = \tfrac{1}{n}\sum\nolimits_{p=1}^{n}\mu_p(i,j)\\[.5em]
& = \tfrac{1}{n}\sum\nolimits_{p=1}^{n}D_p(i,j) – \tfrac{1}{n}\sum\nolimits_{p=1}^{n}N_p(i,j)\\[.5em]
& \approx \tfrac{1}{n}\sum\nolimits_{p=1}^{n}D_p(i,j) – E\left[N(i,j)\right]\\[.5em]
& = \tfrac{1}{n}\sum\nolimits_{p=1}^{n}D_p(i,j)
\end{split}
\label{eq:mu_estimator}
\end{equation}
式中的 $\widetilde{E\,}(\cdot)$ 表示在时域中计算均值（temporal mean），下同。
接着，我们使用 $\widetilde{\mathrm{var}}\left[N(i,j)\right]$ 作为 $\mathrm{var}\left[N(i,j)\right]$ 的近似，从而得到 $\sigma_N^2(i,j)$ 的一致性估计：
\begin{equation}
\begin{split}
\widehat{\sigma}_N^2(i,j) & \approx \widetilde{\mathrm{var}}\left[N(i,j)\right]\\[.8em]
& = \widetilde{\mathrm{var}}\left[D(i,j)\right]\\[.8em]
& = \tfrac{1}{n}\sum\nolimits_{p=1}^{n}\left[D_p(i,j)-\widehat{\mu}(i,j)\right]^2
\end{split}
\label{eq:overall_noise_estimator}
\end{equation}
式中的 $\widetilde{\mathrm{var}}(\cdot)$ 表示在时域中计算样本方差（temporal sample variance），下同。
我使用这种方法对 Nikon D3x 和 SONY A7 两台相机 G 通道的整体噪声方差 $\sigma_N^2(i,j)$ 进行了估计，结果如下图所示。实验时我使用均匀白板作为拍摄对象，并根据 ISO 的不同对光源的强度进行了调解以确保 raw 图像中的响应值尽可能占满有效位深，同时将拍摄张数 $n$ 设定为16。
Nikon D3x 的 $\widehat{\sigma}_N^2(i,j)$
SONY A7 的 $\widehat{\sigma}_N^2(i,j)$
3.3 综合增益系数的估计
接下来我们讨论一下如何通过拍照的方式对相机的综合增益 $g_0$ 进行估计。这部分的内容会稍微 tricky 一些~
由于 $N_I$ 与 $N_C$ 之间不存在相关性，因此整体随机噪声的方差 $\sigma_N^2$ 可被分解为两部分：
\begin{equation}
\begin{split}
\sigma_N^2(i,j) & = \sigma_I^2(i,j) + \sigma_C^2(i,j)\\[.6em]
& = g_0^2\left[K(i,j)U(i,j)I_0 + N_{Th}(i,j)\right] + \sigma_C^2
\end{split}
\label{eq:overall_noise_decomposition}
\end{equation}
接下来要划重点啦~ 对于式 \eqref{eq:decomposition} 中的 $\mu(i,j)$，我们可以对等号两侧同时计算局部空间均值 （local spatial mean），那么有
\begin{equation}
\begin{split}
\overline{E}^w\left[\mu(i,j)\right] & = g_0\left\{\overline{E}^w\left[K(i,j)U(i,j)\right]I_0 + \overline{E}^w\left[N_{FP}(i,j)\right] + \mu_{Th}\right\}\\[.6em]
& + \overline{E}^w\left[D_\mathrm{offset}(i,j)\right]
\end{split}
\label{eq:mu_local_spatial_mean}
\end{equation}
这里的符号 $\overline{E}^w\left[x(i,j)\right]$ 表示以 $(i,j)$ 坐标为中心，在 $w\times{}w$ 的窗口内对矩阵 $\mathbf{x}$ 进行空间均值滤波 （在实际操作时其实就是用一个所有元素都相同的 $w\times{}w$ 核对 $\mathbf{x}$ 进行卷积）。
下面简单说一下这一步空间均值滤波的用处。我们知道，非均匀性函数 $U(i,j)$ 的空间变化整理来说是十分缓慢的（看上图 ），因此当窗口尺寸 $w$ 不是太大时 ，可认为该窗口内的 $U(i,j)$ 是一个常数，不妨就叫作 $\overline{U}^w$ 吧。另一方面，当 $w$ 比较大时 ，这一窗口内包含的像素数就足够多了，因此在一个足够大的窗口内对 $K(i,j)$、$N_{Th}(i,j)$ 以及 $D_\mathrm{offset}(i,j)$ 计算均值，得到的结果也可以认为是位置无关 的，并且根据上文有
\begin{equation}
\left\{
\begin{split}
\overline{E}^w\left[K(i,j)\right] & \approx 1\\[.4em]
\overline{E}^w\left[N_{Th}(i,j)\right] & \approx \mu_{Th}\\[.4em]
\overline{E}^w\left[D_\mathrm{offset}(i,j)\right] & \approx \mu_\mathrm{offset}\\
\end{split}
\right.\qquad{}\mathrm{if}\ w \mathrm{\ is\ large\ enough}
\label{eq:local_spatial_means}
\end{equation}
这里的 $\mu_\mathrm{offset}$ 表示响应值偏置量的期望值（暗电流！），可以从传感器的参数手册里查到，或者用 Dcraw 之类的软件获得。
同样的道理，对于式 \eqref{eq:overall_noise_decomposition}，我们也可以对等号两侧同时同时计算局部空间均值，有
\begin{equation}
\overline{E}^w\left[\sigma_N^2(i,j)\right] = g_0^2\left\{\overline{E}^w\left[K(i,j)U(i,j)I_0\right] + \overline{E}^w\left[N_{Th}(i,j)\right]\right\} + \sigma_C^2
\label{eq:overall_noise_local_spatial_mean}
\end{equation}
所以，只要我们选取一个合适的不大不小的（ 😳 ）窗口尺寸 $w$，然后把式 \eqref{eq:local_spatial_means} 分别代入式 \eqref{eq:mu_local_spatial_mean} 和式 \eqref{eq:overall_noise_local_spatial_mean} 中，就可以得到
\begin{equation}
\left\{
\begin{split}
\overline{\mu}^w(i,j) & = g_0\left[\overline{U}^w(i,j)I_0 + \overline{N}_{FP}^w(i,j) + \mu_{Th}\right] + \mu_\mathrm{offset}\\[.5em]
\overline{\sigma_N^2}^w(i,j) & = g_0\left[\overline{U}^w(i,j)I_0 + \mu_{Th}\right] + \sigma_C^2
\end{split}
\right.
\label{eq:mu_vs_overall_noise}
\end{equation}
观察上式就会发现，如果把 $\overline{\mu}^w(i,j)$ 当作自变量 $x$，把 $\overline{\sigma_N^2}^w(i,j)$ 当作因变量 $y$，那么两者之间其实显然满足 $y=ax+b$ 的直线方程，其中斜率就是 $g_0$，而截距是 $\sigma_C^2 – g_0\mu_\mathrm{offset} – g_0^2\overline{N}_{FP}^w(i,j)$：
\begin{equation}
\overline{\sigma_N^2}^w(i,j) = g_0\overline{\mu}^w(i,j) + \left[\sigma_C^2 – g_0\mu_\mathrm{offset} – g_0^2\overline{N}_{FP}^w(i,j)\right]
\label{eq:mu_vs_overall_noise_formula}
\end{equation}
所以呢，只要我们能够获得足够多的 $(x,y)$ 样本点，那么自然就可以通过直线拟合的方式确定斜率 $g_0$ 了。实际操作时，可以通过改变光源亮度的方式获得 $n$ 组不同的 $D(i,j)$，然后对式 \eqref{eq:mu_estimator} 和式 \eqref{eq:overall_noise_estimator} 进行卷积得到局部均值图像，从而得到 $n$ 组 $\left[\overline{\mu}^w(i,j),\ \overline{\sigma_N^2}^w(i,j)\right]$ 坐标，接着在使用最大似然估计对 $g_0$ 进行计算。
改变光源亮度这一操作对于实验环境的要求比较苛刻，因为大多数标准光源基本上都不允许对亮度进行调节，即使能够调节，也很难保证相对光谱功率不发生变化。此外再澄清一点，图像的张数 $n$ 和窗口尺寸 $w$ 这两个符号在上下文中可能会多次出现，我故意不对它们加以区分，否则公式读起来实在太不友好了…… 大家自行区分吧，它们在不同小节中都是独立存在的。
但是呢，用这种拟合的方式去估计 $g_0$ 会出现一个问题：式 \eqref{eq:mu_vs_overall_noise_formula} 中自变量和因变量都是带有位置坐标 $(i,j)$ 的，因此最终估计出的综合增益 $g_0$ 也必定是位置相关的。而在前文我们说了，$g_0$ 表示的是一个理想的、无任何制造误差的感光单元所对应的综合增益，所以我们还需要对估计出的 $g_0(i,j)$ 求一个平均值以获得最终的综合增益估计值 $\widehat{g}_0$。
对于 Nikon D3x 和 SONY A7 两台相机使用这种方法估计出的综合增益 $g_0$ 如下图所示。实验中我选取的窗口尺寸 $w=15$，光源等级数 $n=8$。
Nikon D3x 的三通道综合增益 $\widehat{g}_0$
SONY A7 的三通道综合增益 $\widehat{g}_0$
3.4 暗电流噪声的估计
在前文中我并没有严格地对「暗电流」这个名词进行定义，这里补充解释一下。
如果我们在完全黑暗（光子数为零）的条件下使用相机进行拍摄，那么式 \eqref{eq:decomposition} 可以改写为
\begin{equation}
\left\{
\begin{split}
\mu_\mathrm{dark}(i,j) & = g_0N_{FP}(i,j) + g_0\mu_{Th} + D_\mathrm{offset}(i,j)\\[.5em]
N_\mathrm{dark}(i,j) & = g_0\Big[N_{Th}(i,j) + N_{S,\mathrm{dark}}(i,j) + N_R(i,j) – \mu_{Th}\Big] + N_Q(i,j)
\end{split}
\right.
\label{eq:decomposition_dark}
\end{equation}
显然，这里的 $\mu_\mathrm{dark}(i,j)$ 对应了相机在未接收到任何光信号 情况下对应的数字响应值期望，因此也通常被称为暗电流（dark current）或者黑电平（black level）。在很多场合我们提起相机的暗电流通常指的是一个具体的数字，比如像下面这张图中 Dcraw 检测到的两台相机 darkness 分别为0和128，但是在需要严格标定的情况下，相机的暗电流噪声应该是一个关于各个像素坐标 $(i,j)$ 的函数，$\mu_\mathrm{dark}(i,j)$，而不是一个简单的常数。在这篇文章中，我认为相机的固定模式噪声 $N_{FP}(i,j)$ 和响应值偏置量 $D_\mathrm{offset}(i,j)$ 共同组成了暗电流噪声，当然，还包括热噪声期望 $\mu_{Th}$ 以及综合增益 $g_0$。
Dcraw 检测出的两台相机的暗电流数值大小 估计暗电流噪声 $\mu_\mathrm{dark}(i,j)$ 当然再简单不过啦，利用 $E\left[N_\mathrm{dark}(i,j)\right]=0$ 这个性质，我们只需要在全黑环境下拍摄 $n$ 幅图像然后取一下均值，就可以得到关于暗电流噪声的一致性估计：
\begin{equation}
\begin{split}
\widehat{\mu}_\mathrm{dark}(i,j) & \approx \widetilde{E\,}\left[\mu_\mathrm{dark}(i,j)\right]\\[.5em]
& = \tfrac{1}{n}\sum\nolimits_{p=1}^{n}D_{\mathrm{dark},p}(i,j) – \tfrac{1}{n}\sum\nolimits_{p=1}^{n}N_{\mathrm{dark},p}(i,j)\\[.5em]
& \approx \tfrac{1}{n}\sum\nolimits_{p=1}^{n}D_{\mathrm{dark},p}(i,j) – E\left[N_\mathrm{dark}(i,j)\right]\\[.5em]
& = \tfrac{1}{n}\sum\nolimits_{p=1}^{n}D_{\mathrm{dark},p}(i,j)
\end{split}
\label{eq:dark_current_estimator}
\end{equation}
下面两幅图是我使用这种方法对两台相机 G 通道估计出的暗电流噪声 $\widehat{\mu}_\mathrm{dark}(i,j)$，这里的拍摄张数 $n$ 我设定为16。由于暗电流噪声中包含了固定模式噪声，因此在下图中其实已经可以明显地看出 $\widehat{\mu}_\mathrm{dark}(i,j)$ 中的 pattern 是存在一定的空间规律的，比如 Nikon D3x 的噪声就呈现出纵向相似性和横向相异性。
Nikon D3x 的 $\widehat{\mu}_\mathrm{dark}(i,j)$
SONY A7 的 $\widehat{\mu}_\mathrm{dark}(i,j)$
3.5 像素响应非均匀性的估计
最后，我们再讨论一下如何对像素响应非均匀性 $K(i,j)$ 进行估计。
这一小节的方法其实和图像处理中的空间滤波有异曲同工之处，因为虽然 $K(i,j)$ 与 $U(i,j)$ 看似结合在了一起共同影响了最终的响应值，但我们可以利用它们空间频率存在差异这一性质对其进行解耦。
为表达方便起见，我首先令 $e(i,j) = K(i,j)U(i,j)I_0$。回想一下前文提到的一些性质，如果我们对 $e(i,j)$ 中各像素在一个 $w\times{}w$ 的邻域内计算均值，那么一方面，当 $w$ 足够小时，可认为在该窗口内的 $U(i,j)$ 为一常数〔因为 $U(i,j)$ 不含高频成分〕，不妨称作 $\overline{U}^w$；另一方面，当 $w$ 足够大时，该窗口内 $K(i,j)$ 的均值必定趋向于1〔因为 $K(i,j)$ 的数值总是在1附近浮动，可参考上文第2节 〕。因此，我们只需要选取一个合适的窗口尺寸 $w$，就可以得到
\begin{equation}
\overline{e}^w(i,j) = \overline{E}^w\left[e(i,j)\right] \approx \overline{U}^w(i,j)I_0
\label{eq:prnu_spatial_mean}
\end{equation}
利用这个近似，就可以对 $K(i,j)$ 进行估计啦：
\begin{equation}
K(i,j) \approx \frac{K(i,j)U(i,j)}{\overline{U}^w(i,j)} = \frac{e(i,j)}{\overline{e}^w(i,j)}
\label{eq:prnu_estimator}
\end{equation}
因此，和第3.3 小节中 $g_0$ 的估计方法一样，只要获得足够多的 $x$ 和 $y$ 坐标〔也就是 $\overline{e}^w(i,j)$ 和 $e(i,j)〕$，我们就可以通过直线拟合的方式估计出一条斜率为 $K(i,j)$ 且经过原点的直线。实际操作时，仍然可通过改变光源亮度的方式获得 $n$ 组 $\left[\overline{e}^w(i,j),\ e(i,j)\right]$ 坐标。
为了避免拍摄物体表面存在的缺陷或污点对直线拟合造成干扰，可以在每次切换光源时稍微移动被摄物体，这样即使在 $n$ 组样本点中存在个别异常点，也很容易发现并剔除，从而保证结果的鲁棒性。
下面再讨论一下 $e(i,j)$ 的计算。
联立式 \eqref{eq:decomposition0} 和式 \eqref{eq:decomposition}，有
\begin{equation}
D(i,j) = g_0\left[e(i,j) + N_{FP}(i,j) + \mu_{Th}\right] + D_\mathrm{offset}(i,j) + N(i,j)
\label{eq:response_reformula}
\end{equation}
再一次地，利用 $E\left[N(i,j)\right]=0$ 这一性质，我们可以通过在同一条件下拍摄 $n$ 幅图像并计算均值的方式消去 $N(i,j)$：
\begin{equation}
\begin{split}
\widetilde{E}\left[D(i,j)\right] & = \tfrac{1}{n}\sum\nolimits_{p=1}^n{}D_p(i,j)\\[.6em]
& \approx g_0\left[e(i,j) + N_{FP}(i,j) + \mu_{Th}\right] + D_\mathrm{offset}(i,j)
\end{split}
\label{eq:response_without_noise}
\end{equation}
所以，联立式 \eqref{eq:response_without_noise} 和式 \eqref{eq:decomposition_dark} 就可以得到 $e(i,j)$ 的估计：
\begin{equation}
\widehat{e}(i,j) \approx \frac{\tfrac{1}{n}\sum\nolimits_{p=1}^n{}D_p(i,j)-\widehat{\mu}_\mathrm{dark}(i,j)}{\widehat{g}_0}
\label{eq:e_estimator}
\end{equation}
下面两幅图是我对两台相机 G 通道中像素响应非均匀性 $K(i,j)$ 估计出的结果，这部分实验中我选用了 $n=16$ 以及 $w=15$。
Nikon D3x 的 $\widehat{K}(i,j)$
SONY A7 的 $\widehat{K}(i,j)$
References
