一种高分辨率遥感图像信噪比曲线自适应获取方法

摘要

本发明一种高分辨率遥感图像信噪比曲线自适应获取方法，包括以下步骤：对整幅图像进行傅里叶变换，通过不同频率区间的掩膜，求得功率谱密度趋于稳定的频率区间；掩膜后的噪声频谱傅里叶反变换到空间域，对信号图像和噪声图像进行非重叠的小块划分，分别求得每个小块的平均值和方差；去除信号平均值和噪声方差的野值，逐步迭代拟合得到噪声方差与信号平均值的线性关系，求得直线的斜率和截距；整个灰度值范围内的信号除以噪声标准偏差得到信噪比曲线。本发明针对单幅遥感图像本身，高精度获取整个灰度值范围内的信噪比曲线。解决了光学遥感卫星信噪比在轨监测需要的均匀区域数量非常多，耗费时间很长，难以实现业务化监测的难题。

权利要求

1.一种高分辨率遥感图像信噪比曲线自适应获取方法，其特征在于步骤如下：

1)在高分辨率遥感图像g(x,y)上建立二维空间直角坐标系；以g(x,y)的左上角为原点、X轴正方向向下、Y轴正方向向右建立二维空间直角坐标系，x＝0,1,…,M-1，y＝0,1,…,N-1，g(x,y)的大小为M×N个像素；M、N为正整数；

2)对高分辨率遥感图像g(x,y)进行傅里叶变换，在得到的傅里叶变换上建立二维频率直角坐标系；以傅里叶变换的左上角为原点、U轴正方向向下、V轴正方向向右建立二维频率直角坐标系，将g(x,y)的傅里叶变换的零频中心变换到二维频率直角坐标系上的(M/2,N/

2)，得到零频中心平移后的傅里叶变换，表示为G(u,v)，u,v为频率变量，u＝0,1,…,M-1，v＝0,1,…,N-1；

3)在二维频率直角坐标系中构造频谱掩膜核函数一K1(u,v)和频谱掩膜核函数二K2(u,v)；其中，K1(u,v)表示为：该函数在区间[R1,Rmax]范围内的值为1，称为非零区间一，在其它区间该函数的值为0；

K2(u,v)表示为：

该函数在区间[R2,Rmax]范围内的值为1，称为非零区间二，在其它区间该函数的值为0；

非零区间一[R1,Rmax]和非零区间二[R2,Rmax]中的R1、R2、Rmax的关系为：

4)对上述零频中心平移后的傅里叶变换G(u,v)分别与K1(u,v)和K2(u,v)进行点乘运算，得到掩膜频谱一和掩膜频谱二，分别将掩膜频谱一和掩膜频谱二的高频中心平移到二维频率直角坐标系的(M/2,N/2)，得到高频中心平移后的掩膜频谱一N1(u,v)和掩膜频谱二N2(u,v)，分别对N1(u,v)和N2(u,v)进行傅里叶反变换，将得到的二维空间直角坐标系上的反变换结果分别作为掩膜后的噪声图像一n1(x,y)和掩膜后的噪声图像二n2(x,y)；

5)对高分辨率遥感图像g(x,y)进行分块，块大小为m×m个像素，共得到 个图像块，计算所有图像块的灰度平均值mg(x,y)(i)和灰度标准方差σg(x,y)(i)，表示上取整；

分别对掩膜后的噪声图像一n1(x,y)和掩膜后的噪声图像二n2(x,y)进行分块，块大小为m×m个像素，都得到 个噪声图像块，计算所有噪声图像块的噪声均值和噪声方差

将噪声方差 与 分别乘以系数λ1、λ2，得到转化的噪声方差

和 表示为：

其中，λ1＝(M×N)/Γ1，λ2＝(M×N)/Γ2，Γ1、Γ2分别表示掩膜频谱一和掩膜频谱二中的非零个数；其中m取值为8或16；

