遥感数字图像处理-07-图像滤波

  1. 1. 概述
  2. 2. 图像噪声、卷积和滤波
    1. 2.1. 空间域滤波
    2. 2.2. 频率域滤波
    3. 2.3. 噪声
      1. 2.3.1. 噪声种类
      2. 2.3.2. 噪声特征
      3. 2.3.3. 遥感图像中常见的噪声
    4. 2.4. 卷积
  3. 3. 图像平滑
    1. 3.1. 均值滤波
    2. 3.2. 中值滤波
    3. 3.3. 比较
    4. 3.4. 高斯低通滤波
    5. 3.5. 梯度倒数加权法
    6. 3.6. 选择式掩膜平滑
  4. 4. 图像锐化
    1. 4.1. 线性锐化滤波器
      1. 4.1.1. 梯度法
      2. 4.1.2. 罗伯特梯度
      3. 4.1.3. Prewitt 和Sobel梯度
      4. 4.1.4. Sobel梯度
      5. 4.1.5. 比较
    2. 4.2. 算子
      1. 4.2.1. Laplacian 算子
      2. 4.2.2. 定向检测
      3. 4.2.3. Canny算子
      4. 4.2.4. 锐化图像的生成
  5. 5. 频率域滤波
    1. 5.1. 滤波器
    2. 5.2. 自定义滤波器
    3. 5.3. 同态滤波
  6. 6. 空间域滤波和频率域滤波

概述

图像滤波也是一种图像增强方法

基于像素的空间邻域关系使用卷积计算或频率解析的方法进行计算,以达到抑制图像噪声、增加图像对比度、增强图像边缘、提取图像特征的图像处理方法

图像噪声、卷积和滤波

空间域滤波

常用方法:卷积运算

缺点:随着采用的模板窗口的扩大,运算量越来越大(例如:7*7窗口)

解决方法:可在频率域中通过简单的乘法计算来实现

频率域滤波

  • 图像中像素值随位置的变化可以用频率来描述,即空间频率。
  • 边缘、线条、噪声等特征,如河流、湖泊的边界、道路、差异较大的地表交界处等具有高的空间频率,即在较短的像素距离内灰度值重现的频率大。
  • 均匀分布的地物或大面积的稳定结构,如植被类型一致的平原、沙漠、海面等具有低的空间频率,即在较长的像素距离内灰度值重现的频率小。
  • 图像中的高频变化对应于局部变化一一·从像素到像素的变化,低频变化对应于区域变化一一从图像的一部分到另一部分的变化。
  • 频率域滤波是在频率域中对输入图像应用滤波函数(滤波器)来改进输出图像的处理方法。
  • 在频率域滤波中,保留图像的低频部分抑制高频部分的处理称为低通滤波,起到平滑的作用。
  • 保留图像的高频部分而削弱低频部分的处理称为高通滤波,起到锐化的作用。

噪声

噪声种类

  • 按产生原因
    • 外部噪声
    • 内部噪声
  • 从统计理论观点
    • 平稳噪声
    • 非平稳噪声
  • 从噪声幅度分布形态
    • 高斯型
    • 瑞利型
  • 按产生过程
    • 量化噪声
    • 椒盐噪声

噪声特征

  • 噪声可看做是对亮度(像素值)的干扰用n(x,y)n(x,y)表示
  • 随机性,用随机过程来描述
  • 由于分布函数或密度函数很难测出或描述,常用统计特征(如均值、方差、总功率)等来描述噪声

加性噪声模型

g(x,y)=f(x,y)+n(x,y)g(x, y)=f(x, y)+n(x, y)

式中:f(x,y)f(x, y)为理想图像,n(x,y)n(x, y)为噪声,g(x,y)g(x, y)为输出图像

乘性噪声模型

g(x,y)=f(x,y)[1+n(x,y)]=f(x,y)+f(x,y)n(x,y)\begin{aligned} g(x, y) &=f(x, y)[1+n(x, y)] \\ &=f(x, y)+f(x, y) n(x, y) \end{aligned}

对于加性噪声而言,n(x,y)n(x, y)和图像光强大小无关,通常表现为高斯噪声或脉冲噪声
对于乘性噪声而言,n(xy)n(x,y)和图像光强大小相关,随亮度的大小变化而变化
乘性噪声或许是图像中最普遍的噪声,其模型和分析计算都比较复杂。通常总是假定信号和噪声互相独立,然后通过对图像做对数变换,将乘性噪声当做加性噪声来处理。

