基于内容的图像检索技术为解决海量遥感影像信息提取和共享难题提供了新的契机。如何克服遥感影像数据在获取过程中由于受角度和分辨率等影响而产生的旋转及尺度差异，是遥感影像检索区别于普通图像检索的特别之处，也是图像检索领域备受关注的前沿问题，而构造旋转与尺度不变特征是解决该问题的关键^[1]。

一类方法是采用全局或局部特征进行检索，如文献[2]提出了一种基于Z搜索的彩色图像检索方法，对每个颜色通道构建Z搜索的共生矩阵；文献[3]提出了尺度不变特征变换(scale invariant feature transform，SIFT)特征，该特征是一种局部描述子，具有平移旋转和尺度不变等性质。另一类方法是利用图像变换进行处理，如文献[4-5]将极坐标（Log-Polar）或Radon变换与小波变换相结合提取不变特征进行旋转纹理图像检索；文献[6]利用Gabor和Gauss滤波系数的归一化能量矢量进行旋转纹理识别；文献[7]使用轮廓波（Contourlet）子带能量分布信息构造了旋转不变的特征；文献[1]结合了Log-Polar变换、Radon变换、非下采样轮廓波变换(nonsubsampled contourlet transform, NSCT)和小波变换，构造了一种多尺度、多方向的纹理旋转不变特征矢量。但上述研究成果存在两个问题：（1）将彩色图像的3个色彩通道单独进行处理，或者将彩色图像灰度化，导致丢失了3个色彩通道之间的关系信息或是直接丢失了颜色信息；（2）仅考虑了如何保持旋转和尺度不变性，没有充分利用数据本身所特有的几何特征，从而影响了遥感影像的检索精度。

文献[8]提出了四元数的概念，彩色图像可以表示成纯四元数的形式来进行处理，既实现了对彩色图像进行整体的处理，又保留了3个颜色通道之间的关联信息，能有效利用颜色信息。文献[9]提出了一种基于四元数的彩色遥感影像水域提取算法，指出了四元数在彩色影像上应用的优势。文献[10]提出了一种四元数旋转不变彩色图像纹理分类方法，但其对尺度变换具有敏感性。本文提出了一种基于四元数变换的多特征遥感影像检索方法，将正交傅里叶梅林矩扩展至四元数域，构造旋转与尺度不变量；使用矩形加楔形的方法对不变量频谱进行特征提取，构造包含纹理信息的旋转和尺度不变的特征量；对图像边缘进行检测，得到边缘色彩图像，并将边缘上的色彩映射到25种视觉颜色中，得到边缘色彩特征。该方法综合利用了色彩、纹理和边缘特征进行检索，对于存在旋转和尺度变换关系的遥感图像检索有较高的精度，而且旋转和尺度不变量具有较好的鲁棒性。

1 四元数及彩色图像的四元数表示

四元数是一种简单的超复数，由实数部分和3个虚数部分i、j、k组成。一般可以写成以下的形式：

$q = a + b{\rm{i}} + c{\rm{j}} + d{\rm{k}}$

(1)

其中，a、b、c和d表示实数；i、j、k满足以下性质：

${\rm{i}} \times {\rm{i}} = {\rm{j}} \times {\rm{j}} = {\rm{k}} \times {\rm{k}} = - 1$

(2)