根据转化得到的噪声方差 和 则图像块的噪声方差 计算公

式如下：

所述R1、R2、Rmax的取值范围为0.52≤R1≤0.56，0.60≤R2≤0.64，0.66≤Rmax≤0.7；

6)根据上述所得到的所有图像块的灰度平均值mg(x,y)(i)与灰度标准方差σg(x,y)(i)剔除不满足条件的图像块，通过逐步迭代拟合得到信号与噪声方差的直线方程，进而得到整个灰度范围内的信噪比曲线，具体步骤如下：(6.1)遍历所有图像块，判断每个图像块的灰度平均值mg(x,y)(i)与标准方差σg(x,y)(i)是否满足关系式 如果不满足则剔除该图像块，否则保留该图像块；

(6.2)将保留下来的每个图像块的灰度平均值与其对应的噪声图像块的噪声方差作为一个点，所有图像块的灰度平均值与其对应的噪声图像块的噪声方差构成点集其中，点集中的点个数小于划分的图像块数 即(6.3)利用最小二乘法对点集 进行直线拟合，并计算点集中的各

个点到直线的距离，如果距离小于设定的阈值，则保留点集中的该点，否则剔除点集中的该点，该步骤中剔除的点集中的点个数记为η；

(6.4)重复步骤(6.2)、步骤(6.3)，直到剔除的点集中的点个数η＝0，利用点集中的剩余的点拟合直线方程，p、q为拟合系数，则得到信号S与噪声方差σ2的关系式可表示为σ2＝p×S+q，信噪比即信号S与噪声标准方差σ的比值，其计算公式如下：

7)对于量化级数为L的遥感图像，其对应的灰度值范围为[0,2L-1]之间的整数，将每一个灰度值作为信号S，通过公式 计算每个灰度值的信噪比，最终得到整个灰度值范围内的信噪比曲线。

说明书

一种高分辨率遥感图像信噪比曲线自适应获取方法

技术领域

[0001] 本发明涉及一种高分辨率遥感图像信噪比曲线自适应获取方法，属于高分辨率遥感卫星成像质量监测与评价技术领域，特别适用于高分辨率遥感卫星在轨成像质量业务化监测。

背景技术

[0002] 卫星在轨图像的信噪比(SNR)作为光学遥感卫星的关键性能指标之一，特别对于高分辨率遥感卫星来说，信噪比直接影响了卫星的图像质量，因此，能够准确获得卫星图像的信噪比具有重要意义。目前，国内外主要通过铺设靶标法和典型目标法定期测量光学遥感卫星的信噪比。

[0003] (1)铺设靶标法：在地面铺设不同反射率的均匀靶标，在不同的太阳高度角成像条件下，卫星在轨对靶标多次成像，经过多次测量插值拟合之后得到信噪比曲线。但是，该方法不可能做到铺设所有反射率的靶标，因此，得不到整个量化灰度范围内的信噪比曲线；同时，受卫星飞行轨道和天气等条件的限制，很难保证获取靶标图像的有效性，导致了信噪比的测量周期长。

[0004] (2)典型目标法：对具有典型特征的地物目标(比如马路、河流、桥梁、圆屋顶等)进行成像，充分利用这些刃边和点目标信息，类似于布设的靶标，通过特征提取分析成像系统的信噪比。该方法对目标有一定条件的限制，由于地物的非绝对均匀性，因此，比较难找到满足测量要求的典型目标。

[0005] 总之，铺设靶标与典型目标法只能通过有限数量的反射率靶标和多个典型目标计算分析得到几个特定灰度值的信噪比值，难以得到整个量化灰度范围内的信噪比曲线。

发明内容

