《测绘学报》

构建与学术的桥梁 拉近与权威的距离

虢英杰, 朱兰艳

昆明理工大学国土资源工程学院, 云南 昆明 650093

收稿日期：2016-11-16； 修回日期：2017-05-19

基金项目：国家自然科学基金（41261043）

第一作者简介：虢英杰(1993—), 男, 硕士, 研究方向为现从事定量遥感、变化信息识别和提取。E-mail:884833945@qq.com

通信作者：朱兰艳, E-mail：KMOY2222@sina.com

摘要：目前，大多数高分辨率遥感影像变化检测方法的效率不高，精度也难以达到实际应用要求，本文引入频域分析的方法进行变化检测预处理，达到提高变化信息检测精度和效率的目标。首先，对影像频域能量的分布特征进行分析，并以此为基础，利用累加距离匹配函数，完成试验影像的合理分块；其次，通过变化性纹理规则度的定义和计算，实现对同能量下变化信息的精确描述，缩小变化区域搜索范围，有效抑制背景信息对目标信息的干扰。最后，通过3种常用的变化检测方法的验证，发现变化检测预处理方法可以使检测算法在现有精度和效率的基础上有进一步提高。

A Preprocess for Change Detection Method from High Resolution Remotely Sensed Imagery in Frequency Domain

GUO Yingjie, ZHU Lanyan

Abstract: At present, the efficiency of change detection from high resolution remotely sensed imagery is generally not very high, and accuracy couldn't meet the practical requirements. A method which introduced frequency domain analysis is put forward in this paper. Firstly, it studies the distribution of image frequency energy. Using it as a basis, the image division is accomplished by summed-up distance matching function. Secondly, the definition and computation of texture regularity of change information(TRCI) clearly and accurately expresses the change information according to equal energy rule of image blocks, greatly narrowing the search area. Thus, the background information could be greatly suppressed so as to improve the efficiency and precision of detection. Finally, through measuring and comparing of the accuracy and efficiency of three different detection algorithms before and after introducing the preprocess, we find that it could provides more accurate results on the existed basis.

Key words: preprocess for change detection frequency domain analysis TRCI optimal image blocking

随着社会与科技的发展，人类开发资源与改造自然的能力不断增强，自然界的变化和人类的各种活动每天都在改变着地表景观[1]。人口的快速增长及城市化的发展，更是加快了这种变化的速度[2]。近年来，随着高分辨率遥感影像的普及化使用和相关研究的不断深入，各种变化检测方法逐渐被发掘和提出，但由于在大范围影像中，不同尺度的地物光谱相互重叠，降低了光谱域的统计可分性，使同一地物在不同时相的影像中呈现光谱、空间特征的差异，严重影响变化检测的效果，从而限制了有较高自动化水平的高分辨率遥感影像变化检测方法在实际工程中的应用[3-5]。有学者指出：通常变化区域所占的比例相对较小，变化检测需要对多时相影像分别进行地物表述，牺牲较大运算量于非变化区域，而获取少量变化信息，不仅效率低下，代价高昂，而且容易受到背景信息的干扰[6]；而这也正是许多变化检测方法之所以难以取得令人满意的精度和效率的原因所在[7]。

根据以上分析，本文尝试从影像地物的视觉特征出发，在保证地表景观完整性的前提下，确定影像中主要地物的纹理周期性和方向性，实现影像分块。以此为基础，结合子影像块的变化性纹理规则度，完成变化信息所在影像块的衡量和标识，从而缩小变化区域的搜索范围，减少非变化信息在空域和频域中对变化检测的影响，进一步达到提高变化检测精度以及工作效率的目的。

1 研究方法1.1 研究思路

本文以高空间分辨率卫星图像为试验数据，经过图像抑噪等预处理后，首先对图像进行二维离散傅里叶变换，通过简单的图像变换和数学推导进行线状特征的频谱分析，证明线状特征和频域变换后的谱线之间的方向关系；然后根据频域能量的径向分布特征和SDMF所表征的周期性，得到图像分块的分块参数；最后，结合论文所提出的变化性纹理规则度，实现有效变化信息的精确表示。完成确定变化信息所在区域的任务，为检测算法在现有的检测精度和效率的基础上做进一步提高提供了保障，论文的大致流程如图 1所示。

图 1 引入频域分析的变化检测预处理过程Fig. 1 The preprocess of change detection based on frequency spectrum analysis

图选项

1.2 影像频域能量的信息表征

为了定量分析遥感影像的空间频率及其对应能量值所反映的地物纹理、细部结构和边缘特征，通过二维离散傅里叶变换对影像进行空域向频域的转换[8]。设影像的大小为M×N，其二维空域表达为f(x,y)，对应的频域表达为F(u,v)，则二维离散傅里叶变换的数学定义式为

