分形的本质特性是自相似性^[7], 而成矿过程中各种尺度地质异常的大小和个数通常具有自相似性^[8], 因而地质异常也具有分形特征。通用的分形模型^[9] 如下。设有数据集合{X_i, i=1, 2, …, N}, 则

N(r)= ∑i=1NX_i=Kr^-D , (1)

式中: r为特征尺度; K为常数; D为分维数; N(r)为尺度大于等于r的X_i的数目。

将式(1)两边取对数可得

lg N(r)=lg K-D lg r 。 (2)

由式(2)可以看出, lg r与lg N(r)呈线性关系。计算时, 可以遥感蚀变分量图像的像元亮度值为r值, 则N(r)为亮度值大于等于r的X_i的数目。

地质异常还具有突变性特征^[6], 变点是指输出序列在某未知时刻发生的突然变化, 因此变点分析可定量反映地质异常的突变特征。通过变点统计模型可有效检验样本序列是否存在变点, 从而判断变点的位置与个数^[10]。采用的均值变点模型^[11]如下。

1)令i=2, …, N, 对每个i将样本分为2段: x₁, x₂, …, x_i-1和x_i, x_i+1, …, x_N。分别计算每段样本的算术平均值( Xi1—和 Xi2—)及方差(S_i), 即

S_i= ∑t=1i-1(X_t- Xi1—)²+ ∑t=iN(X_t- Xi2—)² ; (3)

2)计算整段的算术平均值( X̅)及方差(S), 即

X̅= ∑t=1NX_t/N , (4)

S= ∑t=iN(X_t- X̅)² ; (5)

3)计算期望值(E), 即

E(S-S_i)=E[N^-1(i-1)(N-i+1)( Xi1—- Xi2—)²], i=2, 3, …, N 。 (6)

由上可见, 样本序列如果存在变点, 那么变点会使S和S_i的差值增大; 而S是固定值, 因此最小的S_i值即对应存在的变点。计算时, 可先对lg N(r)与lg r的商取对数, 以构建样本序列X₀; 再计算X₀的统计量S_i。

实验区选择在甘肃北山方山口地区。区内出露地层主要为中泥盆统三个井组(D₂s), 岩性为浅变质中细粒长石石英砂岩、长石砂岩夹粉砂岩和泥质粉砂岩。断裂构造以近EW向压扭性断裂为主, 沿走向局部向NEE向偏离, 形成弧形构造线。区内岩浆活动强烈, 以华力西期为主, 次为印支期, 从基性―酸性均有产出。其中, 辉长岩-辉绿岩呈岩株状产出, 石英闪长岩呈岩基出露于研究区中南部, 二长花岗岩呈NE向不规则侵入于石英闪长岩及三个井组地层中; 花岗斑岩呈岩脉状或不规则岩株状产出, 穿插于石英闪长岩中。区内主要分布热液型矿床, 尤以金矿化为主。区内存在的蚀变类型有硅化、黄铁矿化(褐铁矿化)、绢云母化、铁碳酸盐岩化和青磐岩化等^[12]。

ASTER传感器搭载的Terra卫星发射于1999年12月, 其数据具有3个可见光和近红外波段(VNIR, 0.52~0.86 μ m)、6个短波红外波段(SWIR, 1.6~2.43 μ m)和5个热红外波段(TIR, 8.125~11.65 μ m), 较宽的光谱覆盖范围可供提取多种遥感蚀变矿物信息。

ASTER数据获取时间为2000年10月16日, 数据级别为L1B, 已经过辐射校正和几何纠正, 无云及阴影遮盖, 影像质量较好。本次实验中对该图像数据进行了串扰校正^[13](串扰是成像光谱仪非理想成像特性中的一种, 分为电串扰和光串扰, 广泛存在于成像光谱仪中)、大气校正和重采样处理(图像分辨率重采样至30 m), 并选择对河流、植被及第四系等干扰地物敏感的波段, 通过掩模处理降低其对异常信息提取的影响。

对于遥感图像的混合像元分解问题, 研究人员提出了大量的光谱解混算法。其中, 自适应一致估计法(adaptive coherence estimator, ACE)是一种针对遥感图像非结构化目标探测的算法, 它将像元背景视为符合多元高斯分布的统计模型, 将背景和噪声统一视为背景进行建模, 通过贝叶斯统计计算概率得分来识别目标像元^[14], 其优点是可以基于较少的已知样本就能有效检测目标^[15]。

ACE是广义似然比探测法(generalized likelihood ratio test, GLRT)的一种自适应版本^[14]。假设像元背景是一个混合模型, 且背景b服从多元正态分布b~N(μ ₀, Γ ), 将均值μ ₀从观测值x中移去后, 则有

