岩石有限应变测量作为定量化研究岩石变形的重要手段，在构造地质学研究中扮演着重要的角色。由于其工作量大、工作时间长、枯燥等问题，使得在进行大区域应变分析时，采样密度偏低^[1]。近些年，对颗粒边界进行识别的方法逐渐被引入到岩石有限应变测量工作中。它们主要集中于边缘检测和单个图像的处理上^[2~5]。通过灰度图像边缘检测产生的结果，不能直接应用于应变测量，而且受到特殊的设备^[2]、前期或者后期的图像处理^{[3, 6]}等限制。在单个图像中，当颗粒之间颜色差别比较大时，区分颗粒比较容易；而当边界颜色较暗、与颗粒相近，或在消光位时，将无法把需要测量的颗粒与其他颗粒进行区分。同时矿物存在波状消光也使得单个图像的识别更加困难。多重图像限制的自动种子区域生长算法(Multiply Images Constrained Automated Seeded Region Growing, MICASRG)利用正交偏光显微镜下石英砂岩薄片不同角度的干涉图像进行图像综合增强，然后修改CASRG算法^[1]，对综合图像进行识别，获得了较详细的颗粒信息。此外，该方法的后期工作不需要人工干预，可直接获得各个颗粒的几何参数用于应变测量。

1.   原理与方法实现

多重图像种子约束自动区域增长是将偏光显微镜下同一视域内，不同消光位置的多个图像进行综合，然后根据种子约束自动区域增长的方法对综合后的图像进行识别，达到获得尽可能多的信息的目的。

1.1   图像获取及预处理

正交偏光显微镜下获得的多个图片进行配准，是进行颗粒识别的前提。单偏光显微镜下，相邻颗粒的边界不清晰，许多矿物颗粒具有相似的颜色，以致区分相邻的不同颗粒难度加大。在正交偏光显微镜下，不同的矿物，同种矿物不同的光率体方位、不同的薄片厚度等干涉色均不相同。通过旋转物台还可以获得同一颗粒有规律变化的彩色干涉图像。目前，显微镜下的数字图片较容易获得，这就为多个图像进行叠加综合提供了可能。本文使用了OLYMPUS-BX51型数码相机。由于考虑到程序运行速度和矿物在偏光显微镜正交光下出现消光的问题，笔者选取了两张呈45°的镜下图像。由于在不同的光照条件，同一颗粒的颜色、亮度不相同，所以在进行图像获取时，应使得物台旋转至各个角度时的光照条件相同。

1.2   图像叠加原理

不同方位的颗粒，进行图像差值运算可以使颗粒内部颜色均匀，增大相邻颗粒的颜色信息差别。任意一张正交偏光显微镜下的干涉图像，都是一张包含红色、绿色和蓝色3个标准元色的普通图片。由于石英砂岩薄片内绝大部分矿物颗粒为石英，所以干涉色比较低，颜色信息较少(见图 1)。而将任意矿物的光率体旋转45°后，大多数矿物颗粒的干涉色和亮度会发生变化。对于同一颗粒内部，干涉色的变化是均匀的，即使存在波状消光和双晶；而不同颗粒的两组图像干涉色差值变化较大，这为图像进行差值运算提供了前提。根据以上的算法，提出了合成图像信息的计算公式：

图  1  原始图像

Figure  1.  The original image

下载: 全尺寸图片 幻灯片

上式中， $\overrightarrow {{\rm{New}}} {\rm{(x,y),}}\overrightarrow {{\rm{p1}}} {\rm{(x,y),}}\overrightarrow {{\rm{p2}}} {\rm{(x,y)}}$ 分别为合成图像、原图像、旋转后图像各个点的红、绿和蓝色的颜色信息，ABS为取绝对值函数，用于确定两图像同一位置颜色信息的变化。如果旋转后图像的颗粒边界仍然不清晰或颗粒颜色信息变化范围较小，则对两个图像同一位置的颜色取其最大值，即MAX函数。这样，就可以对颜色较暗且颜色变化不大的颗粒和边界进行有效的区分。本算法可使边界清晰，颗粒内部颜色更加均匀，提高颗粒边界的对比度(见图 2)，其计算公式如下：

图  2  经过多重图像叠加边缘增强后图像

Figure  2.  The image after superposition and boundary enforced

下载: 全尺寸图片 幻灯片

其中，MAX为取两组数据绝对值函数。

1.3   合成图像边缘增强与识别