(1)

根据傅里叶变化位移定理，设影像的变量x,y在二维空间内存在位移x-x0和y-y0，那么二维离散傅里叶变化满足以下关系

(2)

从式(2) 可以看出，时(空)域信号的平移并不能改变其频率和幅度成分，发生改变的仅仅是相位谱部分。由于空间域内单个条带图像必定带有直流电平的属性[9]，故其所对应的谱线必然会通过频域平面的中心点，并集中反映在垂直于原特征线方向的谱线上，这一现象称之为傅氏变换幅度谱的自配准性质[10]。

为了进一步探讨影像中地物特征与能量谱线之间的关系，需要将傅里叶频谱能量分析和统计方法结合起来考虑。本文采用文献[11]中所提出的环状采样和楔状采样的方法，假设(r,θ)为(u,v)平面上的极坐标，以采样间隔为θi≤arctan(v/u)＜θi+1，在以频谱图中心为原点出发的扇面上对r进行求和，其关系为

(3)

则称Er(θi)为影像频谱能量的角向分布特征(楔特征)，如图 2(a)所示。

图 2 能量谱的环状采样和楔状采样Fig. 2 Radius sampling and angle sampling of the frequency spectrum

图选项

在高分辨率遥感影像中，地物特征线会随机分布在很多方向上[12]，但根据傅里叶变换能量谱的自配准性，相同方向上的特征线无论分布如何，最终其能量总会叠加在一起[13]。而楔状采样作为对影像频域能量分布方向性的描述方法，其角向分布特征曲线的峰值将也必然会出现在与其能量谱的谱线相垂直的方向上。

统计分析两期影像频域能量的角向分布后，即可得到两期试验影像的角向分布特征曲线[14]，从而表征研究区域主要地物的排列方向，进一步可确定影像分块的旋转参数α。

若以采样间隔为ri2≤u2+v2＜ri+12，在以频谱图中心为圆心的一个环面上对θ求和，其关系为

(4)

则称Eθ(ri)为影像频谱能量的径向分布特征(环特征)，如图 2(b)所示。

应该注意的是，在一幅遥感影像中往往包含有各种周期性成分、非周期成分以及噪声等信息[15]。在径向采样中，除了直流分量外，如果有某一波峰特别突出，往往说明原图像中包含以该波峰频率为主频率的周期成分存在，而且峰值的大小表征了周期性成分的强弱程度[16-18]。

1.3 影像分块参数的确定

大多数变化检测方法需要对地表景观进行完整的描述，这就要求在对影像进行分块的同时，尽量保证影像中主要地物的完整性，不会因影像分块导致地物破碎，从而避免影响地物提取精度以及变化检测精度[19]。传统影像分块方式及其效果见图 3。

图 3 传统影像分块方式及分块效果Fig. 3 The tradition way of image division and effectiveness

图选项

在进行影像分块的时候，为避免分块线穿过主要地物，应根据影像的频谱能量角向分布曲线的峰值θ，对影像进行顺时针旋转，使影像中主要地物尽量按照经纬方向进行排列，以期在合适的步长下完成对地物的完整分块。

为了使图像既有足够的纹理基元组合，又尽可能形成单一纹理，应针对不同影像中主要地物的纹理基元的周期性，设计影像分块步长。论文选取累加距离匹配函数(summed-up distance matching function, SDMF)对图像纹理基元的周期性进行描述[20]，通过分析影像中纹理基元的周期性，完成分块步长的选取。

文献[20]所提出的累加距离匹配函数，是在牺牲少量计算精度的代价下，提高传统单窗口灰度共生矩阵(cray level co-occurrence matrix, GLCM)对纹理基元周期性的执行效率，其旨在计算影像不同方向间不同步长的亮度差值的累积量，从而判定纹理周期。行方向和列方向上的SDMF分别可以表达为

(5)

(6)

式中，M、N分别表示图像的尺寸大小；r和c分别表示行数和列数；δ表示像素间的空间距离。通常，可以将δ的范围取作0~N-1，以获取步长为1~N下的SDMF值，从而使SDMF的纹理周期计算具有全局性，进一步确定最佳分块步长。

