对海洋磁测资料的分析、处理和解释，通常需要根据磁异常特征对不同区域的异常进行分区，再结合地质情况作进一步的处理分析.对磁异常特征进行异常分区有助于了解洋壳基底属性、揭示岩浆活动的影响范围、划分构造单元、了解海底扩张过程等.海洋条带磁异常反映了海底扩张的过程，记录了海底扩张的时代及构造信息.海盆内的条带磁异常常因断裂、海山及后期岩浆活动的改造而形成不同的磁异常纹理特征，给条带磁异常的追踪和识别带来困难.磁异常分区是依据磁异常平面图中磁异常的走向、形态、分布范围、幅值、梯度、极值等纹理特征进行分区，通常是由经验丰富的专家进行，但人工分区过程较为繁琐，对区块边界的划分时常受主观影响.近年来，随着计算机技术的发展，图像的纹理分析识别技术发展迅速，可利用图像的纹理分析识别技术对海洋磁异常资料的纹理特征进行提取、识别和分区，为海洋磁异常的分区解释提供定量可供参考的依据，从而减轻磁异常分区工作的负担.

在磁异常资料解释中，通常基于异常的极值点和零点信息，用边界检测和边缘强化技术来推断确定异常的平面边界，如利用垂向导数、水平导数、解析信号振幅或者它们的组合来实现异常边界的提取和增强.Wijin等^[1]提出Theta图法，以提高对深部弱异常的增强识别能力.Ferreira等^[2]利用水平导数求tilt梯度的方法，增强了异常边界的横向分辨率.严良俊等^[3]将小子域滤波方法用于磁异常边界的识别.Cooper等^[4]提出归一化标准差方法，利用极大值位置确定异常边界.马国庆等^[5]利用不同阶水平导数线性组合算法进行位场边界的识别.王彦国等^[6]基于迭代差分的视导数进行边界识别，提高了边界识别的抗干扰能力.王万银等^[7]对常用的边界识别方法建立模型进行对比分析，对不同方法的效果进行了对比总结.这些方法均是通过提取识别孤立场源体的边界来对异常进行边界识别.而对于较大区域的海洋磁异常资料的处理解释，通常需要对具有相同形态、走向、幅值等特征的区域进行分块处理，对可能具有相同起源和构造特征的异常区域划分为一个区块，以便作进一步的处理解释.在区域磁异常的分区中，要求解释人员具有丰富的经验，而且因个人标准不同，可能造成不同的分区结果.因此，许多学者对异常的分区开展了定量化的自动识别技术，如姜效典^[8]探索应用人工神经网络方法来对南海磁异常进行分区；Zhang等^[9]利用改进的Randon变换和梯度实现对重磁图像线性特征的检测.陈永良等^[10]将图像方法用于重磁极值位置的自动提取，也取得较好的效果.

本文采用Gabor滤波器和灰度共生矩阵相结合的图像纹理特征分析方法，提取海洋磁异常的纹理特征信息，然后通过聚类实现不同纹理特征磁异常场的分区，给磁异常的分区解释提供定量参考.

1 方法原理

Gabor滤波器和灰度共生矩阵在图像纹理特征信息提取分析中具有较好的效果，Gabor滤波器通过不同尺度和方向的Gabor小波提取磁异常局部到区域不同尺度的时间域频率域信息，对具有方向性的线性纹理信息的提取灵敏度较高，而灰度共生矩阵可以提取异常的空间依赖性，刻画异常的空间展布和结构特征.因此，本文采用Gabor滤波器和灰度共生矩阵相结合的图像纹理分析方法提取海洋磁异常的纹理特征信息，然后利用OKMS聚类算法对提取的纹理特征信息进行聚类分区，从而实现对不同纹理特征磁异常的分区.

1.1 Gabor滤波器

磁异常图像可以看成是一个二维离散序列，Gabor滤波器可由二维Gabor函数构成.二维Gabor函数是复正弦调制高斯函数，其表达式为

$\begin{align} &g\left( x,\text{ }y \right)=\left( \frac{1}{2\pi {{\sigma }_{\text{x}}}{{\sigma }_{y}}} \right)\text{exp}\left[ -\frac{1}{2}\left[ \frac{{{x}^{2}}}{{{\sigma }_{x}}}+ \right. \right. \\ &\quad \quad \quad \left. \left. \frac{{{y}^{2}}}{{{\sigma }_{y}}} \right]+2\text{ }\!\!\pi\!\!\text{ j}\omega x \right] \\ \end{align}$

(1)

式中：σ_x, σ_y分别为g(x, y)沿x和y方向的标准方差；ω=2π/λ为高斯函数的复调制频率，λ为波长.