H₀: x=b, H₁=S· a+σ ²· b , (7)

式中: H₀和H₁为2种假设地物; H₀和H₁的分布分别为x~N(0, Γ )和x~N(S· a, σ ²· Γ ), 即背景在2种假设下具有相同的协方差结构, 但有着不同的方差, 该方差直接关系到目标物体占像元区域的比例。在上述情况下, ACE算法可表达为

D_ACE(x)= xTΓ-1S(STΓ-1S)-1STΓ-1xxTΓ-1x> H1< H0η _ACE , (8)

式中: Γ 为图像协方差矩阵; x为光谱向量; S为目标向量。

本次实验提取了褐铁矿、绢云母和绿泥石3种蚀变矿物异常信息。实验中均采用ACE目标检测算法, 以USGS矿物标准波谱库为基础, 通过ENVI软件提取3种蚀变异常图像, 并将提取结果的亮度值拉伸到0~255范围; 然后应用分形模型计算异常图像的双对数图, 以分析其自相似性; 并在此基础上采用均值变点法分析异常的突变性特征, 即通过计算获取最小的S_i值。

对于褐铁矿异常, 在i=59时, S_i值最小, 对应的像元亮度值r=60, 故r=60可作为地质背景和褐铁矿Ⅲ 级异常的临界点; 继续利用上述方法对r≥ 60的数据进行计算, 分别得到r=111时的Ⅱ , Ⅲ 级异常的临界点和r=119时的Ⅰ , Ⅱ 级异常的临界点。图1(a)展示了地质背景和蚀变异常的自相似性分布特征以及各自的自相似性拟合直线, 图1(b)展示了Ⅰ , Ⅱ , Ⅲ 级异常的临界阈值及各自的自相似性拟合直线, 不同颜色线段汇交点为各级异常临界阈值。

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 褐铁矿蚀变分形特征分析图Fig.1 Diagram of fractal characteristics analysis for limonite alteration

从图1可以看出, 褐铁矿蚀变异常段中3段直线方程的拟合率均较好, 表明研究区内褐铁矿蚀变异常区域具有多重分形特征, 不同程度的蚀变异常对应不同的分维值D(图1(b))。根据确定的临界阈值对蚀变异常图像进行密度分割, 即可制作出褐铁矿蚀变异常分级图。

对于绢云母异常, 在i=73时, S_i值最小, 对应的像元亮度值r=74, 故r=74可作为地质背景和绢云母Ⅲ 级异常的临界点; 继续利用上述方法对r≥ 74的数据进行计算, 分别得到r=125时的Ⅱ 和Ⅲ 级异常的临界点以及r=132时的Ⅰ 和Ⅱ 级异常的临界点。图2为绢云母蚀变异常灰度值分形特征分析图, 可以看出异常段中3段直线方程的拟合率均超过0.91, 表明研究区内绢云母蚀变异常也具有多重分形特征, 不同程度的蚀变异常对应不同的分维值D(图2(b)), 不同颜色线段汇交点为各级异常临界阈值。根据确定的临界阈值可制作出绢云母蚀变异常分级图。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 绢云母蚀变分形特征分析图Fig.2 Diagram of fractal characteristics analysis for sericite alteration

对于绿泥石异常, 在i=169时, S_i值最小, 对应的像元亮度值r=170, 故r=170可作为地质背景和绿泥石Ⅲ 级异常的临界点; 继续利用上述方法对r≥ 170的数据进行计算, 分别得到r=221时的Ⅱ , Ⅲ 级异常的临界点和r=229时的Ⅰ , Ⅱ 级异常的临界点。图3为绿泥石蚀变异常灰度值分形特征分析图, 可以看出异常段中3段直线方程的拟合率均超过0.98, 表明研究区绿泥石蚀变异常同样具有多重分形特征, 不同程度的蚀变异常对应不同的分维值D(图3(b)), 不同颜色线段汇交点为各级异常临界阈值。根据确定临界阈值可制作出绿泥石蚀变异常分级图。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 绿泥石蚀变分形特征分析图Fig.3 Diagram of fractal characteristics analysis for chlorite alteration

图4为根据上述临界阈值制作的实验区褐铁矿、绢云母和绿泥石遥感蚀变异常分布图。

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 实验区遥感蚀变异常分布图Fig.4 Distribution of remote sensing alteration anomalies in experimental area

从图 4可以看出, 蚀变矿物信息呈近EW向线性展布, 与区域构造线方位基本一致, 区内已知热液型金属矿(化)点几乎都落在遥感蚀变异常区内或其边缘; 同时, 遥感蚀变异常的分布与已有的化探异常和重砂异常有较好的吻合性。