为了规范化探讨SDMF在具备规则性纹理的影像中的特性，通过10个大小为43×344像元的矩形构建人工规则纹理图像(见图 4(a))，其中，设黑色亮度为0，白色亮度为1。对规则纹理图像的SDMF做进一步的分析后发现，SDMF会产生大小、尺寸不一的波峰(见图 4(b))，借助SDMF的一阶、二阶导数，计算各次波峰的峰值和谷值，通过设定阈值，判断各次峰值和谷值的显著性，获取各显著峰的尺寸，进而得到纹理图像的周期。为了使显著性判定阈值具有良好的自适应性，经过规则纹理图像实验验证，论文采用文献[21]中的观点，使用公式(7) 计算显著性判定阈值

图 4 人工规则纹理图像及其SDMF分析图Fig. 4 The artificial image and statistical analysis of SDMF in the grid row

图选项

(7)

式中，SDMF的第一显著峰值是全局最优的平均稳定周期[22]，第一显著峰的尺寸即可表征纹理基元的尺寸。但对于规则性纹理图像，纹理周期是指基元特征重复出现的尺寸大小[23]。因此，SDMF的所有显著峰均可以表示纹理图像的周期，但对于步长的选取而言，需要综合考虑各种纹理基元的多样性，又要照顾到地物的相对独立性[24]。进而，论文规定步长的选取公式为

(8)

式中，S为影像分块步长；μ为单位步长，即SDMF所反映的纹理周期；m为步长倍数。

应该格外注意的是，经过反复试验，在研究区域中变化信息对整体信息的贡献率低于52%的情况下，由两期SDMF所确定的步长大小基本可以保证同一，此时将其最或是值作为单位步长即可满足分析需求；而当影像中地物发生剧烈变化，导致变化信息占比超过52%时，将难以保证两期影像的SDMF所反映出的地物周期性一致。此时，应将两期影像中较大的地物周期作为单位步长，并进行适当扩大步长倍数m，这样虽然了牺牲了少量运算效率，但总体上可以保证变化信息不会因影像分块而丢失，并可以较好地反映地物周期性。

1.4 有效变化信息所在影像块的确定

纹理感知的精神心理学认为：人类识别纹理特征的能力主要来自对纹理基元重复模式的判断，其中3个重要的特征是结构元素的规则度、重复性和方向性[25]。当研究区域地物发生变化时，必然会导致其纹理信息和边缘特征的变化，从而带来纹理规则度的改变，借此可以衡量图像变化的程度，进而缩小变化区域的搜索范围。

首先，利用纹理图像频谱的径向特征曲线E(r)统计其能量分布情况，对径向特征曲线进行归一化处理；其次，由其下包络曲线B(r)对径向特征曲线进行拟合，计算E(r)和B(r)的能量比值e；再次，将B(r)中各处的点值乘以e，得到同能量下包络曲线T(r)；最后，求取E(r)对T(r)的方差σEi，即确定图像纹理规则度[26]。

通过两期影像的纹理规则度，计算变化性纹理规则度ΔσE为

(9)

该算法可以确保下包络线和径向特征曲线所映射的幅度谱具有相同的能量，而且具有单调递减的特性，从而有效减少多期影像间因光谱能量、辐照度等条件的不同对变化检测所产生的干扰。

根据以上指标的定义过程，可以得到结论：图像中结构性纹理所占的比例越大，图像的纹理规则度就越高；图像中包含的变化信息越大，变化性纹理规则度则越高。这也就是说，变化性纹理规则度ΔσE可以反映影像块产生变化的剧烈程度，通过设定对变化信息分割阈值η，即可实现对变化信息所在影像块的准确判断。

2 试验和结果分析2.1 试验

为了验证本文所提出方法的有效性，选择2012年5月17日和2015年8月11日获取的云南某地区QuickBird卫星遥感影像作为研究对象。QuickBird影像包含了4个多光谱波段和一个全色波段，其中全色波段波谱范围为0.45~0.90 μm，全色波段空间分辨率0.61 m。截取全色波段中2800×2800像素大小的图像作为研究数据。

首先，进行图像预处理。论文采用文献[27]中所提出的适合于去除高斯噪声和椒盐噪声的偏微分方程模型进行图像噪声抑制处理。处理后发现，遥感数据亮度响应不明显，为使图像更好地呈现地物的细节结构特征，将试验数据进行直方图拉伸，尽量使拉伸后的图像直方图分布于0~255全部灰度级，以增加对比度，有利于地物类型的识别。

然后，计算影像分块的两个参数：旋转参数α和分块步长s。通过对两期影像频谱能量的角向分布特征进行统计分析，得其角向分布特征曲线，见图 5和图 6。

图 5 2012年试验影像的角向分布特征曲线Fig. 5 The image distribution of wedges feature curve in 2012

图选项

图 6 2015年试验影像的角向分布特征曲线Fig. 6 The image distribution of wedges feature curve in 2015