以g(x, y)为母小波函数，进行适当的尺度变换和旋转变换就可以得到一组不同频率、不同方向的多通道滤波器，即为Gabo小波，其表达式如下：

${{g}_{mn}}={{a}^{-m}}g\left( {{x}^{'}},\text{ }{{y}^{'}} \right),\text{ }a＞1$

(2)

式中:$\left[ {\begin{array}{*{20}{c}} {x'}\\ {y'} \end{array}} \right] = {a^{ - m}}\left[ {\begin{array}{*{20}{c}} {{\rm{cos}}{\theta _n}} & {{\rm{sin}}{\theta _n}}\\ { - {\rm{sin}}{\theta _n}} & {{\rm{cos}}{\theta _n}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} x\\ y \end{array}} \right],{\theta _n} = n{\rm{ \mathsf{ π} /}}K$，n为方向数, n∈[0, K－1], K代表Gabor总的方向数；a^－m为尺度因子, m为尺度数.

Gabor滤波器的位置由n方向数和ω高斯函数的复调制频率2个参数决定.一般取n=4，即4个方向分别为：0°，45°，90°，135°，基本可以搜索到磁异常全方位的纹理特征信息.ω高斯函数的复调制频率的选择即波长λ的选择，采用Jain等^[11]给出的算法：对于异常数据，大小为N_r×N_c，则λ取${{N}_{\text{c}}}\sqrt{2}/{{2}^{2}},{{N}_{\text{c}}}\sqrt{2}/{{2}^{3}},\cdots ,{{N}_{\text{c}}}\sqrt{2}/{{2}^{\log _{2}^{\left( {{N}_{\text{c}}}/4 \right)}}}$, 其中N_r，N_c分别为数据的行数和列数.根据以上设置，可得到4个不同方向上log₂(N_c/4)1个不同尺度的Gabor滤波器组.

设磁异常值I(x, y)，利用不同尺度和方向的Gabor滤波器g_mn(i, j)对图像I(x, y)进行特征提取，其中(i, j)为滤波器元素，则

${{\mathit{\boldsymbol{W}}}_{mn}}\left( x,\mathit{y} \right)=\sum\limits_{i=1}^{M}{\sum\limits_{j=1}^{N}{I\left( i,\mathit{j} \right){{g}_{mn}}\left( i,\mathit{j} \right)}}$

(3)

通过一组m个不同尺度、n个不同方向的多通道Gabor滤波器对图像进行特征提取，得到纹理特征在每个像素点都有一个m×n维的特征向量W.

1.2 灰度共生矩阵

灰度共生矩阵(GLCM)是图像中2个像素灰度级联合分布的统计形式，是统计2个像素点在一定距离和一定方向上的联合概率分布，能较好地反映纹理灰度级空间相关性的规律.灰度共生矩阵定义为^[12]：设图像某一区域有N个灰度值，则对应区域的灰度共生矩阵为N×N阶的矩阵P_N×N，在矩阵位置P(i, j)(i, j=1, …, N)处元素是从灰度为i的像元离开某个固定位置关系δ=(D_x, D_y)处像元灰度j的概率，D_x和D_y称为偏移距离，δ称为位移量，在实际的图像纹理特征提取中δ通常选取0°，45°，90°，135°四个方向.

本文利用灰度共生矩阵提取纹理特征，在该过程中，采用滑动窗口技术，将磁异常平面图分成大小相等、互不重叠的小块，然后计算每块的灰度共生矩阵.在实际图像纹理特征提取中，通常不直接使用灰度共生矩阵进行纹理分析，而是利用灰度共生矩阵计算得到的纹理特征统计量.Haralick等^[13]提出14种由灰度共生矩阵计算出的纹理特征统计量，本文选用最常用的4个：能量(E_n)、相关性(C_or)、对比度(C_on)、逆差矩(I_dm)，其表达式如下：

${{E}_{n}}=\sum\limits_{i=0}^{N-1}{\sum\limits_{j=0}^{N-1}{P{{\left( i,\text{ }j \right)}^{2}}}}$

(4)

${{C}_{\text{or}}}=\frac{\sum\limits_{i=0}^{N-1}{\sum\limits_{j=0}^{N-1}{ijP\left( i,\text{ }j \right)-{{\mu }_{\text{x}}}{{\mu }_{y}}}}}{\sigma _{\text{x}}^{2}\sigma _{y}^{2}}$

(5)

${{C}_{\text{on}}}=\sum\limits_{i=0}^{N-1}{\sum\limits_{j=0}^{N-1}{{{\left( i-j \right)}^{2}}P\left( i,\text{ }j \right)}}$