[0006] 本发明解决的技术问题是：克服现有技术的不足，提供一种高分辨率遥感图像信噪比曲线自适应获取方法，针对任意高分辨率光学遥感卫星图像，不仅能够得到某个灰度值的信噪比，同时可以准确获取整个灰度值范围内的信噪比曲线。解决了目前光学遥感卫星信噪比在轨监测依赖铺设靶标或寻找典型目标，而导致的信噪比监测周期长、难以实现业务化监测以及难以获得信噪比曲线的问题。

[0007] 本发明的技术方案是：一种高分辨率遥感图像信噪比曲线自适应获取方法，步骤如下：

[0008] 1)在高分辨率遥感图像g(x,y)上建立二维空间直角坐标系；以g(x,y)的左上角为原点、X轴正方向向下、Y轴正方向向右建立二维空间直角坐标系，x＝0,1,…,M-1，y＝0,1,…,N-1，g(x,y)的大小为M×N个像素；M、N为正整数；

[0009] 2)对高分辨率遥感图像g(x,y)进行傅里叶变换，在得到的傅里叶变换上建立二维频率直角坐标系；以傅里叶变换的左上角为原点、U轴正方向向下、V轴正方向向右建立二维频率直角坐标系，将g(x,y)的傅里叶变换的零频中心变换到二维频率直角坐标系上的(M/2,N/2)，得到零频中心平移后的傅里叶变换，表示为G(u,v)，u,v为频率变量，u＝0,1,…,M-
1，v＝0,1,…,N-1；

[0010] 3)在二维频率直角坐标系中构造频谱掩膜核函数一K1(u,v)和频谱掩膜核函数二K2(u,v)；其中，K1(u,v)表示为：

[0011]

[0012] 该函数在区间[R1,Rmax]范围内的值为1，称为非零区间一，在其它区间该函数的值为0；

[0013] K2(u,v)表示为：

[0014]

[0015] 该函数在区间[R2,Rmax]范围内的值为1，称为非零区间二，在其它区间该函数的值为0；

[0016] 非零区间一[R1,Rmax]和非零区间二[R2,Rmax]中的R1、R2、Rmax的关系为：

[0017] 4)对上述零频中心平移后的傅里叶变换G(u,v)分别与K1(u,v)和K2(u,v)进行点乘运算，得到掩膜频谱一和掩膜频谱二，分别将掩膜频谱一和掩膜频谱二的高频中心平移到二维频率直角坐标系的(M/2,N/2)，得到高频中心平移后的掩膜频谱一N1(u,v)和掩膜频谱二N2(u,v)，分别对N1(u,v)和N2(u,v)进行傅里叶反变换，将得到的二维空间直角坐标系上的反变换结果分别作为掩膜后的噪声图像一n1(x,y)和掩膜后的噪声图像二n2(x,y)；

[0018] 5)对高分辨率遥感图像g(x,y)进行分块，块大小为m×m个像素，共得到个图像块，计算所有图像块的灰度平均值mg(x,y)(i)和灰度标准方差表示上取整；

[0019] 分别对掩膜后的噪声图像一n1(x,y)和掩膜后的噪声图像二n2(x,y)进行分块，块大小为m×m个像素，都得到 个噪声图像块，计算所有噪声图像块的噪声均值和噪声方差

[0020] 将噪声方差 与 分别乘以系数λ1、λ2，得到转化的噪声方差和 表示为：

[0021]

[0022] 其中，λ1＝(M×N)/Γ1，λ2＝(M×N)/Γ2，Γ1、Γ2分别表示掩膜频谱一和掩膜频谱二中的非零个数；其中m取值为8或16；

[0023] 根据转化得到的噪声方差 和 则图像块的噪声方差σn(x,y)2(i)计算公式如下：

[0024]

[0025] 6)根据上述所得到的所有图像块的灰度平均值mg(x,y)(i)与灰度标准方差σg(x,y)(i)剔除不满足条件的图像块，通过逐步迭代拟合得到信号与噪声方差的直线方程，进而得到整个灰度范围内的信噪比曲线，具体步骤如下：