遥感图像中常见的噪声

  • 高斯噪声:噪声的像素值分布可以使用高斯概率密度来描述
  • 脉冲噪声:脉冲噪声又称椒盐噪声,构成噪声的像素在空间分布上是随机的,但像素值是固定的。随机改变一些像素值,在二值图像上表现为使一些像素点变白(用b表示),一些像素点变黑(用a表示)。
  • 周期噪声:获取过程中受成像设备影响产生的。这是唯一的一种空间依赖型噪声。周期噪声可通过频率域滤波进行压抑。

脉冲噪声的概率密度函数可表示为:

p(z)={paz=apbz=b0 其他 p(z)=\left\{\begin{array}{ll} p_{a} & z=a \\ p_{b} & z=b \\ 0 & \text { 其他 } \end{array}\right.

PaPaPbPb相等时,脉冲噪声值将类似于随机分布在图像上的胡椒和盐粉微粒。由于这个原因,双极脉冲噪声也称为椒盐噪声。

空间分布都是随机,但高斯噪声的值成高斯分布。椒盐噪声的值是常数

卷积

i1=n1n1j1=m1m1fi+i1,j+j1wij\sum_{i 1=-n 1}^{n 1} \sum_{j 1=-m 1}^{m 1} f_{i+i 1, j+j 1} \bullet w_{i j}

邻域定义、邻域大小、卷积核的值

图像平滑

  • 受传感器和大气影响图像上会存在噪声
  • 表现为:亮点或亮度过大的区域
  • 图像平滑目的:
    • 抑制噪声、改善图像质量
    • 突出图像主体
    • 增强面的信息

均值滤波

均值滤波是最常用的线性低通滤波器,它均等地对待邻域中的每个像素。
取邻域像素值的平均作为该像素的新值
计算公式为:

g(x,y)=k=1nl=1m(f(k,l)h(k,l))g(x, y)=\sum_{k=1}^{n} \sum_{l=1}^{m}(f(k, l) \cdot h(k, l))

  • 优点:算法简单,计算速度快
  • 缺点:造成图像模糊,削弱了边缘和细节信息
  • 算法改进:引进阈值T

g(x,y)={g(x,y), 当 f(x,y)g(x,y)>Tf(x,y), 当 f(x,y)g(x,y)Tg(x, y)=\left\{\begin{array}{l} g(x, y), \text { 当 }|f(x, y)-g(x, y)|>T \\ f(x, y), \text { 当 }|f(x, y)-g(x, y)| \leq T \end{array}\right.

中值滤波

将窗口内的所有像素值按大小排序后,取中值作为中心像素的新值。

在对图像应用滤波器进行过滤时,边界问题是一个需要处理的问题。一般来说,有3种处理的方法。

  1. 不做边界处理
    不对图像的边界作任何处理,在对图像进行滤波时,滤波器没有作用到图像的四周,因此图像的四周没有发生改变。
  2. 填充0
    对图像的边界做扩展,在扩展边界中填充0,对于边长为2k+1的方形滤波器,扩展的边界大小为k,若原来的图像为[m, n],则扩展后图像变为[m+2k, n+2k]。进行滤波之后,图像会出现一条黑色的边框。
  3. 填充最近像素值
    扩展与 填充0 的扩展类似,只不过填充0的扩展是在扩展部分填充0,而这个方法是填充距离最近的像素的值。

比较

随机噪声(高斯噪声、脉冲噪声):均值滤波优于中值滤波
脉冲噪声干扰的椒盐噪声:中值滤波有效

对于椒盐噪声,中值滤波效果比均值滤波效果好。
原因

  • 椒盐噪声是幅值近似相等但随机分布在不同位置上,图像中有干净点也有污染点。中值滤波是选取合理的邻近像素值来替代噪声点,所以只适合于椒盐噪声的去除,不适合高斯噪声的去除。
  • 因为噪声的均值不为0,所以均值滤波不能很好地去除噪声点。

对于高斯噪声,均值滤波效果比中值滤波效果好
原因:

  • 高斯噪声是幅值近似正态分布,但分布在每点像素上。
  • 因为图像中的每点都是污染点,所中值滤波选不到合适的干净点。
  • 因为正态分布的均值为0,所以根据统计数学,均值可以消除噪声。(注意:实际上只能减弱,不能消除**(因为均值滤波是通过滑动窗口进行局部“均值”进行滤波,而局部的高斯噪声均值不等于0
    )**)

高斯低通滤波

在一个二维空间,高斯分布函数定义为:

f(x,y)=12πσ2ex2+y22σ2f(x, y)=\frac{1}{2 \pi \sigma^{2}} e^{-\frac{x^{2}+y^{2}}{2 \sigma^{2}}}

按照MATLAB软件文档,空间域的高斯低通滤波器可定义为:

hg(x,y)=ex2+y22σ2h_{g}(x, y)=e^{-\frac{x^{2}+y^{2}}{2 \sigma^{2}}}

归一化后为:

H(x,y)=hg(x,y)hg(x,y)H(x, y)=\frac{h_{g}(x, y)}{\sum h_{g}(x, y)}

高斯滤波器对于高斯噪声的滤除非常有效,平滑程度是由σ来控制,σ值越大,平滑程度越高,滤波后的图像越模糊。

梯度倒数加权法

原理:在离散图像内部相邻区域的变化大于区域内部的变化,在同一区域中中间像素的变化小于边沿像素的变化。

目标:平滑其贡献重要来自区域内部的像素,平滑后的图像边沿和细节不会受到明显损害

建立归一化的权重矩阵作为模板

h=[w(k1,l1)w(k1,l)w(k1,l+1)w(k,l1)w(k,l)w(k,l+1)w(k+1,l1)w(k+1,l)w(k+1,l+1)]h=\left[\begin{array}{lll} w(k-1, l-1) & w(k-1, l) & w(k-1, l+1) \\ w(k, l-1) & w(k, l) & w(k, l+1) \\ w(k+1, l-1) & w(k+1, l) & w(k+1, l+1) \end{array}\right]

规定,w(k,l)=1/2w(k,l)=1/2,其余8个加权系数之和为1/2
定义其它权重矩阵元素为

w(k,l)=1f(x+k,y+l)f(x,y)w(k, l)=\frac{1}{|f(x+k, y+l)-f(x, y)|}

通过下式将总和归一化为1/2

w(k,l)=w(k,l)2k=11l=11w(k,l)w(k, l)=\frac{w(k, l)}{2 \sum_{k=-1}^{1} \sum_{l=-1}^{1} w(k, l)}

选择式掩膜平滑

  • 凡含有尖锐边缘的区域,方差必定比平缓区域大
  • 完成滤波操作,不破坏区域边界的细节

取5×5窗口,在窗口内 以中心像素为基准点,制作4个五边形、4个六边形、1个边长为3的正方形共9个掩模

计算各掩模的均值aia_{i}及方差kik_{i}

ai=[j=1Qfj(x+k,y+l)]/Qki=j=1Q{[fj(x+k,y+l)]2ai2}\begin{aligned} a_{i} &=\left[\sum_{j=1}^{Q} f_{j}(x+k, y+l)\right] / Q \\ k_{i} &=\sum_{j=1}^{Q}\left\{\left[f_{j}(x+k, y+l)\right]^{2}-a_{i}^{2}\right\} \end{aligned}

式中,i=129i=1,2,…,9
QQ为掩模对应的像素个数;
kkll为掩模内像素相对于中心像素的位移量。

kik_{i}排序,最小方差kmink_{min}所对应的掩模的灰度级均值aia_{i}作为图像f的平滑输出gg
用同样的方法作用于每一像素,即可完成整个图像的平滑。

图像锐化

  • 突出图像中的地物边缘、轮廓或线状目标
  • 提高了边缘与周围像素之间的反差,也被称为边缘增强

线性锐化滤波器

原始图像减去低通图像
梯度图像=原始图像-均值滤波的图像

梯度法

梯度实际上反映了相邻像素之间灰度的变化率,图像中的边缘,例如河流、湖泊的边界、道路等处灰度的变化率较大,因此在边缘处一定有较大的梯度值;而大面积的平原、海面灰度变化较小,一定具有较小的梯度值;对于灰度级为常数的区域,梯度值为0。

因此,以梯度值替代像素的原灰度值生成梯度图像,在梯度图像上梯度值较大的部分就是边缘。

非线性的锐化滤波器
梯度

gradf(x,y)=[fxfy]=[f(x,y)xf(x,y)y]\operatorname{gradf}(x, y)=\left[\begin{array}{c} f_{x}^{\prime} \\ f_{y}^{\prime} \end{array}\right]=\left[\begin{array}{l} \frac{\partial f(x, y)}{\partial x} \\ \frac{\partial f(x, y)}{\partial y} \end{array}\right]

梯度的模为

gradf(x,y)=((x,y)x)2+((x,y)y)2|\operatorname{grad} f(x, y)|=\sqrt{\left(\frac{\partial(x, y)}{\partial x}\right)^{2}+\left(\frac{\partial(x, y)}{\partial y}\right)^{2}}

离散的数字图像

gradf(x,y)f(x,y)f(x+1,y)+f(x,y)f(x,y+1)|g r a d f(x, y)| \cong f(x, y)-f(x+1, y)|+| f(x, y)-f(x, y+1) \mid

grad f=fh1+fh2\text {grad } f=f^{*} h 1+f^{*} h 2

罗伯特梯度

采用交叉差分的方法

lgradf(x,y)=f(x,y)f(x+1,y+1)+f(x+1,y)f(x,y+1)\operatorname{lgradf}(x, y)|=| f(x, y)-f(x+1, y+1)|+| f(x+1, y)-f(x, y+1) \mid

Prewitt 和Sobel梯度

Prewitt算法扩大了模板,从2×2扩大到3×3来进行差分

f(x1,f(x,yf(x+1,yy1)1)1)f(x1,yf(x,y)f(x+1,y)f(x1,f(x,yf(x+1,yy+1)+1+1)\begin{array}{|c|c|c|} \hline f(x-1, & f(x, y & f(x+1, y \\ y-1) & -1) & -1) \\ \hline f(x-1, & & \\ y & f(x, y) & f(x+1, y) \\ \hline f(x-1, & f(x, y & f(x+1, y \\ y+1) & +1 & +1) \\ \hline \end{array}

Sobel梯度

  • Sobel梯度是在Prewitt算法的基础上,对4邻域采用加权方法进行差分,因而对边缘的检测更加精确
  • 在Prewitt和模板中,h1突出水平方向的地物,h2突出垂直方向的地物。在图像分割中,经常使用sobel模板。
  • 需要注意的是,梯度模板的总和为0值,而前面用于图像平滑的模板总和不为0。上述模板不适用于含有大量噪声的图像。

比较

对于具有明显分布方向的地物,Sobel算子锐化效果更好一些。

算子

算子是方便计算而使用的符号(运算符),指从一个函数空间到另一个函数空间(或它自身)的映射。算子接受一个函数并按照特定的规则得到另一个函数

常用的算子
&-D&微分算子,gradgrad梯度算子,2\nabla^{2}拉普拉斯算子

g(x,y)=Df(x,y)+2f(x,y)g(x, y)=D f(x, y)+\nabla^{2} f(x, y)

g(x,y)=2f(x,y)g(x, y)=\nabla^{2} f(x, y)指图像的拉普拉斯算子,即使用拉普拉斯的卷积核对f(x,y)进行卷积计算,其结果是图像g(x,y)。

Laplacian 算子

Laplacian算子是线性二阶微分算子,即

2f=2fx2+2fy2\nabla^{2} f=\frac{\partial^{2} f}{\partial x^{2}}+\frac{\partial^{2} f}{\partial y^{2}}

对于离散的数字图像,推导出Laplacian算子的表达式为

2f(x,y)=f(x+1,y)+f(x1,y)+f(x,y+1)+f(x,y1)4f(x,y)\nabla^{2} f(x, y)=f(x+1, y)+f(x-1, y)+f(x, y+1)+f(x, y-1)-4 f(x, y)

即取某像素的上下左右四个相邻像素的值相加的和减去该像素的四倍,作为该像素新的灰度值。

因为仅用拉普拉斯算子进行处理,会把图像中的噪声放大。通常是先将图像用高斯函数进行平滑处理后,再进行拉普拉斯滤波。这样处理符合了人类视觉的生理特性,这种滤波器称为高斯-拉普拉斯滤波器(LOG)

定向检测

为了有目的的提取某一特定方向的边缘或线性特征,可以选用特定的模板进行卷积运算。

Canny算子

边缘是图像的重要特征,基于梯度可以提取图像的边缘信息。canny算子对于受白噪声影响的阶跃型边缘是最优的,是当前影响最大的图像边缘检测方法。该算子的设计目标为:
①不丢失重要的边缘,不产生虚假的边缘,减少噪声的影响
②检测到的边与实际的边具有最小的偏差,强调结果的正确性
③将多个响应降低为单边缘的响应,限制单个边缘点对亮度变化的定位。

锐化图像的生成

计算各个像素的梯度值,根据不同的需求生成锐化图像。(拉普拉斯算子、sobel算子等生成的都是梯度图)

 锐化图像=原图像(或平滑后图像)+ / - 梯度图像 \text { 锐化图像=原图像(或平滑后图像)+ / - 梯度图像 }

也可以组合使用不同的规则产生锐化后图像:
以梯度值代替其原灰度值,即

g(x,y)=gradf(x,y)g(x, y)=|\operatorname{grad} f(x, y)|

适当选取 TT ,使梯度值\geqslant TT 的各点的灰度等于该点的梯度值,其他则保留原灰度值,形成背景,即

g(x,y)={gradf(x,y)gradf(x,y)Tf(x,y) 其他 g(x, y)=\left\{\begin{array}{l|l} |\operatorname{grad} f(x, y)| & \operatorname{grad} f(x, y) \geqslant T \\ f(x, y) & \text { 其他 } \end{array}\right.

根据需要指定一个灰度级 LG,L_{G}, 例如,令 LG=255L_{G}=255 。以 LGL_{G} 表示边缘,其他保留原背景值, 即

g(x,y)={LGgradf(x,y)Tf(x,y) 其他 g(x, y)=\left\{\begin{array}{ll} L_{G} & |\operatorname{grad} f(x, y)| \geqslant T \\ f(x, y) & \text { 其他 } \end{array}\right.

注:可利用密度分割确认阈值T
指定一个灰度级 LBL_{B} 表示背景,例如,令 LB=0,L_{B}=0, 形成黑背景,保留边缘梯度变化,即

g(x,y)={gradf(x,y)+cgradf(x,y)TLB 其他 g(x, y)=\left\{\begin{array}{ll} |\operatorname{grad} f(x, y)|+c & \operatorname{grad} f(x, y) \geqslant T \\ L_{B} & \text { 其他 } \end{array}\right.

式中,c 为常数,取 0 或更大的值,用于更好的区分梯度和背景。
将梯度与灰度图像分别以灰度级 LGL_{G}LBL_{B} 表示,例如,255 表示边缘, 0 表示背景,形成二值图像,即

g(x,y)={LGgradf(x,y)TLB 其他 g(x, y)=\left\{\begin{array}{ll} L_{G} & |\operatorname{grad} f(x, y)| \geqslant T \\ L_{B} & \text { 其他 } \end{array}\right.

注:梯度大的保留,小的压制,相当于提取了边

频率域滤波

基于频率域图像使用滤波器一直特定范围频率的图像处理方法

g(u,v)=h(u,v)F(u,v)g(u, v)=h(u, v)^{\star} F(u, v)

滤波的关键是正确选择了滤波器并且确定了合适的通或阻的频率。与滤波类型相比,选择合适的截止频率更难一些。遥感图像与物理信号不同,是不同频率地物的混杂,分离出特定地物的频率并据此进行滤波,往往需要多次的实验才能得到满意的结果。但是,对于图像中的周期性信号(主要是图像扫描产生的噪声)的去除,频率域滤波更容易得到满意的结果。

按照信号处理理论,根据滤除的频率的特征,滤波有三种:

  • 低通滤波:削弱或抑制高频部分而保留低频部分
  • 高通滤波:突出图像的边缘和轮廓,进行图像锐化
  • 带通滤波:仅仅保留指定频率范围的滤波
  • 带阻滤波:仅仅保留指定频率范围的滤波

滤波器

概念:对特定频率的频点或该频点意外的频率进行滤除,得到一个特定频率的信号,或消除一个特定频率信号的硬件设备或函数

理想滤波器:能够完全截止和滤波的滤波器

自定义滤波器

在实际工作中,最常用的是用户定义滤波器,可以根据频率的分布灵活选择滤波的频率范围,包括线、矩形、椭圆、多边形等。

滤波器选择和设计是关键,如果有偏差或错误,就不会得到合适的结果。要进行多次的尝试比较才能最后确定需要的滤波器。

尚没有明确的指导如何确定滤波器的半径
频率是由整个图像的数据计算产生的,与局部没有对应关系,截止频率(半径)并不对应于特定的局部区域。

遥感图像处理中,特别是光学遥感图像处理中,频率域滤波不如空间域滤波的应用范围广。

同态滤波

减少低频(背景)增加高频,从而减少光照变化并锐化边缘或细节的图像滤波方法

一幅图像f(x,y)f(x, y)可以用照射分量和反射分量来模拟,即

f(x,y)=i(x,y)×r(x,y)f(x, y)=i(x, y) \times r(x, y)

空间域滤波和频率域滤波