图选项

从影像角向分布特征曲线分析，虽然两期影像的频域能量分布形式略有不同，但可以清晰地表征出2012年和2015年试验影像块的频域能量明显集中于47°和137°，其余的方向并未出现明显的峰值，表现为一条较为平稳的曲线。这表明原影像中包含大量位于该能量峰值方向上的边缘或直线信息，而主要地物的优势边缘是按照47°的方向进行排列，而其余地物特征未表现出明显方向性。

为了使主要地物按照经纬方向排列，尽量避免由于分块不当而带来的景观破碎，取试验影像角向特征分布曲线峰值方向的余角α作为旋转参数，将影像顺时针旋转47°，得到适合分块的影像分布方式，见图 7。

图 7 旋转处理后的两期试验影像Fig. 7 Two experiment images of rotation in 2012 and 2015

图选项

对2012年试验影像的SDMF计算值进行统计，计算其SDMF的一阶导数f′(x)、二阶导数f″(x)，分析各次显著峰的波峰与波谷的位置，从而准确反映出影像的纹理周期单位，该SDMF统计图及其分析见图 8。

图 8 2012年试验影像的SDMF统计分析结果Fig. 8 The SDMF statistical results of experiment image in 2012

图选项

通过对一阶、二阶导数的判断，得到2012年试验影像SDMF计算结果中所有峰值和谷值。经计算，确定本实验的阈值为4.606×107，统计各次显著波峰值、谷值及其所对应步长，见表 1。

表 1 2012年SDMF的各次显著波峰值、谷值统计表Tab. 1 The peak and valley value of SDMF and the associated significant wave statistical list in 2012

表选项

为了进一步考察在变化信息所占比例相对较小的两期影像中SDMF的分布情况，统计2015年影像SDMF计算结果中各次显著波峰值、谷值及其所对应步长，见表 2。

表 2 2015年SDMF的各次显著波峰值、谷值统计表Tab. 2 The peak and valley value of SDMF and the associated significant wave statistical list in 2015

表选项

结合表 1和表 2可以看出，两期影像中各次特征点所对应的步长值仅出现微小扰动，但并没有固定的偏离方向和大小，呈现出偶然误差的特性。通过计算整条显著波内各周期步长大小的最或是值，发现两期遥感影像中SDMF所反映的单位步长分别为：219.655、219.725。这进一步说明在变化区域占比相对较小的情况下，两期影像的SDMF所反映出的周期性规律基本具备同一性。

为了保证分块后的子影像中包含尽量多的地物种类，结合试验验证，以两倍显著峰尺寸作为步长，即以440像元为步长，将影像分成81幅子图像，其结果可基本保证影像块中地物的完整性。

最后，在完成影像分块后，分别统计各子图像块的频谱能量径向分布，计算其同能量下包络线，获取各影像块的纹理规则度，最终得到对应影像块的变化性纹理规则度，并按其行列分布统计为表 3。

表 3 试验图像的变化性纹理规则度Tab. 3 The value of TRCI about exemplificative images

表选项

通过对各影像块的变化性纹理规则度进行统计分析后发现，除影像信息空白区的变化性纹理规则度为0外，其余54幅子影像块均存在大小不一的变化信息；进一步的，通过人工判断表中不同纹理规则度在提取变化信息能力上的强弱，确定变化信息分割阈值η=0.03，通过对变化信息的描述，有效变化信息所在的影像块已被准确标识(见表中加粗部分)，以子图像块(6, 2) 为例，验证变化性纹理规则度对变化信息描述的效果，见图 9和图 10。

图 9 2012年试验影像块的归一化径向分布曲线及同能量下包络曲线Fig. 9 The normalized radial distribution curve and co-energy lower enveloping curve in image block in 2012

图选项

图 10 2015年试验影像块的归一化径向分布曲线及同能量下包络曲线Fig. 10 The normalized radial distribution curve and co-energy lower enveloping curve in image block in 2015

图选项

结合两期验证子影像块发现，在2012年试验影像块的归一化径向分布曲线中，径向分布曲线和下包络曲线的频谱能量差之和为1.696，纹理规则度σE1=0.014；2015年试验影像块的径向分布曲线中存在明显波峰，其径向分布曲线和下包络曲线的频谱能量差之和为2.432，纹理规则度达到σE2=0.062。

两期试验影像块的变化性纹理规则度ΔσE=0.048>0.03。反映出两期影像块的结构性纹理发生较剧烈的变化，通过对试验影像块的目视解译，可以发现该区域确实产生大量的地物变化，验证了此方法的可靠性和有效性。

2.2 评价与分析

