np.mean计算噪声均值输出异常的问题排查与原因咨询

背景说明 #

需对锐利边缘图像信号，估算其噪声标准差与噪声均值，首先按以下步骤生成带噪边缘图像：

import numpy as np
from scipy import special

定义边缘描述函数：

esf_p = lambda x, sigma, mu, I1, I0: I0  + (0.5)*(I1 - I0)*(1 + special.erf((x-mu)/(np.sqrt(2)*sigma)))

生成边缘图像：

x = np.mgrid[0:100,0:100][1]
im = esf_p(x, 0.62, 50 ,100, 0)

无噪图像效果如下：

给图像添加噪声：

noisefunc = np.random.normal(loc = 0,  scale = esf_p(x, 0.62, 50, 10,  10))
im_noise = im + noisefunc

带噪图像效果如下：

问题复现 #

按以下方式估算噪声：

noise_est = im_noise - np.mean(im_noise, axis=0)
noise_mean_estimate = np.mean(noise_est, axis=0)

绘制图像最后一行像素的噪声对比：

plt.plot(noise_est[-1], label = 'estimated')
plt.plot(noisefunc[-1], label= 'original noise')
plt.legend()

估算噪声与ground truth噪声对比效果如下：

可见二者匹配度较高，但后续噪声均值计算结果完全不符合预期：

noise_mean = np.mean(noise_est, axis=0)
plt.plot(np.sqrt(100)*noise_mean, label = 'mean of noise from noise estimate')
plt.legend()

估算噪声的均值曲线如下：

该结果远低于预期（设置的噪声标准差为10），乘以sqrt(100)是为了抵消100行平均的缩放效应，还原噪声本身的均值量级。
若计算时排除最后一行像素，结果发生13个数量级的大幅变化：

noise_mean = np.mean(noise_est[:-1, :], axis=0)
plt.plot(np.sqrt(99)*noise_mean, label = 'mean of noise from noise estimate')
plt.legend()

排除最后一行的估算噪声均值曲线如下：

此时结果有所改善，但仍比ground truth低一个数量级。
直接计算ground truth噪声的均值无该问题：

plt.plot(np.sqrt(100)*np.mean(noisefunc, axis=0), label = 'mean of ground truth noise')
plt.plot(noisefunc[-1], label = 'noise in last row of pixels')
plt.legend()

ground truth噪声均值与最后一行噪声对比效果如下：

异常规律验证 #

测试估算噪声均值的标准差随参与计算行数的变化规律：

test_std = np.zeros(100)
for i in range(100):
    test_std[i] = np.sqrt(i)*np.std(np.mean(noise_est[:i+1, :], axis=0), axis=0)
plt.plot(test_std, label = 'std of estimated noise mean')
plt.legend()

估算噪声均值的标准差随行数变化曲线如下：

行数较少时标准差接近预期的10，随行数增加反而大幅下降，100行时仅为1e-13量级，与“行数越多精度越高”的预期完全相反。
ground truth噪声的均值标准差随行数变化时保持稳定：

test_std = np.zeros(100)
for i in range(100):
    test_std[i] = np.sqrt(i)*np.std(np.mean(noisefunc[:i+1, :], axis=0), axis=0)
plt.plot(test_std, label = 'std of ground truth noise function')
plt.legend()

ground truth噪声均值的标准差随行数变化曲线如下：

疑问 #

1. np.mean出现该异常的原因是什么？
2. 为什么仅删除一行数据就会导致噪声均值结果发生数量级的变化？
3. 为什么匹配度很高的估算噪声和ground truth噪声，用np.mean计算的结果差异巨大？

该问题的核心是噪声估算逻辑存在数学错误，与np.mean的实现无关：

• 你执行noise_est = im_noise - np.mean(im_noise, axis=0)操作时，本质是对每一列的所有像素做了去均值处理。根据均值运算的基本性质，处理后得到的noise_est每一列的算术平均值天然为0，你观测到的1e-13量级的计算结果就是浮点数运算的正常误差，完全符合理论预期。
• 删除最后一行后计算结果发生数量级变化的原因是：当你仅用前99行计算列均值并做去均值处理时，前99行的noise_est均值依然为0，但你此时计算的均值没有被全列的去均值操作强制归零，对应结果是采样误差的量级，所以数值会大幅上升。
• 单行列的估算噪声与原始噪声匹配度高的原因是：列去均值操作仅移除了列内的直流分量（包含固定的边缘信号梯度和噪声的列均值），单一行噪声的波动形态没有被破坏，所以视觉上匹配度高，但列层面的均值已经被操作强制修正，自然和ground truth的均值计算结果存在巨大差异。

修正方案：若要正确估算噪声，需要先拟合得到无噪的边缘信号im_fit，再通过noise_est = im_noise - im_fit计算噪声估计值，后续再做均值、标准差计算即可得到符合预期的结果。

内容的提问来源于stack exchange，提问作者wilziax