CSD SOFTWARE DESIGN FOR FINITE-STRAIN DETERMINATION WITH INVERSE SURFOR WHEEL
-
摘要: 用计算机实现反向轮法测岩石有限应变的过程, 可提高有限应变测量的准确性和效率.开发出了用Visual Basic 5.0编写的CSD软件, 利用岩石薄片图像或者矿物颗粒轮廓图测出岩石有限应变的大小及应变椭圆长轴方向, 同时测出矿物颗粒分布的优选方位.其操作相当简便, 首先统计出图像中0°~180°各方向上矿物颗粒的边界数目; 然后在方位-边界坐标系中投点, 并利用最小二乘法进行数据点的多项式曲线拟合; 最后求出曲线的极值点坐标, 并根据坐标绘制相应的应变椭圆.软件运行中, 统计出的边界数据和拟合曲线以及应变椭圆图形都是可视的, 并能进行相应的保存.Abstract: Computer simulation of finite-strain determination with the inverse SURFOR wheel may improve the preciseness and efficiency of the strain determination. The CSD softwre, programmed in Visual Basic 5.0 and very simple for its friendly interface, can be used to measure the size of the finite strain and the direction of the long axis of the finite strain ellipse in deformed rocks by means of section images or mineral grain sketches. At the same time, this software can also be used to measure the preferred orientation of the mineral grain distribution. Firstly, the CSD software collects boundaries of mineral grains in the image at all directions from 0°-180°. Secondly, the boundary data dots are distributed in the orientation-boundary coordinate system, and the polynomial curve fitting is made by means of the least-squares. Lastly, the extreme point coordinate of the fitted curve is obtained, and the corresponding strain ellipse is drawn up in line with the coordinate. During the operation of the CSD software, all the collected data of mineral grain boundaries and the resulting fitted curve, and the strain ellipse images can be both visible and saved in the corresponding files.
-
随着计算机和计算机图形学的发展, 地质构造的一些测量已经由手工向计算机自动测量发展.如AMIN用计算机实现Fry法应变测量点图(见文献[1]), Allard等[2]用计算机在岩石薄片图像上实现优选方位的自动测量.本文应变自动测量基于的反向轮法[3]是由Panozzo在1987年提出的, 由于当时是手工进行测量, 所以未得到充分利用.其基本理论是: 岩石在未变形时各种组构是随机的, 各向同性的, 但变形后将出现一定的优选方位, 在给定截面上可得有限应变椭圆.假设变形前岩石中矿物颗粒随机分布, 则给定线段在任一方向与矿物边界相交机率不变, 即
式中: N为变形前单位线段与颗粒边界交点个数, φ为线段与参考轴x之间的夹角.岩石变形后, 所选线段与矿物颗粒边界交点个数不变, 但线段长度以及交点间的距离发生了改变, x轴变为x′轴, 所给线段与x′轴的夹角为φ′.当应变椭圆长轴与x′轴不平行时, 有一夹角φi, 则单位长度线段与颗粒边界交点数为[3]:
其中
和 是倒易主长度比, 这实际上是一个椭圆公式, 于是, 应变椭圆轴比可表示成:应变椭圆长轴方向为: φi=φ′min.
实际手工操作时, 在一个圆内作一系列平行线制成轮状, 测量时将轮心固定在x′轴原点, 每次反向转动一个角度(图 1), 平行线与x′轴的交角为φ′, 统计对应不同φ′角时平行线与矿物颗粒边界交点总数, 直到转动180°, 这便得到一系列数据, 在平面直角坐标系中, 作n′(φ′)对φ′的投点, 并进行曲线拟合.曲线最高点对应n′(φ′)max, 最低点对应n′(φ′)min.
CSD软件就是用计算机来实现反向轮法测量岩石有限应变的这一系列手工操作, 以实现应变测量的自动化和提高数据统计的精度和测量效率.
1. 软件设计
1.1 边界统计
本软件的主要步骤之一就是统计图像上0°~180°各方向颗粒边界数目.在该步骤中怎样判别矿物边界是关键问题.当前, 关于计算机判别颗粒对象边界的文献较多[4, 5].这些文献主要是针对颗粒对象进行边界提取, 而本文中的边界判断是对任一直线方向上的矿物颗粒边界进行判断, 其原理如下:
岩石矿物颗粒一般是紧密接触的, 在薄片中通常也是如此.由于相邻矿物的种类和光性不一样, 所以在正交显微镜下相邻矿物的颜色是不同的, 这一点就为我们提供了判别边界的依据.如图 2, 计算机沿直线AB方向读取点(像素点)的颜色可以利用颜色分解函数把读取到的颜色分解为RGB三原色, 如果在直线AB方向上前后相邻两点(像素点)的RGB三原色值中存在一个大于10的差值(这是为了排除矿物颗粒内部干涉色不均引起的误差, 实际相邻矿物的RGB三原色值一般至少会存在一个大于10的差值), 即判断存在一个矿物颗粒边界.当判断出一个边界后, 须忽略掉5个点(像素点)的判断结果, 即为了避免重复判断同一边界, 实际颗粒边界一般是小于5个像素点.但值得注意的是: 沿直线方向读点, 不能用数学上三角函数来计算像素点的坐标, 因为那样会产生较多间隙, 很多点会读取不到.本文采用了Bresenham直线算法[6]来读取像素点的坐标, 该算法可以连续画点而绘出直线, 但是作者并没有用其画点功能, 而是把画点改进为读点的坐标, 以达到不间断地读取直线上像素颜色的目的.
为了统计出图像上0°~180°各方向颗粒边界的数目, 须进行以下几步: (1)建立坐标系统(图 3), 原点(0, 0)定在反向轮中心, 半径R由矩形图像短边长度的一半确定. (2)用一些简单的三角函数便可从反向轮的半径和直线数目算出各直线与x′轴平行时的初始端点坐标. (3)由各直线初始端点坐标可求出反向轮中平行线反向转动φ′角度后的端点坐标. (4)有了φ′角度下各直线端点的坐标, 便可以利用改进的Bresenham算法函数, 比较相邻两点在RGB三原色上的差值来判断是否为边界.最后累加各直线方向上的矿物边界数目, 得到φ′角度对应的矿物边界总数n′(φ′). (5)从0°到180°之间递增变化, 且循环(3)和(4)步骤便可得出0°到180°各方向上的边界统计结果.
为了实现统计意义, 软件中设置了确定反向轮直线数目和每次反向转动角度大小的对话框, 可以根据输入的薄片图像大小(该软件可以接受WINDOWS的BMP格式、JPEG格式和GIF格式的图像)、颗粒粗细来设置反向轮直线数, 一般图像越大颗粒越细, 则直线数少; 图像越小颗粒越粗, 则直线数多.
1.2 最小二乘法曲线拟合
Hart等[7]对各向异性极点数据进行过椭圆拟合; Thompson等[8]用最小二乘法进行过多项式曲线拟合.本文是对0°~180°各方向上边界数据进行投点、最小二乘法多项式曲线拟合, 并求出曲线的最大值与最小值坐标.下面是该软件测出的一组数据(表 1)和对它进行的最小二乘法多项式拟合曲线(图 4).
表 1 CSD软件统计出的边界数据和方向Table Supplementary Table The boundary data and orientation collected by the software CSD关于边界数目和方向的增广矩阵A(N)是从最小残差平方和理论[7]
求得, 其中(xi, yi)是数据点中的第i点, 如M > N则有拟合多项式
.其中, M表示数据点的数目, N表示多项式拟合次数.得出关于边界数目和方向的增广矩阵A(N)后, 利用高斯-亚当解方程组法便可解出方程组, 得到多项式的系数.但值得注意的是, 该过程中运算数据的精度要为双精度型, 这是由于xi是从0°到180°之间变化, A(N)中的有些高次项可能会很大, 产生溢出错误.本文则采取了一种先缩小再放大的方法来解决数据过大溢出的问题.在计算A(N)中每一项时, 都先把xi缩小10倍, 即xi从0°到18°取值, 这样便不会产生溢出问题.当在方向-边界坐标系中投点及拟合曲线时, 再把xi放大10倍, 让xi还原到从0°到180°.由上面的拟合曲线也可看出先缩小再放大的方法是可行的, 有效的.当拟合出多项式曲线后, 便可根据求出的多项式
求出边界数目最大与最小坐标.从拟合曲线可看出矿物的边界分布情形及矿物分布的优选方位.值得注意的是, 曲线的拟合次数在该软件中是低次拟合效果较好(3~5次).反向轮法应变测量中, 判断拟合次数好的标准是, 拟合曲线接近于正弦状, 曲线最高与最低(波峰与波谷)点角度相差近于90°, 这是因为数据点理论上服从椭圆公式.1.3 应变椭圆的绘制
从反向轮法原理中, 可知应变椭圆主轴比是
, 长轴方向是φi=φ′min, 所以在曲线拟合步骤中求出主轴比e和长轴方向ω, 下一步便是怎样把它显示到计算机的屏幕上的问题.椭圆的一般参数方程为:但是该方法只能绘出水平放置的椭圆, 不能满足实际需要.作者的解决方法是: 椭圆上的每一点都绕原点反时针旋转ω角度, 这样才能得到实际需要的应变椭圆.实现旋转的方法如下:
按(xi, yi)坐标点连线绘出的椭圆就完全符合实际应变椭圆形状.按上述设计步骤采用可视化编程工具Visual Basic 5.0开发的CSD软件运行良好, 其功能菜单及程序结构如图 5所示.
图 6 砂岩矿物颗粒轮廓图(a)及手工拟合曲线(b)[3]Fig. 6. Outlines of grain boundaries in a quartzite and its curve fitted by hand2. 检验
该软件适用于矿物颗粒的轮廓图和岩石薄片整体图像, 但整体图像要进行预处理.由于图像上单个矿物颗粒的颜色一般不均一, 所以要在图像处理软件如Adobe PhotoShop等中, 把矿物颗粒颜色处理均匀, 并把解理去掉.作者利用Panozoo的砂岩矿物颗粒轮廓图[3]中的一个矩形区域和其颗粒充添颜色后的图像(图 7)来检验该CSD软件的正确性.
图 6是Panozoo的砂岩矿物颗粒轮廓图[3]及其手工拟合曲线[3].对比其结果和CSD软件得出的结果(如图 8及图 9)显示了由CSD软件得到的应变椭圆, 两者十分相近, 这可以证明CSD软件的正确性和有效性.
3. 结论
CSD软件可实现有限应变的自动测量, 同时可得出矿物分布优选方位且这些结果可以保存.CSD软件实现了反向轮法原理和最小二乘法对数据点进行多项式曲线拟合, 以及椭圆图形绕原点的旋转.这些功能的实现是地质和计算机图形学相结合而完成的.该软件能够增强反向轮法的实用性和扩大反向轮法的应用范围, 减少岩石有限应变测量中人为因素的干扰, 提高工作效率.
-
图 6 砂岩矿物颗粒轮廓图(a)及手工拟合曲线(b)[3]
Fig. 6. Outlines of grain boundaries in a quartzite and its curve fitted by hand
表 1 CSD软件统计出的边界数据和方向
Table 1. The boundary data and orientation collected by the software CSD
-
[1] Ghaleb A R, Fry N. CSTRAIN: A Fortran77 program to study Fry' s plots in two-dimensional simulated models[ J]. Computers & Geosciences, 1995, 21(7): 825~831. [2] Allard B, Benn K. Shape preferred-orientaition analysis using digitized images on a microcomputer[ J]. Computers & Geosciences, 1989, 15(3): 441~448. [3] Panozzo R. Two-dimensional strain determination by the inverse SURFOR wheel[ J]. Journal of Structural Geology, 1987, 9(1): 115~119. doi: 10.1016/0191-8141(87)90049-6 [4] GoodchildJ S, Fueten F. Edge detection in petrographic images using the rotating polarizer stage[ J]. Computers & Geosciences, 1998, 24(8): 745~751. [5] Knappertsbusch M W. A simple Fortran77 program for outline detection[ J]. Computers & Geosciences, 1998, 24(9): 897~900. [6] 刘振安, 苏仕华. C语言图形设计[M]. 北京: 人民出版社, 1995.16~18. [7] Hart D, Rudman A J. Least-squres fit of an ellipse to anisotropic polar data: application to azimuthal resistivity surveys in karst regions[ J]. Computers & Geosciences, 1997, 23(2): 189~194. [8] Thompson G T, Balch S J. An efficient algorithm for polynomial curve fitting[ J]. Computers & Geosciences, 1988, 14(5): 547~556. -