专利名称：一种基于边缘和灰度的遥感图像变化检测方法
技术领域：
本发明涉及一种遥感图像变化检测方法，特别是一种基于边缘和灰度的遥感图像 变化检测方法。
背景技术：
文献“一种新的多波段遥感影像变化检测方法，中国图象图形学报，2009， Vol. 14(4)，p572-578”公开了一种基于模糊C均值聚类(FCM)算法的遥感影像变化检测方 法。该方法首先利用一种改进的时间复杂度降低的模糊C均值(FCM)算法进行遥感影像非 监督分类处理，然后引入多波段综合变化掩膜的方法进行变化检测。进行遥感影像非监督 分类处理即进行每个像素是否变化的判决，属于像素级而非特征级的变化检测，而人工目 标包括道路、桥梁、机场、铁路、房屋类建筑等，是一类非常重要线状目标。该方法在对包含 此类目标的遥感图像进行变化检测时，因为没有考虑其线性特征，极易受到噪声的影响，造 成检测到目标不完整，最终变化检测的总体精度仅能达到87. 75%，所以文献公开的方法应 用于处理包含人工目标的遥感图像时具有较大的局限性。
发明内容
为了克服现有的遥感图像变化检测方法精度低的不足，本发明提供一种基于边缘 和灰度的遥感图像变化检测方法，该方法利用基于双边滤波的Carmy算法进行多时相图像 边缘特征提取，然后对灰度差值图像进行OSTU阈值分割和边缘提取，获得灰度特征。再将 所提取的边缘和灰度特征进行综合，检测遥感图像的变化区域。由于充分利用了图像中的 线性特性的同时，利用灰度差值图像弥补了因配准误差造成的断线，可以提高遥感图像变 化检测方法的精度。本发明解决其技术问题所采用的技术方案一种基于边缘和灰度的遥感图像变化 检测方法，其特点是包括下述步骤(a)设去噪算子为Dh，建立双边滤波模型Dh(x) = lw(.x,y)^ti(y)dy(i = ι，2)(ι)
Ω式中，Xti表示ti传感器对于地面上同一位置在ti时获得的遥感图像，Xti (y)表 示图像Xti的位置y处的灰度值，Ω表示积分区域；在低通滤波模型中LPF(x)= jwD(x,y)Xti(y)dy(2)
Ω权函数w(x，y)只和空间分布有关，处理后图像中，每个像素χ处的灰度值LPF(x) 是在其小邻域Ω内的像素上作灰度加权运算得到；式中，wD(x, y)是衡量像素X和y在几 何上邻近关系的权函数；以低通滤波模型为基础，双边滤波模型增加了衡量该像素和其邻域内像素的灰度 值接近程度的权因子BF(x)= jwD(x,y)wR(x,y)Xti(y)dy(3)
Ω
式中， (χ，y)是衡量像素χ和y之间灰度值相近程度的权因子，取
Il^-Jll
wD(x,y) = V-~e 2σ (4)
ZD(x)
(Α-, (X)-Xfi (y)}2
2^(5)
ZR(x)
式中，Ω是积分区域，取为8连通邻域；和σ〗分别是距离方差和灰度方差；Zd (χ)
'D^ ^kj R-
ll^-yllW^yx1Wt
和Zk(χ)分别是归一化因子，计算方法是ZD(x) = \e 2σ° dy ,ZR(x) = fe 2^ dy。
ΩΩ以前相图的边缘为基准，邻域搜索后相图中与之一致的边缘并去除，将两副图像 保留的边缘结果合并，小区域去除后得到边缘变化检测结果；(b)设有像元向量Xtl和Xt2分别为传感器对于地面上同一位置在tl和t2时刻获 得的遥感图像，地面目标的变化表现为向量Xtl和Xt2之间的差值δ = Xtl-Xt2(6)当δ的灰度> (图像均值*比例系数)，则认为当前δ中该像素代表的区域为变 化区域，比例系数设为1. 0 ；当δ的灰度> Thresh时，则认为当前δ中该像素代表的区域为变化区域；OTSU阈值采用
thresh255Thresh = max( ^ % ^ hist(i)^(mx-m2)2)(7)
/=1i=thresh+\计算；式中，hist是图像中灰度为δ的像素数，即直方图中灰度为i时的值；
thresh255
^ i * hist(i) [ i * hist(i)mx = ^-m2 = ^^--(8)
^ hist(J)^ hist{i)
i=li=thresh+\利用形态学滤波对灰度差值图像进行边缘检测，作为灰度特征提取结果；(c)以边缘变化检测结果为基准作边缘跟踪，对于无法闭合的边缘的端点，在灰度 变化结果的边缘图中作邻域搜索，采用“或”逻辑对边缘变化检测结果进行补充，得到闭合 的变化检测边缘及变化区域。
0028]本发明的有益效果是由于利用基于双边滤波的Carmy算法进行多时相图像边缘 特征提取，然后对灰度差值图像进行OSTU阈值分割和边缘提取，获得灰度特征。再将所提 取的边缘和灰度特征进行综合，检测遥感图像的变化区域。由于充分利用了图像中的线性 特性的同时，利用灰度差值图像弥补了因配准误差造成的断线，遥感图像变化检测方法的 精度由背景技术的87. 75%提高到90. 32%，提高了检测准确性。下面结合具体实施方式
对本发明作详细说明。
具体实施例方式
1、边缘特征的提取。
本发明中的边缘检测算法采用改进的Carmy算法。Carmy算子能在噪声抑制和边 缘检测之间取得较好的平衡。传统的Carmy算子采用Gauss滤波在一定程度上减弱了边缘 特性。在本发明中，采用双边滤波(Bilateral filter)而不是Gauss滤波，能更好的获得 图像边缘特性。双边滤波方法如下记去噪算子为Dh，滤波模型可写成如下的一般形式Dh(X)=
Ω式中，Xti表示ti传感器对于地面上同一位置在ti时获得的遥感图像，Xti (y)表 示图像Xti的位置y处的灰度值，Ω表示积分区域。在传统的低通滤波模型中，权函数w(x， y)只和空间分布有关，处理后图像中，每个像素χ处的灰度值LPF(x)是在其小邻域Ω内的 像素上作灰度加权运算得到，表达式为LPF(X)= jwD(x,y)Xli(y)dy(2)
Ω式中，wD(x，y)是衡量像素χ和y在几何上邻近关系的权函数。在此基础上，双边 滤波模型增加了衡量该像素和其邻域内像素的灰度值接近程度的权因子。其表达式如下BF{x)= \wD{x,y)wR{x,y)Xti{y)dy(3)
Ω式中，wK(x，y)是衡量像素χ和y之间灰度值相近程度的权因子，通常取
1 Ik-yllWD(X>y) = -—-e 2σ (4)
Z0(X)
.(Xri(X)-XJy))2^ 2σ (5)
ZR(x)式中，Ω是积分区域，取为8连通邻域，σ 和<7〗分别为距离方差和灰度方差，
Ik-Ml
这里分别取为1和4，Zd(X)和Zk(X)分别为归一化因子，计算方法为Zz3(JC)= P 2σ dy ,
Ω
{X.jxyX^y)}2
ZR(x)= je 2σ1 办。
Ω双边滤波模型在平滑图像过程中，既考虑了几何上的邻近关系，又考虑了灰度上 的相似性，根据邻域内像素灰度值的分布情况，作不同权重的加权运算，很好地体现了算法 各向异性的性质。在图像的边界处，通常的邻域卷积运算会导致图像边界变得模糊。加上 表示灰度值差异的权因(χ，y)后，在图像边界处灰度值有一个比较大的落差，此时权值较 小，在边界一边的灰度分布不会影响到另一边的灰度分布。这样就可以在平滑噪声的同时， 可以有效的保持边界。双边滤波模型最大的优点是其非迭代性，该算法描述简洁易于程序实现，且在很 多情况都达到了很好的去噪效果，这使得去噪模型的计算量有了大幅度的降低，而且该算 法也有很好的并行性。在本发明中对不同时期同一地区拍摄的尺寸与比例均相同的两幅配准图像作边 缘检测，利用边缘检测提取边缘信息，比较并检测出变化的线性特征；同时对两幅差值灰度图像做变化检测得到变化区域的轮廓；将检测的线性特征和变化区域的轮廓特征综合得到 最终的变化结果。2、灰度特征的提取。把图像差值δ作为变化检测的处理图像，从理论上来理解是最直观的，同时它对 于并行处理不同时相之间的变化检测也是非常有效的。在应用这种方法时，必须考虑选取 阈值将有变化像元和无变化像元区分开来，阈值的选择在基于灰度差值的变化检测算法中 非常重要。本发明采用大津法(OTSU)算法进行自适应阈值选择。设有像元向量Xtl和Xt2分别为传感器对于地面上同一位置在tl和t2时可获得的 遥感图像。不考虑其它因素，如大气条件、传感器的影响，即在理想的情况下，地面目标的变 化表现为向量Xtl和Xt2之间的差值δ = Xtl-Xt2(6)通常，tl和t2为不同年的相同或相近的月或日期，以尽量排除季节不同等因素引 起的变化等。算法说明如下以变化前后图像的差值作讨论，当δ的灰度> (图像均值*比例 系数)，则认为当前δ中该像素代表的区域为变化区域，比例系数设为1.0。OTSU方法判断δ变化区域说明如下当δ的灰度> Thresh时，则认为当前δ中 该像素代表的区域为变化区域。OTSU的阈值计算方法是
thresh255Thresh = max( J]/zwi(z)* Σ hist(i)*(jnx-m2)2)(7)
/=i=thresh+l式中，hist是图像中灰度为δ的像素数，即直方图中灰度为i时的值。
thresh255
^z*hist(i)^ /*histif)二·^-m2--(8)
^ hist{i)I hist{i)
/=1i=thresh+l利用形态学滤波对灰度差值图像进行边缘检测，作为灰度特征提取结果。3、边缘和灰度特征的融合。对边缘变化检测结果进行边缘跟踪，扫描无法闭合的边缘端点，方法如下扫描一 个点，如果其周围8个点所有相邻两个点的差的绝对值的和为2X255，则其为端点。以边缘端点为基准，在灰度变化结果的边缘图中做邻域搜索，采用“或”逻辑对边 缘变化检测结果进行补充，得到闭合的变化检测边缘及变化区域。经测试，本发明方法检测遥感图像变化的精度达到了 90. 32%，提高了检测准确性。
权利要求
一种基于边缘和灰度的遥感图像变化检测方法，其特征在于包括下述步骤(a)设去噪算子为Dh，建立双边滤波模型 <mrow><msub> <mi>D</mi> <mi>h</mi></msub><mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo></mrow><mo>=</mo><munder> <mo>∫</mo> <mi>Ω</mi></munder><mi>w</mi><mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo></mrow><msub> <mi>X</mi> <mi>ti</mi></msub><mrow> <mo>(</mo> <mi>y</mi> <mo>)</mo></mrow><mi>dy</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>=</mo> <mn>1,2</mn> <mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>1</mn> <mo>)</mo></mrow> </mrow>式中，Xti表示ti传感器对于地面上同一位置在ti时获得的遥感图像，Xti(y)表示图像Xti的位置y处的灰度值，Ω表示积分区域；在低通滤波模型中 <mrow><mi>LPF</mi><mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo></mrow><mo>=</mo><munder> <mo>∫</mo> <mi>Ω</mi></munder><msub> <mi>w</mi> <mi>D</mi></msub><mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo></mrow><msub> <mi>X</mi> <mi>ti</mi></msub><mrow> <mo>(</mo> <mi>y</mi> <mo>)</mo></mrow><mi>dy</mi><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>2</mn> <mo>)</mo></mrow> </mrow>权函数w(x，y)只和空间分布有关，处理后图像中，每个像素x处的灰度值LPF(x)是在其小邻域Ω内的像素上作灰度加权运算得到；式中，wD(x，y)是衡量像素x和y在几何上邻近关系的权函数；以低通滤波模型为基础，双边滤波模型增加了衡量该像素和其邻域内像素的灰度值接近程度的权因子 <mrow><mi>BF</mi><mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo></mrow><mo>=</mo><munder> <mo>∫</mo> <mi>Ω</mi></munder><msub> <mi>w</mi> <mi>D</mi></msub><mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo></mrow><msub> <mi>w</mi> <mi>R</mi></msub><mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo></mrow><msub> <mi>X</mi> <mi>ti</mi></msub><mrow> <mo>(</mo> <mi>y</mi> <mo>)</mo></mrow><mi>dy</mi><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>3</mn> <mo>)</mo></mrow> </mrow>式中，wR(x，y)是衡量像素x和y之间灰度值相近程度的权因子，取 <mrow><msub> <mi>w</mi> <mi>D</mi></msub><mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo></mrow><mo>=</mo><mfrac> <mn>1</mn> <mrow><msub> <mi>Z</mi> <mi>D</mi></msub><mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo></mrow> </mrow></mfrac><msup> <mi>e</mi> <mrow><mo>-</mo><mfrac> <mrow><mo>|</mo><mo>|</mo><mi>x</mi><mo>-</mo><mi>y</mi><mo>|</mo><mo>|</mo> </mrow> <mrow><mn>2</mn><msubsup> <mi>σ</mi> <mi>D</mi> <mn>2</mn></msubsup> </mrow></mfrac> </mrow></msup><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>4</mn> <mo>)</mo></mrow> </mrow> <mrow><msub> <mi>w</mi> <mi>R</mi></msub><mrow> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo></mrow><mo>=</mo><mfrac> <mn>1</mn> <mrow><msub> <mi>Z</mi> <mi>R</mi></msub><mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo></mrow> </mrow></mfrac><msup> <mi>e</mi> <mrow><mo>-</mo><mfrac> <msup><mrow> <mo>(</mo> <msub><mi>X</mi><mi>ti</mi> </msub> <mrow><mo>(</mo><mi>x</mi><mo>)</mo> </mrow> <mo>-</mo> <msub><mi>X</mi><mi>ti</mi> </msub> <mrow><mo>(</mo><mi>y</mi><mo>)</mo> </mrow> <mo>)</mo></mrow><mn>2</mn> </msup> <mrow><mn>2</mn><msubsup> <mi>σ</mi> <mi>R</mi> <mn>2</mn></msubsup> </mrow></mfrac> </mrow></msup><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>5</mn> <mo>)</mo></mrow> </mrow>式中，Ω是积分区域，取为8连通邻域；和分别是距离方差和灰度方差；ZD(x)和ZR(x)分别是归一化因子，计算方法是以前相图的边缘为基准，邻域搜索后相图中与之一致的边缘并去除，将两副图像保留的边缘结果合并，小区域去除后得到边缘变化检测结果；(b)设有像元向量Xt1和Xt2分别为传感器对于地面上同一位置在t1和t2时刻获得的遥感图像，地面目标的变化表现为向量Xt1和Xt2之间的差值δ＝|Xt1 Xt2|(6)当δ的灰度＞(图像均值*比例系数)，则认为当前δ中该像素代表的区域为变化区域，比例系数设为1.0；当δ的灰度＞Thresh时，则认为当前δ中该像素代表的区域为变化区域；OTSU阈值采用 <mrow><mi>Thresh</mi><mo>=</mo><mi>max</mi><mrow> <mo>(</mo> <munderover><mi>Σ</mi><mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn></mrow><mi>thresh</mi> </munderover> <mi>hist</mi> <mrow><mo>(</mo><mi>i</mi><mo>)</mo> </mrow> <mo>*</mo> <munderover><mi>Σ</mi><mrow> <mi>i</mi> <mo>=</mo> <mi>thresh</mi> <mo>+</mo> <mn>1</mn></mrow><mn>255</mn> </munderover> <mi>hist</mi> <mrow><mo>(</mo><mi>i</mi><mo>)</mo> </mrow> <msup><mrow> <mo>(</mo> <msub><mi>m</mi><mn>1</mn> </msub> <mo>-</mo> <msub><mi>m</mi><mn>2</mn> </msub> <mo>)</mo></mrow><mn>2</mn> </msup> <mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>7</mn> <mo>)</mo></mrow> </mrow>计算；式中，hist是图像中灰度为δ的像素数，即直方图中灰度为i时的值； <mrow><msub> <mi>m</mi> <mn>1</mn></msub><mo>=</mo><mfrac> <mrow><munderover> <mi>Σ</mi> <mrow><mi>i</mi><mo>=</mo><mn>1</mn> </mrow> <mi>thresh</mi></munderover><mi>i</mi><mo>*</mo><mi>hist</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo></mrow> </mrow> <mrow><munderover> <mi>Σ</mi> <mrow><mi>i</mi><mo>=</mo><mn>1</mn> </mrow> <mi>thresh</mi></munderover><mi>hist</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo></mrow> </mrow></mfrac> </mrow> <mrow><msub> <mi>m</mi> <mn>2</mn></msub><mo>=</mo><mfrac> <mrow><munderover> <mi>Σ</mi> <mrow><mi>i</mi><mo>=</mo><mi>thresh</mi><mo>+</mo><mn>1</mn> </mrow> <mn>255</mn></munderover><mi>i</mi><mo>*</mo><mi>hist</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo></mrow> </mrow> <mrow><munderover> <mi>Σ</mi> <mrow><mi>i</mi><mo>=</mo><mi>thresh</mi><mo>+</mo><mn>1</mn> </mrow> <mn>255</mn></munderover><mi>hist</mi><mrow> <mo>(</mo> <mi>i</mi> <mo>)</mo></mrow> </mrow></mfrac><mo>-</mo><mo>-</mo><mo>-</mo><mrow> <mo>(</mo> <mn>8</mn> <mo>)</mo></mrow> </mrow>利用形态学滤波对灰度差值图像进行边缘检测，作为灰度特征提取结果；(c)以边缘变化检测结果为基准作边缘跟踪，对于无法闭合的边缘的端点，在灰度变化结果的边缘图中作邻域搜索，采用“或”逻辑对边缘变化检测结果进行补充，得到闭合的变化检测边缘及变化区域。FSA00000284735400016.tif,FSA00000284735400017.tif,FSA00000284735400018.tif,FSA00000284735400019.tif
全文摘要
本发明公开了一种基于边缘和灰度的遥感图像变化检测方法，用于解决现有的遥感图像变化检测方法精度低的技术问题。技术方案是利用基于双边滤波的Canny算法进行多时相图像边缘特征提取，然后对灰度差值图像进行OSTU阈值分割和边缘提取，获得灰度特征。再将所提取的边缘和灰度特征进行综合，检测遥感图像的变化区域。由于充分利用了图像中的线性特性的同时，利用灰度差值图像弥补了因配准误差造成的断线，遥感图像变化检测方法的精度由背景技术的87.75％提高到90.32％，提高了检测准确性。
文档编号G06T7/00GK101968885SQ201010292888
公开日2011年2月9日 申请日期2010年9月25日 优先权日2010年9月25日
发明者仝小敏, 张艳宁, 朱宇, 林增刚, 袁琪, 郗润平 申请人:西北工业大学