为了能够更好地区分颗粒边界，笔者对图像的红色、绿色和蓝色数值分别进行卷积运算，增强颗粒内部与边界的对比度。首先用3×3的高斯滤波器分别对图像各个颜色分量进行滤波，以去除图像中的噪声；然后对图像的每个像素，计算其梯度的大小G和方向θ^[7~8]；最后对梯度图像进行非极大抑制，二次阀值运算，将两次阀值的图像和红色、绿色和蓝色分量的各个图像进行合成，加强边缘(见图 2)，得到合成后图像边缘的精确位置 ${\overrightarrow {{\rm{Px}}} }$ 和 ${\overrightarrow {{\rm{Py}}} }$ 。笔者修改了CASRG算法^[1]，并将修改后算法应用于合成后的图像识别。该方法应用欧式距离和自动阀值进行颗粒的识别，并且可以在一定程度上绕过边缘增强方法在颗粒内部产生的裂隙。具体计算公式如下：

2.   结果验证与比较

为了验证该方法得到的颗粒参数能否应用于岩石有限应变测量，本文分别从几何学参数和各种岩石有限应变测量方法上，将MICASRG方法与传统的手工描图的结果进行了比较(见图 3、图 4)。为了能够取得更加可信的结果，本文选择了156个颗粒进行分析对比。

图  3  颗粒自动识别与手工描图结果对比

Figure  3.  Contrast of results between automatic identified and handcrafted

下载: 全尺寸图片 幻灯片

图  4  自动识别与手工描图几何参数的对比

Figure  4.  The contrast of geometric parameter between MICASRG and hand crafted

下载: 全尺寸图片 幻灯片

2.1   几何学方法的对比

矿物颗粒的几何学参数包括颗粒面积，质心坐标，长、短轴，方位等特征，这些参数可以描述矿物颗粒的外形、大小、方位以及颗粒经历的变形历史和强度。通过比较MICASRG方法与传统的手工方法得到的颗粒参数，观察它们之间的拟合程度，从而正确反映该方法的可行性和精度。

对比颗粒长轴长度与短轴长度的频数统计直方图(见图 5)可得，它们的分布相对较集中。长轴数值主要分布于20~40之间；短轴主要分布于10~30之间。分布的相对稳定表明该区的变形较稳定，也从一个侧面反映了使用156个颗粒的岩石有限应变测量方法是可行的。

图  5  自动识别颗粒长轴长、短轴长频数统计直方图

Figure  5.  Histogram of length of the major axis and axis by MICASRG

下载: 全尺寸图片 幻灯片

在方位识别方面，两种方法的相关性系数较低，仅为0.65。引起这么大差距的原因可能有2点：① 颗粒方位上-90°与90°是同一方位，即位于坐标轴4个角的点，可以投影到同一区域内。经变换后，相关系数提高到0.85；② 从图 5可知该区变形强度较弱，矿物没有明显的拉伸与挤压，这增大了判断颗粒方位的主观性，尤其是在传统的手工描图方法上，这种影响表现的更为突出。

2.2   岩石有限应变测量结果的比较

目前存在很多岩石有限应变测量的方法^[9~12]，本文应用R_f/φ法(R_f为变形后椭圆轴率，φ为初始椭圆与变形椭圆之间的夹角)，分别对2组数据进行了岩石有限应变测量，并得出了相似的结果(见表 1)。表 1中长轴方向差别比较大，这是由颗粒本身应变量比较小引起的。长短轴差别较小，使得手工描图无法准确地判断长轴方向，从而造成长轴方位识别上的困难。

表  1  两组数据的R_f/ $\varphi $ 法岩石应变测量结果

Table  1.  Results of finite strain estimates (R_f/ $\varphi $ ) based on two sets of data

参数	手工描图	自动识别
Rs	1.34	1.37
长轴方向	14.25	6.84
Ri	1.66	1.73
注：Rs为应变椭圆轴率，Ri为初始椭圆轴率

下载: 导出CSV 

| 显示表格

3.   结论

传统的手工描绘进行岩石有限应变测量的方法需要耗费大量的人力、物力和时间，而且存在很大的偶然误差。前人利用边缘检测和单个图像的处理获得矿物颗粒的信息，但往往受到一定的限制。本文利用正交偏光显微镜下石英砂岩薄片不同角度的干涉图像，通过修改CASRG算法对图像进行识别，获得更多和更详细的信息。同时，后期工作不需要进行人工干预，可直接获得各个颗粒的几何参数，进行岩石有限应变测量。

将本文测得的156个颗粒的结果与传统方法进行几何参数和岩石应变测量结果进行对比，得到了近似的结果。质心坐标、面积、长短轴识别的相关系数较高，达到0.92以上；岩石有限应变测量的结果也很相似。在几何参数和岩石应变测量结果中长轴方位的相关系数较低，这可能是由于本次试验中所采用的石英砂岩岩石薄片变形较小，很多颗粒的长短轴接近等长所引起的。与其他的岩石有限应变测量自动化方法相比，本方法综合多个图像，不受特殊设备的限制，利用彩色信息对各个颗粒进行识别，取得了令人满意的效果。在具体操作中笔者为了快速检验算法的可行性，采用了Matlab语言，造成识别速度比较慢。若采用C语言或者C++等语言，这一问题可以避免。