[0026] (6.1)遍历所有图像块，判断每个图像块的灰度平均值mg(x,y)(i)与标准方差σg(x,y)(i)是否满足关系式 如果不满足则剔除该图像块，否则保留该图像块；

[0027] (6.2)将保留下来的每个图像块的灰度平均值与其对应的噪声图像块的噪声方差作为一个点，所有图像块的灰度平均值与其对应的噪声图像块的噪声方差构成点集(mg(x,y)(j),σn(x,y)2(j))，其中，点集中的点个数小于划分的图像块数 即

[0028] (6.3)利用最小二乘法对点集(mg(x,y)(j),σn(x,y)2(j))进行直线拟合，并计算点集中的各个点到直线的距离，如果距离小于设定的阈值，则保留点集中的该点，否则剔除点集中的该点，该步骤中剔除的点集中的点个数记为η；

[0029] (6.4)重复步骤(6.2)、步骤(6.3)，直到剔除的点集中的点个数η＝0，利用点集中的剩余的点拟合直线方程，p、q为拟合系数，则得到信号S与噪声方差σ2的关系式可表示为σ2＝p×S+q，信噪比即信号S与噪声标准方差σ的比值，其计算公式如下：

[0030]

[0031] 7)对于量化级数为L的遥感图像，其对应的灰度值范围为[0,2L-1]之间的整数，将每一个灰度值作为信号S，通过公式 计算每个灰度值的信噪比，最终得到整个灰度值范围内的信噪比曲线。

[0032] 所述R1、R2、Rmax的取值范围一般为0.52≤R1≤0.56，0.60≤R2≤0.64，0.66≤Rmax≤0.7。

[0033] 本发明与现有技术相比有益效果为：

[0034] (1)本发明提供一种高分辨率遥感图像信噪比曲线自适应获取方法，较传统的遥感卫星依赖铺设靶标或寻找典型目标进行信噪比在轨监测的方法，该方法仅需要单幅遥感图像就可以获得某个灰度值的信噪比，同时可以高精度获取图像整个灰度值范围内的信噪比曲线，实现遥感卫星的信噪比业务化监测。

[0035] (2)本发明提供一种基于高频掩膜的遥感图像信噪比曲线自适应获取方法，通过自适应高频掩膜，除了获取高斯分布的白噪声之外，还能获取与信号相关的泊松分布的散弹噪声。

附图说明

[0036] 图1是本发明高分辨率遥感图像信噪比曲线自适应获取方法流程图；

[0037] 图2是本发明信号能量与噪声示意图；

[0038] 图3是本发明构造的频谱掩膜核函数一和频谱掩膜核函数二的示意图；

[0039] 图4a是本发明实施例采用的遥感图像，图4b是对遥感图像进行模糊的MTF卷积核，图4c是对图4a进行MTF卷积模糊后的图像；

[0040] 图5a是本发明频谱掩膜区间一对应的掩膜核图像，图5b是掩膜后的噪声图像，图5c是频谱掩膜区间二对应的掩膜核图像，图5d是掩膜后的噪声图像；

[0041] 图6是本发明所有图像块的信号均值与噪声方差关系图；

[0042] 图7是本发明基于最小二乘法去野值后的信号均值与噪声方差，以及信号均值与噪声方差拟合得到的直线方程；

[0043] 图8a是本发明实施例中加入随机噪声的100次试验，每次试验的测量值与理论值的2阶范数误差图，图8b与图8c给出了第一次、第二次试验的测量值与理论值曲线图。

具体实施方式

[0044] 本发明一种高分辨率遥感图像信噪比曲线自适应获取方法，具体步骤如图1所示，该方法由以下步骤实现：

[0045] 1、在高分辨率遥感图像g(x,y)上建立二维空间直角坐标系；以g(x,y)的左上角为原点、X轴正方向向下、Y轴正方向向右建立二维空间直角坐标系，x＝0,1,…,M-1，y＝0,1,…,N-1，g(x,y)的大小为M×N个像素；M、N为正整数。