(6)

${{I}_{\text{dm}}}=\sum\limits_{i=0}^{N-1}{\sum\limits_{j=0}^{N-1}{P{{\left( i,\text{ }j \right)}^{2}}}}/\left[ 1+{{\left( i-j \right)}^{2}} \right]$

(7)

式中:μ_x, μ_y是均值，${{\mu }_{x}}=\sum\limits_{i=0}^{N-1}{\sum\limits_{j=0}^{N-1}{iP\left( i,\text{ }j \right)}},{{\mu }_{y}}=\sum\limits_{i=0}^{N-1}{\sum\limits_{j=0}^{N-1}{jP\left( i,j \right)}}$, σ_x, σ_y是方差, ${{\sigma }_{x}}=\sum\limits_{i=0}^{N-1}{\sum\limits_{j=0}^{N-1}{{{(i-{{\mu }_{x}})}^{2}}}}P\left( i,\text{ }j \right)\text{, }{{\sigma }_{y}}=\sum\limits_{i=0}^{N-1}{\sum\limits_{j=0}^{N-1}{{{(j-{{\mu }_{y}})}^{2}}P\left( i,\text{ }j \right)}}$.

1.3 纹理特征聚类分区

利用Gabor滤波器和灰度共生矩阵提取磁异常的纹理特征信息后，则可对提取的纹理特征数据进行聚类分区.本文采用算法成熟且对大型数据处理效率较高的K-均值聚类方法，聚类数的选择则采用OKMS算法^[14-16]，该算法自动搜索最佳的聚类数，对初始的聚类中心，使提取的特征信息矩阵到其所属类别的中心的距离平方和最小以达到最佳聚类.Gabor滤波器提取的不同方向和尺度的综合纹理信息可表示为特征向量{G_1j, G_2j, …, G_lj}，灰度共生矩阵提取的能量、相关性、对比度和逆差距可表示为特征向量{T_1j, T_2j, …, T_lj}，其中l为纹理特征向量数据集的长度，j=1, 2, …, m_s为每个数据集中纹理特征量的个数.聚类过程对Gabor滤波器和灰度共生矩阵提取的纹理特征向量数据集{G_1j, G_2j, …, G_lj; T_1j, T_2j, …, T_lj}进行，其计算步骤如下：

(1) 对Gabor滤波器和灰度共生矩阵提取的纹理特征向量{G_1j, G_2j, …, G_lj; T_1j, T_2j, …, T_lj}，选k个样本作为初始聚类中心(Z₁, Z₂, …, Z_k).

(2) 对每个样本找到离它最近的聚类中心Z_v(v=1, 2, …, k)，并将其分配到Z_v所标明的类U_v(v=1, 2, …, k).

(3) 采取平均的方法重新计算分类后的各聚类中心.

(4) 计算样本到聚类中心的距离$D=\min {{\sum\limits_{i=1}^{n}{\left[ d{{\left( {{T}_{\text{i}}},{{Z}_{v}} \right)}^{2}}+d{{\left( {{G}_{\text{i}}},{{Z}_{v}} \right)}^{2}} \right]}}_{v=1,2,\cdots ,k}}$, d为欧式距离.

(5) 如果D值满足收敛条件，则将聚类结果输出并终止计算，否则转至步骤(2) 继续计算.

其中聚类数k事先无法确定, 这里采用OKMS算法来确定最佳的聚类数, 其算法如下：

(1) 选择聚类数的搜索范围[k_min, k_max]，其中k_min=2, ${{k}_{\max }}=\operatorname{int}\left( \sqrt{n} \right)$.

(2) 对于不同的聚类数k从k_min到k_max，按照上述K均值聚类算法，得到聚类结果，并计算聚类结果的Silhouette指标.

(3) 比较Silhouette指标，指标值达到最大时所对应的k即为最佳聚类数k_opt.

(4) 将最佳的聚类数k_opt作为上述K-均值聚类算法的聚类数k.

2 模型试验

海洋磁条带异常记录了海底扩张的时代及构造信息.海盆内的条带磁异常常会受到断裂、海山及后期岩浆活动的改造，形成不同的磁异常纹理特征，给磁异常的追踪和识别带来困难.采用纹理特征提取方法可将磁条带异常区和受岩浆活动干扰的区域圈定出来，以便对磁异常进行进一步的处理分析.为了检验该方法对不同磁异常区的识别圈定能力，设置的模型如图 1a，模型包含走向存在差异的磁条带异常区、受海山污染的磁条带异常区及火山岩分布异常区.模型参数为：网格数为101×101，网格间距1 km；顶底界面埋深3~7 km，磁化强度幅值为1~10 A·m^-1.图 1a表示模型场源的位置，图 1b为模型产生的理论磁异常.