$ \left\{ {\begin{array}{*{20}{c}} {{\rm{i}} \times {\rm{j}} = - {\rm{j}} \times {\rm{i}} = {\rm{k}}}\\ {{\rm{j}} \times {\rm{k}} = - {\rm{k}} \times {\rm{j}} = {\rm{i}}}\\ {{\rm{k}} \times {\rm{i}} = - {\rm{i}} \times {\rm{k}} = {\rm{j}}} \end{array}} \right. $

(3)

四元数的加法和减法与一般的复数相似，乘法满足分配率，但是对于交换率却不满足，即q₁×q₂和q₂×q₁不相等。

根据四元数的性质，可以将一幅彩色图像表示成纯四元数的形式，彩色图像的每个像素都可以视为一个实部为零的纯四元数:

$ f\left( {x,y} \right) = {f_R}\left( {x,y} \right){\rm{i}} + {f_G}\left( {x,y} \right){\rm{j}} + {f_B}\left( {x,y} \right){\rm{k}} $

(4)

式中，f (x，y)表示彩色图像在(x，y)位置的像素值；f_R (x，y)、f_G (x，y)、f_B (_x，y)分别表示该像素的红、绿、蓝分量。由式（4）可知，彩色图像可以表示为一个纯四元数矩阵。

2 基于四元数正交傅里叶-梅林矩的遥感影像检索方法

正交傅里叶-梅林矩(orthogonal Fourier-Mellin moments, OFMMs)是一种正交不变矩，对图像的形状特征有着很好的表达能力，还具有对图像噪声不敏感的特性。本文将正交傅里叶梅林矩扩展至四元数域，提出了一种基于四元数正交傅里叶-梅林矩(quaternion orthogonal Fourier -Mellin moments，QOFMMs)的遥感影像检索方法。

2.1 四元数正交傅里叶梅林矩

QOFMMs继承了OFMMs模值具有旋转不变性和良好抗噪性的特性，此外，还能实现对彩色图像的特征描述，保证颜色通道之间的相关信息不会丢失。由于四元数不满足乘法交换律，所以四元数正交傅里叶-梅林变换分为左边变换、右边变换和双边变换，因左边变换更加方便计算，本文对四元数左边正交傅里叶-梅林矩极坐标下的定义G_pq(f)为：

$ {{G}_{pq}}\left( f \right)=\frac{1}{2\text{ }\!\!\pi\!\!\text{ }{{a}_{p}}}\overset{2\text{ }\!\!\pi\!\!\text{ }}{\mathop{\underset{0}{\mathop{\mathop{\int }^{}}}\,}}\,\overset{1}{\mathop{\underset{0}{\mathop{\mathop{\int }^{}}}\,}}\,{{Q}_{p}}\left( r \right){{\text{e}}^{-\mu q\theta }}f\left( r,\theta \right)r\text{d}r\text{d}\theta $

(5)

式中，a_p为归一化常量；Q_p (r)是关于r的p次多项式，在单位圆上正交；p、q为非负整数；μ为单位纯四元数；f (r，θ)为极坐标形式下彩色图像的四元数表示；r表示极坐标系中的极径；θ为极角。f (r，θ)定义为：

$ f\left( r,\theta \right)={{f}_{R}}\left( r,\theta \right)\text{i}+{{f}_{G}}\left( r,\theta \right)\text{j}+{{f}_{B}}\left( r,\theta \right)\text{k} $

(6)

式中，f_R (r，θ)、f_G (r，θ)、f_B (r，θ)分别为对红、绿、蓝3个色彩分量进行极坐标变换后得到的极坐标分量；i、j、k为四元数虚部。a_p、Q_p (r)、μ的定义为：

$ {{a}_{p}}=\frac{1}{2\left( p+1 \right)} $

(7)

$ {Q_p}\left( r \right) = \sum\limits_{s = 0}^p {{{\left( { - 1} \right)}^{p + s}}\frac{{\left( {p + s + 1} \right)!}}{{\left( {p - s} \right)!s!\left( {s + 1} \right)!}}} {\mkern 1mu} $

(8)

$ \mu =-\frac{\text{i}+\text{j}+\text{k}}{\sqrt[{}]{3}} $

(9)

式中，s ∈ [ 0，p ]。

设旋转角度为α，容易证明QOFMMs的模值具有旋转不变性，即：

$ \left| {{G}_{pq}}\left( f\left( r,\theta +\alpha \right) \right) \right|=\left| {{G}_{pq}}\left( f\left( r,\theta \right) \right) \right| $

(10)

针对正交多项式Q_p (r)计算的复杂性，文献[11]提出了一种基于Q_p (r)的多项式之间的递推公式，有效提高了正交傅里叶-梅林矩的计算效率。基于Q_p (r)的递推式为：

$ {{Q}_{p}}\left( r \right)=\left( {{L}_{1}}r+{{L}_{2}} \right){{Q}_{p-1}}\left( r \right)+{{L}_{3}}{{Q}_{p-2}}\left( r \right) $

(11)

其中，Q₀ (r) = 1，Q₁ (r) = 3r - 2

2.2 利用四元数正交傅里叶-梅林矩构造彩色图像不变量

虽然QOFMMs的模值不具有缩放不变性，但可以通过代数表达式来构造彩色图像的旋转和缩放不变量。

假设f (r，θ)为极坐标下表示的一幅图像，经过缩放λ倍、旋转α角度变换后图像变为f (λr，θ + α)，使用f ′表示变换后的图像，则旋转和尺度不变量的构造过程如下

$ \begin{align} & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{G}_{pq}}\left( {{f}'} \right)= \\ & \ \ \ \frac{1}{2\text{ }\!\!\pi\!\!\text{ }{{a}_{p}}}\overset{2\text{ }\!\!\pi\!\!\text{ }}{\mathop{\underset{0}{\mathop{\mathop{\int }^{}}}\,}}\,\overset{1}{\mathop{\underset{0}{\mathop{\mathop{\int }^{}}}\,}}\,{{Q}_{p}}\left( r \right){{\text{e}}^{-\mu q\theta }}f\left( \lambda r,\theta +\alpha \right)r\text{d}r\text{d}\theta = \\ & \ \ \ \frac{1}{2\text{ }\!\!\pi\!\!\text{ }{{a}_{p}}}\overset{2\text{ }\!\!\pi\!\!\text{ }}{\mathop{\underset{0}{\mathop{\mathop{\int }^{}}}\,}}\,\overset{1}{\mathop{\underset{0}{\mathop{\mathop{\int }^{}}}\,}}\,{{Q}_{p}}\left( \frac{r}{\lambda } \right){{\text{e}}^{-\mu q\left( \theta -\alpha \right)}}f\left( r,\theta \right)\frac{r}{\lambda }\text{d}r\text{d}\theta = \\ & \frac{1}{2\text{ }\!\!\pi\!\!\text{ }{{a}_{p}}}{{\text{e}}^{-\mu q\alpha }}{{\lambda }^{-2}}\overset{2\text{ }\!\!\pi\!\!\text{ }}{\mathop{\underset{0}{\mathop{\mathop{\int }^{}}}\,}}\,\overset{1}{\mathop{\underset{0}{\mathop{\mathop{\int }^{}}}\,}}\,{{Q}_{p}}\left( \frac{r}{\lambda } \right){{\text{e}}^{-\mu q\theta }}f\left( r,\theta \right)r\text{d}r\text{d}\theta = \\ & \ \ \ \ \ \ \ \frac{1}{2\text{ }\!\!\pi\!\!\text{ }{{a}_{p}}}{{\text{e}}^{-\mu q\alpha }}{{\lambda }^{-2}}\overset{p}{\mathop{\underset{m=0}{\mathop{\mathop{\sum }^{}}}\,}}\,\overset{p}{\mathop{\underset{k=m}{\mathop{\mathop{\sum }^{}}}\,}}\,{{\lambda }^{-k}}{{c}_{pk}}{{d}_{km}}{{G}_{mq}}= \\ & \ \ \ \ \ \ \ \ \ \ \frac{1}{2\text{ }\!\!\pi\!\!\text{ }{{a}_{p}}}{{\text{e}}^{-\mu q\alpha }}\overset{p}{\mathop{\underset{m=0}{\mathop{\mathop{\sum }^{}}}\,}}\,\overset{p}{\mathop{\underset{k=m}{\mathop{\mathop{\sum }^{}}}\,}}\,{{\lambda }^{-\left( k+2 \right)}}{{c}_{pk}}{{d}_{km}}{{G}_{mq}} \\ \end{align} $

(12)

其中，

设旋转和尺度不变量的表达式为:

$ {{M}_{pq}}=\frac{1}{2\text{ }\!\!\pi\!\!\text{ }{{a}_{p}}}{{\text{e}}^{-\mu q\alpha }}\overset{p}{\mathop{\underset{m=0}{\mathop{\mathop{\sum }^{}}}\,}}\,\overset{p}{\mathop{\underset{k=m}{\mathop{\mathop{\sum }^{}}}\,}}\,{{\lambda }^{-\left( k+2 \right)}}{{c}_{pk}}{{d}_{km}}{{G}_{mq}} $

(13)

由G₀₀(f ') = λ^-2G₀₀(f)，可得$\lambda ={{\left| {{G}_{00}}(f) \right|}^{-\frac{1}{2}}}$，同理可得α = tan^-1(G₁₁(f))

综上，可以得到QOFMMs的旋转与尺度的不变量:

$ {{\text{e}}^{-\mu q{{\theta }_{f}}}}\overset{p}{\mathop{\underset{m=0}{\mathop{\mathop{\sum }^{}}}\,}}\,\overset{p}{\mathop{\underset{n=m}{\mathop{\mathop{\sum }^{}}}\,}}\,{{\Gamma }^{-\left( n+2 \right)}}{{c}_{pn}}{{d}_{nm}}{{G}_{mq}} $

(14)

其中，m、n为非负整数；θ_f = arctan (G₁₁)；$\mathit{\Gamma }=\sqrt{\left| {{G}_{00}} \right|}$；G_mq为m阶q重复度的四元数正交傅里叶梅林矩；c_pn、d_nm分别为：

2.3 基于四元数正交傅里叶-梅林矩的图像特征提取

傅里叶变换的频谱图能在一定程度上反映图像的纹理信息，文献[12]利用长方环傅里叶周向谱能量百分比法提取图像纹理特征，结果证明了其方法比进行环形采样有着更高的精度。文献[13]基于四元数傅里叶-梅林矩的频谱提出了环形与楔形的特征提取方法，得到了旋转不变的纹理特征，在纹理分类中有着不错的表现，但其频谱仅仅具有旋转不变的特性，且对方形频谱进行环形采样时会导致信息利用不全面以及存在几何误差等问题。

本文中采用矩形和楔形相结合的方法来对四元数正交傅里叶-梅林矩的不变量M_pq（p = 1，2…N，q = 1，2…N）的频谱图进行特征提取。对于N×N的M_pq频谱图。首先将其均匀地分为n个方形环，最外层的环编号为1，最内层为n，每个环的宽度为N/（2n）。然后将每个环按照楔形分为8个梯形（三角形），一共得到8n个采样区域。如图 1所示。

将每个采样区域中的M_pq数值进行累加，然后除以所有采样区域之和，便可以得到每个采样区域的能量占比，由此共可以获得8n个特征值。对于一幅N × N的数字图像，首先计算它的N × N大小的旋转和尺度不变量，然后再通过图 1分块累加计算便可以计算出一组大小为8n的特征向量。该特征向量在一定程度上反映了图像的纹理特征，而且具有旋转和尺度的不变性。定义E_iu为第i层的第u个采样区的特征向量（1 ≤ u ≤ 8），则对应第i层的8个采样区域的计算公式为：

QOFMMs的模可以用来对图像进行边缘检测^[13]。假设彩色图像中的某个像素点为单位圆的圆心，且图像的某个边缘在这个圆中，则可以建立边缘检测的模型。如图 2所示，L表示彩色图像中的某条边缘，θ表示该边缘与Y轴的夹角，l是原点与边缘的垂直距离，h表示原背景值，h+k表示跃迁后的背景值，k为跃迁高度。将图像顺时针旋转θ角后，由于模值保持不变，可计算得出：

$ \theta =\text{ta}{{\text{n}}^{-1}}\frac{I\left( 2{{G}_{01}}+{{G}_{11}} \right)}{R\left( 2{{G}_{01}}+{{G}_{11}} \right)} $

(15)

其中，I ( )表示虚部值；R ( )表示实部值。

根据QOFMMs的矩参数，可以计算出k和l的值：

$ l=\frac{3\left( 4G{{\text{ }\!\!'\!\!\text{ }}_{10}}+G{{\text{ }\!\!'\!\!\text{ }}_{20}} \right)}{5\left( 2G{{\text{ }\!\!'\!\!\text{ }}_{01}}+G{{\text{ }\!\!'\!\!\text{ }}_{11}} \right)} $

(16)

$ k=\frac{2G{{\text{ }\!\!'\!\!\text{ }}_{01}}+G{{\text{ }\!\!'\!\!\text{ }}_{11}}}{2{{(1-{{l}^{2}})}^{\frac{1}{2}}}} $

(17)

其中，G′_pq为图像旋转θ后的p阶q次左边正交傅里叶梅林矩变换值。

使用QOFMMs的核函数在单位圆上积分构造4个7×7的矩模板，通过模板与四元数表示的彩色图像进行四元数卷积运算后，对图像进行边缘检测，计算出l、θ的值和边缘处跃迁高度k，设置合适阈值，计算得到边缘的实际位置(x'，y')为：

$ \left[ \begin{matrix} {{x}'} \\ {{y}'} \\ \end{matrix} \right]=\left[ \begin{matrix} x \\ y \\ \end{matrix} \right]+\frac{2l}{7}\left[ \begin{matrix} \text{cos}\theta \\ \text{sin}\theta \\ \end{matrix} \right] $

(18)

将检测出的边缘图像与原图像相乘，可以得到边缘色彩图像，使得图像特征不仅仅受图像颜色的影响，还受到图像中主体的轮廓信息的影响，使具有相似颜色不同类别的图像更加具有区分性。为了提取边缘色彩图像的特征，本文基于25种视觉颜色构造边缘色彩直方图来进行特征匹配。文献[14]给出了25种视觉颜色查色表，本文中使用欧氏距离来将边缘色彩图像中的颜色映射到查色表中，对于边缘色彩图像中的每一种颜色C (C _R，C_G，C_B)到查色表中第i种颜色C_i(C_iR，C_iG，C_iB)的距离D_i定义为：

$ {{D}_{i}}=\sqrt[{}]{{{({{C}_{R}}-{{C}_{{{i}_{R}}}})}^{2}}+{{({{C}_{G}}-{{C}_{{{i}_{G}}}})}^{2}}+{{({{C}_{B}}-{{C}_{{{i}_{B}}}})}^{2}}} $

(19)

根据式（18），将边缘色彩图像中的每一种颜色分配到查色表中与它距离最近的颜色中，可以得到具有25个元素的特征向量，即图像的边缘色彩特征。

2.4 相似性度量

选择频谱环的分割的层数n = 16，可以得到含有16×8个元素的特征向量E，则两幅图像的特征向量E₁和E₂的相似度距离d₁定义为：

$ {{d}_{1}}=\frac{1}{128}\sum\limits_{i=1}^{16}{\sum\limits_{j=1}^{8}{{{\left( {{E}_{{{1}_{q}}}}-{{E}_{{{2}_{q}}}} \right)}^{2}}}} $

(20)

两幅图像的边缘色彩特征向量C₁、C₂间的距离d₂定义为：

$ {{d}_{2}}=\frac{1}{25}\sum\limits_{i=1}^{25}{\frac{{{\left( {{C}_{{{1}_{i}}}}-{{C}_{{{2}_{i}}}} \right)}^{2}}}{\max \left( {{C}_{1,i}}{{C}_{{{2}_{i}}}} \right)}} $

(21)

综合以上两种特征间的距离，使用参数λ来调和特征间的距离，λ代表距离d₁的权重，两个距离间的权重应满足和为1。两幅图像间相似度距离定义为：

$ d=\lambda {{d}_{1}}+\left( 1-\lambda \right){{d}_{2}} $

(22)

距离越近，说明两幅图像越相似。

3 算法时间复杂度分析

本文提出的基于四元数变换的彩色遥感影像检索方法，计算复杂度主要考虑3个方面：四元数正交傅里叶-梅林矩计算、构造基于正交傅里叶-梅林矩的彩色图像不变量和边缘色彩特征提取。实验根据多项式之间的递推关系，使用低阶多项式的值来计算高阶多项式的值，降低算法计算复杂度。

具体分析如下：对四元数左边傅里叶-梅林矩计算原理进行分析，该部分算法时间复杂度与计算阶数p有关，且算法时间复杂度为O(p)。对QOFMMs的旋转与尺度不变量M_pq进行分析，算法时间复杂度可分为d_nm、c_pn、G_mq三部分。对于G_mq部分，算法计算复杂度为O(p)。对于d_nm和c_pn部分，其算法复杂度根据式(13)和该不变量的构造过程易得c_pn的时间复杂度为O(p)，d_nm的时间复杂度为O(n)（0≤m≤n≤p），由于c_pn、d_nm的实现方式为嵌套，所以，c_pnd_nm的算法时间复杂度为O(p²)，即该部分算法时间复杂度为O(p²)。根据边缘色彩特征提取实现原理，易得该部分算法时间复杂度为O(N²)。综上，对于有N×N个像素点的图像，基于四元数正交傅里叶梅林矩的图像特征提取总的算法时间复杂度为O(N²p²)。

4 基于四元数变换的遥感影像检索实验

本文从美国地质勘探局（United States Geological Survey，USGS）的UCMerced_LandUse的遥感影像数据库中选取了10类遥感图像^[15]，每类中含有40幅图像，将原始图像通过4次旋转、4次缩放，共得到6 400幅图像，将原始400幅图像作为待检索的图像，6 400幅图像作为数据库图像。

纹理特征和边缘特征在进行相似性度量时权值λ为经验值。图 3所示为λ取不同权值的实验对比结果，本文实验选择λ。为了验证本文方法对图像旋转和尺度的不变性，对每幅待检索图像，检索到的相关图像N（令N = 16）一定时，设此时一共检索出的图像为M，则检索的准确率P表示为：P = N/M。由于视觉词袋模型算法(bag-of-visual words mode，BoVW）^[16]在遥感影像中已得到广泛应用，其所采用的SIFT特征具有良好的旋转和尺度的不变性，同时对图像具有很强的描述能力，本文使用BoVW模型（视觉词典维度d=300）来进行对比。实验结果如表 1所示，从表 1中可以看出，使用本文提出的综合特征对于每一类图像都具有较好的旋转和尺度的不变性，而SIFT特征对于大多数的图像虽然也具有良好的不变性，但对于某些灰度变换不明显的图像则具有较差的鲁棒性。实验结果说明了本文提出的不变量在图像的旋转和尺度变换方面比SIFT特征具有更好的鲁棒性。

为进一步评价本文方法的检索效果，采用P-R曲线来评价算法的检索精度。P表示检索的准确率，R表示检索的召回率。P和R的定义为:

$ \left\{ \begin{matrix} P=N/M \\ R=N/K \\ \end{matrix} \right. $

(23)

式中，N表示检索到的相关影像；M表示所有检索到的影像；K表示图像库中所有相关影像。

使用本文方法与视觉词袋模型、文献[11]提出的基于四元数傅里叶-梅林变换的算法、文献[17]提出的边缘直方图方法以及非四元数的OFMMs不变量加边缘色彩特征的算法进行了对比。实验结果如图 4所示，从图 4（a）、4（b）中可以看出，在这两类形状轮廓较为明显的图像库中，本文提出的算法具有更高的检索精度。从图 4（c）可以看出，本文方法对于主体形状或轮廓不太明显的图像检索精度会相对下降。从图 4（d）中可以看出，本文提出的方法的平均检索精度比其他4种方法更高，而且使用四元数对彩色图像进行整体的处理比将彩色图像进行灰度化然后再进行特征提取的方法有着更高检索精度。这说明使用四元数可以保存3个颜色通道之间的关联信息，从而使检索结果具有更高的精度。

图 5给出了对一幅遥感影像进行检索的结果，列出了相似度最高的前40幅图像, 其中第一幅图像表示待检索的图像。从检索出的前16幅图像中可以看出，本文方法对于旋转和尺度变化具有良好的鲁棒性，所使用的综合特征对于一般的图像也具有较好的检索精度。

将本文方法与BoVW进行对比，本文方法对于每幅图像的平均检索时间为0.097 s，相较于BoVW的平均检索时间0.38 s有明显缩短。由表 2可知，本文方法生成检索特征的时间比BoVW短，而BoVW生成时间受所使用的聚类算法影响较大，实验采用了最为常用的K-means算法。

5 结语

本文提出了一种基于四元数正交傅里叶-梅林变换的遥感影像检索方法。针对已有工作中存在的基于四元数傅里叶-梅林矩所提取纹理特征只具有旋转不变性的问题，构造出具有旋转和尺度不变的纹理特征，同时使用QOFMMs对图像进行边缘检测，得到边缘色彩特征。结合以上几种特征对遥感图像进行检索，实验结果证明其对存在旋转和尺度缩放关系的彩色遥感影像检索具有较高的检索精度，高于对色彩通道单独进行处理的结果，且旋转和尺度不变量具有较好的鲁棒性，适用于遥感图像的检索。