[0046] 2、对高分辨率遥感图像g(x,y)进行傅里叶变换，在得到的傅里叶变换上建立二维频率直角坐标系；以傅里叶变换的左上角为原点、U轴正方向向下、V轴正方向向右建立二维频率直角坐标系，将g(x,y)的傅里叶变换的零频中心变换到二维频率直角坐标系上的(M/2,N/2)，得到零频中心平移后的傅里叶变换，表示为G(u,v)，u,v为频率变量，u＝0,1,…,M-
1，v＝0,1,…,N-1。将得到的傅里叶坐标变换进行零频坐标平移后，为了确保移动后的坐标为整数，要求M和N为偶数。

[0047] 根据遥感图像的退化原理，在轨获得的高分辨率遥感图像g(x,y)由清晰图像f(x,y)与模糊核h(x,y)进行模糊卷积，再加上噪声n(x,y)所得，即：

[0048] g(x,y)＝h(x,y)*f(x,y)+n(x,y)

[0049] 因此，上述退化原理在频域的傅里叶变换关系表示为：

[0050] G(u,v)＝H(u,v)F(u,v)+N(u,v)

[0051] 对于通常的卫星图像，尤其是高分辨率遥感卫星图像，其模糊核H(u,v)的高频部分是非常小的，遥感图像通过模糊核作用过后的H(u,v)F(u,v)在高频部分的能量非常弱，可以忽略不计。如图2所示，在高频R>R0处，噪声部分的能量N(u,v)占主导作用。根据这一点可以设计噪声方差σ2的快速获取算法，从g(x,y)中获取噪声n(x,y)的方差σ2。

[0052] 3、在二维频率直角坐标系中构造频谱掩膜核函数一K1(u,v)和频谱掩膜核函数二K2(u,v)；其中，K1(u,v)表示为：

[0053]

[0054] 该函数在区间[R1,Rmax]范围内的值为1，称为非零区间一，在其它区间该函数的值为0。

[0055] K2(u,v)表示为：

[0056]

[0057] 该函数在区间[R2,Rmax]范围内的值为1，称为非零区间二，在其它区间该函数的值为0。

[0058] 所述频谱掩膜核函数一K1(u,v)的非零区间[R1,Rmax]和频谱掩膜核函数二K2(u,v)的非零区间[R2,Rmax]中的R1、R2、Rmax的关系约束为：

[0059]

[0060] 根据经验值，R1、R2、Rmax的取值范围一般为0.52≤R1≤0.56，0.60≤R2≤0.64，0.66≤Rmax≤0.7。

[0061] 图3是本发明构造的频谱掩膜核函数一和频谱掩膜核函数二的示意图，选取两个高频区间的频谱掩膜核函数，其主要目的是防止选取的不够准确的高频下限值R1或R2会使得G(u,v)·K1(u,v)或G(u,v)·K2(u,v)在频率高频处衰减得不够厉害，噪声部分的能量N(u,v)占主导作用的假设不正确，从而导致获取出来的噪声方差往往要比实际的σ2值要偏大。因此，为了解决这一问题，采取了选择两个不同下界的策略[R1,Rmax]和[R2,Rmax]；而设定上界Rmax的主要目的是为了避免边界效应带来的影响。

[0062] 4、对上述零频中心平移后的傅里叶变换G(u,v)分别与K1(u,v)和K2(u,v)进行点乘运算，得到掩膜频谱一和掩膜频谱二，分别将掩膜频谱一和掩膜频谱二的高频中心平移到二维频率直角坐标系的(M/2,N/2)，得到高频中心平移后的掩膜频谱一N1(u,v)和掩膜频谱二N2(u,v)，分别对N1(u,v)和N2(u,v)进行傅里叶反变换，将得到的二维空间直角坐标系上的反变换结果分别作为掩膜后的噪声图像一n1(x,y)和掩膜后的噪声图像二n2(x,y)，噪声图像的大小都等于M×N个像素。