利用Gabor滤波器对模型磁异常进行纹理特征信息的提取，Gabor滤波器尺度为381×381, 197×197, 95×95, 49×49，分别提取0°, 45°, 90°, 135°四个方向的纹理特征.图 2为Gabor滤波器提取的磁异常不同尺度、不同方向的纹理特征信息(从左到右依次为0°, 45°, 90°, 135°).

利用灰度共生矩阵方法对模型磁异常进行纹理特征信息的提取，采用窗口大小为64×64, 32×32, 16×16, 8×8的滑动窗口计算灰度共生矩阵，其中每个窗口计算0°, 45°, 90°, 135°四个方向的灰度共生矩阵, 然后取平均值，由得到的灰度共生矩阵计算能量、相关性、对比度、逆差矩等纹理特征向量.图 3为灰度共生矩阵4个不同窗口提取的纹理特征(从左到右依次为能量、相关性、对比度、逆差矩).

根据Gabor滤波器和灰度共生矩阵提取的纹理特征向量，采用上述方法对纹理特征进行聚类分区. 图 4a为纹理分析计算得到的模型纹理分区结果，图 4b为分区结果在模型磁异常图上的投影.从分区结果可以看出，模型磁异常被分为了4个区块，这4个区块分别对应岩浆活动区(A)、受岩浆活动污染的磁条带异常区(B)、磁条带异常区(C)、磁条带异常区(D).分区结果显示磁条带异常区和岩浆活动区被准确地圈定出来，区块位置与模型场源位置吻合较好.

3 实际应用

加罗林海域位于太平洋板块、菲律宾海板块和印度澳大利亚板块交汇处，是西太平洋的重要组成部分.加罗林板块与周围板块之间存在张扭性、挤压性等复杂边界，海盆磁条带异常受欧里皮克海隆、西加罗林海脊两翼和断裂带干扰的影响，给识别造成困难^[17].应用本文提出的磁异常图像纹理分区方法，对加罗林海域磁异常进行分区.从图 5分区结果图上可以看出，加罗林海域磁异常被分割为了11个区块.将分区结果投影到加罗林海域磁异常图上，如图 6，可以看出磁条带异常区被圈定出来.表 1给出了各分区区块的位置及磁异常纹理特征.从分区结果可以看出，该方法可较准确地将磁条带异常区圈定出来，如西加罗林海盆磁条带异常区Ⅱ，Ⅲ和东加罗林海盆磁条带异常区Ⅳ，Ⅴ，其中西加罗林海盆磁条带异常以西加罗林海槽为界分为Ⅱ，Ⅲ两部分，分别位于海槽两侧；东加罗林海盆磁条带异常以科尔斯哥德海槽为界分为Ⅳ，Ⅴ两部分，分别位于海槽两侧.欧里皮克海隆异常区(Ⅰ)是后期岩浆沿岩石圈断裂喷发形成，东加罗林海盆南部异常区(Ⅵ)可能与后期岩浆活动有关，无磁条带异常.加罗林海脊异常区(Ⅶ)磁异常特征与太平洋板块向马里亚纳海沟俯冲引起的断裂有关.阿玉海槽异常区(Ⅷ)南北向带状磁异常与海槽的扩张有关.新几内亚岛弧异常区(Ⅸ)磁异常呈低幅值的片状，与邻近的线性磁异常明显不同.俾斯麦海北部磁异常区(Ⅹ)呈斑状，这里有许多海山出露.太平洋西部异常区(Ⅺ)磁异常呈团块状，此处岩浆活动强烈，有许多海山.从图 6加罗林海域磁异常分区图可以看出，本文提出的磁异常图像纹理分区方法可以将磁条带异常区、岩浆活动区等不同纹理特征的磁异常区圈定出来，说明该方法对海洋磁异常的分区具有较好效果，具有良好的应用前景.

4 结论

(1) 模型试验表明，利用Gabor滤波器和灰度共生矩阵相结合的图像纹理分区方法可以有效提取海洋磁异常图像的纹理特征，尤其是对条带状磁异常的识别圈定效果较好.

(2) 利用OKMS算法自动搜索聚类数，可以给出较合理的分区结果，从而避免因人为给出聚类数的不同而使得磁异常的分区结果出现较大差异.

(3) 采用本文的纹理分区方法对加罗林海域磁异常进行分区，结果表明加罗林海域磁异常可分为11个磁异常区块，相应圈定了磁条带异常分布区和岩浆活动区，为磁条带异常的追踪、识别和后续的处理工作奠定了基础.