目前变化检测的主流方法大致可以分为面向像元级、面向对象级、面向目标级3类，试验将采用3类检测手段中较为常用的方法，通过对比使用预处理方法前后的检测精度及效率的变化情况，完成对变化预处理方法的评价。

对于像元级变化检测方法，采用影像比值运算法，是指将两期遥感影像中对应的像元值做比值运算，获取差异影像，通过选取适当的阈值，实现对变化区域的提取。

对于对象级变化检测方法，采用MeanShift-CVA法，是指通过MeanShift图像分割算法进行地理对象的提取，然后通过变化矢量分析法(change vector analysis, CVA)计算两个不同时期的影像差异，最后通过阈值判断，完成变化信息的提取。

对于目标级变化检测方法，以建筑物为例，通过傅里叶变换实现影像从空域向频域的变换，统计其频谱能量分布，获取建筑物的纹理区域和边缘特征，构建相适应的滤波器，完成建筑物的二值化提取，最后通过对二值化信息进行像素级的变化检测即实现变化信息的提取。

通过以上3种代表方法在同一硬件环境、同一运算平台下，对试验影像进行变化检测，统计其执行变化检测预处理前后各算法的检测精度和效率(计入变化检测预处理时间)，见表 4。

表 4 3种变化信息检测方法执行变化检测预处理前后检测精度和效率变化表Tab. 4 The comparisons of accuracy and efficiency before and after the preparing preprocess by three different change detection methods

表选项

针对统计结果，从宏观层面看，3个级别的变化检测方法在进行预处理后，检测精度和运行效率均有所提高，这充分地反映了本文所提出的变化检测预处理方法对提高检测精度及其执行效率的有效性。

从微观层面看，检测预处理对3种检测级别起到的作用有所差异。其中，像素级比值法精度较未进行预处理前仅提高3.39%，运行时间减少5.79 s，运算效率提高10.80%，是3种方法中收效最差的，其根本原因在于像素级变化检测方法对高分辨率遥感影像的不适用性，导致检测预处理对其效果较差。对象级MeanShift-CVA法检测精度提高9.54%，节省运行时间45.38 s，运算效率提高53.57%；在此基础上，作为对检测结果纹理信息同质性的描述，其紧质度也有一定的提高。频域下目标级变化检测法精度提高12.50%，节省运行时间69.69 s，运算效率提高54.53%，在3种方法中，是收效最大的，这是因为在目标级变化检测中，需要对两期影像分别进行目标地物的提取，其往往会受到背景地物在光谱和频域上的影响，通过对影像分块，很大程度的抑制了目标地物的提取误差，从而提高变化检测精度。

进一步的，从算法性质上来讲，由于全自动方法在很多实际应用中会产生精度不足、算法不够健全和稳定等问题，尚无法适用于多种复杂情况[28]；故本文提出了一种人工干预的半自动化处理方法，借助复杂情况下的人工判断，对检测过程加以修正，最终达到提高检测精度的要求，而恰当的干预次数和干预位置就成为半自动化算法实现高精度、高效率的关键所在。论文在选取步长倍数和变化信息分割阈值两处操作上进行人工干预。经过试验验证，仅使人工干预次数在原算法的基础上平均增加34.12%，这意味着，所进行的人工干预不但可以对复杂情况下的参数进行弹性设置，而且不会对算法原有的自动化程度造成过多的影响，从而侧面验证了本实验的科学性。

应该注意的是，目标级变化检测的精度很大程度上是由目标地物的提取精度决定，变化检测预处理对提高目标级变化检测精度的有效性，进一步从侧面反映出，论文中的影像分块方法并没有引起大规模的地物破碎，而造成变化信息丢失、检测精度下降，再次证明了影像分块方法的科学性和合理性。

3 结语

本文以频域能量守恒定律和傅立叶变换幅度谱的自配准性质为依据，通过合理的影像分块，结合对应子影像块的变化性纹理规则度，确定有效变化信息的范围；研究表明，高分辨率遥感影像变化检测的频域分析预处理方法可以帮助原检测算法实现更高层次的检测目标，对于减少背景信息的干扰、缩小变化区域搜索范围具有一定的参考价值。但应该注意到的是，检测精度及效率所能够提高的程度主要由原检测算法的本质及其性质决定，如何弱化其对预处理方法的较大影响，使之更具普适性，是一个有待于进一步研究和解决的问题。

【引文格式】虢英杰，朱兰艳。 高分辨率遥感影像变化检测的频域分析预处理方法[J]. 测绘学报，2017，46(6)：743-752. DOI: 10.11947/j.AGCS.2017.20160570

更多精彩内容：