[0063] 对G(u,v)分别利用频谱掩膜核函数一K1(u,v)和频谱掩膜核函数二K2(u,v)进行点乘运算，即对应元素相乘，表示为G(u,v)·K1(u,v)＝[H(u,v)F(u,v)+N(u,v)]·K1(u,v)和G(u,v)·K2(u,v)＝[H(u,v)F(u,v)+N(u,v)]·K2(u,v)，由于在高频部分H(u,v)F(u,v)对应的能量非常弱，本发明中频域中的圆环频谱掩膜函数半径的区间范围[R1,Rmax]和[R2,Rmax]都是高频区间，所以[H(u,v)F(u,v)]·K1(u,v)≈0和[H(u,v)F(u,v)]·K2(u,v)≈0。本部分也证明了遥感图像在高频部分占主导作用的主要是噪声部分。

[0064] 5、对高分辨率遥感图像g(x,y)进行分块，块大小为m×m个像素，共得到个图像块，计算所有图像块的灰度平均值mg(x,y)(i)和标准方差表示上取整；m一般取值8或16。

[0065] 分别对掩膜后的噪声图像一n1(x,y)和掩膜后的噪声图像二n2(x,y)进行分块，块大小为m×m个像素，都得到 个噪声图像块，计算所有噪声图像块的噪声均值和噪声方差

[0066] 由于噪声方差仅在频谱掩膜区间得到的，在频谱掩膜函数为零的区间被舍弃，因此，需要转化到整个区间，即将噪声方差 与 分别乘以系数λ1、λ2，得到转化的噪声方差 和 表示为

[0067]

[0068] 其中，λ1＝(M×N)/Γ1，λ2＝(M×N)/Γ2，Γ1、Γ2分别表示掩膜频谱一和掩膜频谱二中的非零个数；其中m取值为8或16。

[0069] 根据转化得到的噪声方差 和 则图像块的噪声方差σn(x,y)2(i)计算公式如下：

[0070]

[0071] 为了防止向外插值引起的不稳定现象，对上述计算得到的σn(x,y)2(i)的大小加以2
限制，要求 即最终得到的σn(x,y) (i)被限制
在 的范围内。

[0072] 6、根据上述所得到的所有图像块的灰度平均值mg(x,y)(i)与灰度标准方差σg(x,y)(i)剔除不满足条件的图像块，通过逐步迭代拟合得到信号与噪声方差的直线方程，进而得到整个灰度范围内的信噪比曲线，具体步骤如下：

[0073] (1)遍历所有图像块，判断每个图像块的灰度平均值mg(x,y)(i)与标准方差σg(x,y)(i)是否满足关系式 如果不满足则剔除该图像块，否则保留该图像块；

[0074] (2)将保留下来的每个图像块的灰度平均值与其对应的噪声图像块的噪声方差作为一个点，所有图像块的灰度平均值与其对应的噪声图像块的噪声方差构成点集(mg(x,y)(j),σn(x,y)2(j))，其中，点集中的点个数小于划分的图像块数 即

[0075] (3)利用最小二乘法对点集(mg(x,y)(j),σn(x,y)2(j))进行直线拟合，并计算点集中的各个点到直线的距离，如果距离小于设定的阈值δ，则保留点集中的该点，否则剔除点集中的该点，该步骤中剔除的点集中的点个数记为η，阈值δ取值范围为0.2≤δ≤0.5；

[0076] (4)重复步骤(2)、步骤(3)，直到剔除的点集中的点个数η＝0，利用点集中的剩余的点拟合直线方程，p、q为拟合系数，则得到信号S与噪声方差σ2的关系式可表示为σ2＝p×S+q，信噪比即信号S与噪声标准方差σ的比值，其计算公式如下：

[0077]

[0078] 遥感图像的噪声部分主要包括两部分，一部分为与信号大小相关的满足泊松分布的散弹噪声(系数为p，该系数与信号大小成正比例关系)，另外一部分为与信号大小无关的满足高斯分布的白噪声q。总的噪声方差满足：

[0079] σ2＝p×S+q

[0080] 对光学遥感卫星获取的任意图像，在不需要典型均匀场景的情况下，通过对设定的高频区域进行非零掩膜，而其他频率区间进行零掩膜，利用多个非重叠的分块计算得到与信号无关的高斯分布白噪声和与信号相关的泊松分布散弹噪声，线性拟合分析得到任意一幅图像的信噪比曲线。

[0081] 7、对于量化级数为L的遥感图像，其对应的灰度值范围为[0,2L-1]之间的整数，则每一个灰度值作为信号S，通过公式 计算每个灰度值的信噪比，最终得到整个灰度值范围内的信噪比曲线。

[0082] 例如，对于量化级数为8级的遥感图像，可得到灰度范围为[0,255]区间上的信噪比曲线；对于量化级数为10级的遥感图像，可得到灰度范围为[0,1023]区间上的信噪比曲线。

[0083] 实施例

[0084] 本实施例利用大小为2048×2048个像素的真实高分辨率遥感图像(图4a所示)。为得到不含噪声的遥感图像，需要进行如下处理：首先，对图像通过MTF卷积核模糊(频域的MTF卷积核如图4b所示)，卷积后的图像如图4c所示，大小为2048×2048个像素；然后，对图4c采用4倍的降采样处理，得到512×512个像素的降采样图像，此时的图像可看作是不含噪声的图像。

[0085] 为得到实际噪声的大小，本实施例在降采样图像中人为地加入泊松分布的信号噪声，再加入高斯分布的加性白噪声，得到的图像即为带有噪声的图像(大小为512×512个像素)，并给出理论的信噪比曲线。

[0086] 1、对上述人为加入噪声的图像进行傅里叶变换。

[0087] 2、设定R1＝0.52，R2＝0.60，Rmax＝0.68，则频谱掩膜核函数一的非零区间为[0.52,0.68]，频域中对应的图像如图5a所示，掩膜后的噪声图像如图5b所示。频谱掩膜核函数一的非零区间为[0.60,0.68]，频域中对应的图像如图5c所示，掩膜后的噪声图像如图5d所示。

[0088] 3、对人为加入噪声的图像图像进行分块，每块大小为8×8个像素，共得到64×64块图像，计算所有块图像的灰度平均值。

[0089] 对上述得到的噪声图像分别进行分块，每块噪声图像大小也为8×8个像素，计算所有块噪声图像的均值和方差。所有块图像的信号均值与噪声方差的关系如图6所示。

[0090] 4、根据所有图像块的灰度平均值与标准方差剔除不满足条件的图像块，通过逐步迭代拟合得到灰度平均值与噪声方差的直线方程，即得到信号与噪声方差的线性关系。本实施例取阈值δ＝0.5，图7是本发明基于最小二乘法剔除不满足条件的图像块后的信号的灰度均值与噪声方差，以及信号的灰度均值与噪声方差拟合得到的直线方程：σ2＝0.00168×S+0.36。

[0091] 5、本实施例量化级数为8级的遥感图像，根据拟合得到的信号与噪声方差的关系式σ2＝0.00168×S+0.36，进一步得到[0,255]灰度范围内的信噪比曲线，图8a是本发明实施例中加入随机噪声的100次试验，每次试验的测量值与理论值的2阶范数误差图，图8b与图8c给出了第一次、第二次试验的测量值与理论值曲线图。

[0092] 用以上的方案进行了测量值与理论值的对比分析，测量的曲线与理论上的曲线非常逼近，获取的噪声大小的精度优于90％。

[0093] 本发明未详细说明部分属本领域技术人员公知常识。