刘世金, 刘大利

(湖北国土资源职业学院,湖北 武汉　430090)

1　数学地质与民生地质工程

1.1　基本概念

数学地质(mathematical geology)是由数学、地质学和计算机科学结合而成的一门新兴边缘地质科学[1]。它利用数学和计算机科学的理论和方法,定量研究和解决地质科学问题。自20世纪50年代末期开始的计量运动以来,数学地质方法已在矿产资源定量预测、地质多元统计、环境地质评价等民生地质工程领域得到了广泛的应用[2]。

地质学的定量化离不开数学地质和电子计算[3]。数学地质中的统计分析方法,是现代地质学发展史上计量革命的主要成果之一,具体包括相关分析方法、回归分析方法、方差分析方法、时间序列分析方法、主成分分析方法、聚类分析方法、趋势面分析方法和马尔可夫预测方法等。

1.2　民生数学地质研究内容与方法

“民生地质”是指依靠地质科学的技术、方法、手段来研究和开展地质环境保护、地质安全保障和地质资源利用,以保障人民生产生活,促进经济社会发展。民生地质涉及地质灾害防治、区域地质调查、城市地质、农业地质、水文地质、环境地质等具体地质工程。民生数学地质就是数学地质在民生地质工程中的应用,其研究包括矿产资源定量预测(如蒙特卡罗模拟、Weng旋回模型法等)、地质多元统计分析、地质过程数学模拟等数学地质内容,涉及到的数学方法除地质统计分析方法外,还有小波分析方法、模糊数学方法等等[4]。

数学地质方法是研究和解决民生地质工程的主要方法之一,然而其繁琐的数学理论推导与复杂的计算机编程往往让有些地质工作者望而却步。MATLAB是当今国际公认的最优秀科技应用软件,具有强大的数值、符号和矩阵计算能力,同时它将计算、可视化和编程等功能集于一易于掌握和使用的环境[5]。本文以数学地质方法中的趋势面分析为例,调用MATLAB软件提取民生地质工程中相关有用信息。

2　基于MATLAB的趋势面分析

2.1　概念

图1　趋势面示意图Fig.1　Schematic diagram of trend surface

2.2　原理

在数学地质中,用于趋势面分析的常用数学模型是多项式和傅里叶级数(本文以低维多项式趋势面分析为例)。

多项式趋势面的一般形式为:

z=a0+a1x+a2y+a3x2+a4xy+a5y2+…

(1)

式中：z为地质变量；x、y为观测点的地理坐标。确定上述多项式曲面，即根据有限个观测值Mi(xi,yi,zi)(i=1,2,…,n)确定多项式中的各系数。设a1,a2,a3,…的估计值为b1,b2,b3,…,则可得近似趋势面方程为：

(2)

将各观测值Mi(xi,yi,zi)代入上式后,即可得关于估计值b1,b2,b3,…的线性方程组：

(3)

实现原理:通过回归分析原理,运用最小二乘法拟合一个二维非线性函数,模拟地质变量在空间上的分布规律,展示其在地域空间上的变化趋势。

由最小二乘法原理,使得残差平方和为：

(4)

而Q是关于b1,b2,b3,…的二次函数,且Q≻0,故有：

(5)

即解正规方程组：

A·B=C

(6)

其中,A=XTX,B=(b1b2…bp)T,C=XTZ,Z=(z1z2…zn)T,

解上述p阶正规方程组,可得b1,b2,…bp,从而确定趋势面方程。

2.3　方法

根据最小二乘法原理,进行多项式趋势面分析的难点就是求解正规方程组式(6)及其可视化分析。因为其中涉及到矩阵求逆等复杂的数学运算,当阶数较高时计算量较大。而这一短板,恰恰是MATLAB矩阵计算和绘图功能的强项所在。MATLAB是一款集数值分析、矩阵计算和图形图像处理等强大的可视化功能于一体的优秀科学和工程应用软件,而且计算精确、使用简捷。

2.4　应用实例

表1为某地区某月份12个气象站的平均降水量信息和观测站地理坐标位置数据。为揭示降水量的空间地理分布规律,下面以降水量为因变量z,地理位置的横、纵坐标分别为自变量x、y,将降水量在地理空间位置上的分布情况做趋势面分析。

表1　某地区降水量及观测点地理位置数据Table 1　Precipitation and geographical locating data of observation points

根据最小二乘法原理,将表1中数据代入正规方程组式(6),运用MATLAB软件计算得二次、三次趋势面多项式分别如下:

z=5.998+17.438x+29.787y-3.558x2
+0.357xy-8.070y2

(7)

z=-48.810+37.557x+130.130y+8.389x2-
33.166xy-62.740y2-4.133x3+6.138x2y
+2.566xy2+9.785y3

(8)

应用MATLAB进行低维趋势面分析时,可根据上面二、三次趋势面方程调用MATLAB库函数surf和contour,分别绘制出二、三次趋势面及其对应的等值线图，如图2、图3所示。

图2　二次趋势面及其等值线图Fig.2　Quadratic trend surface and contour map

图3　三次趋势面及其等值线图Fig.3　Cubic trend surface and contour map

对比二、三次趋势面等值线图可知,本例中二次趋势面的回归方程较三次趋势面显著。

为了进一步对比本例中的二次趋势面与三次趋势面拟合程度,下面分别调用plot3和griddata绘制出各观测值及其插值后的空间分布，分别如图4-(a)和图4-(b)所示。

图4　各观测值及其插值后的空间分布图Fig.4　Observation values and spatial distribution map

由图4-(a)可见,各观测值在空间的分布表现为12个离散点,从中无法分析降水量的变化趋势情况。为此,选用MATLAB库函数griddata中的′V4′插补方法对观测值进行插值处理,插值后数据变化趋势明显增强,如图4-(b)所示。同时,对比图2、图3和图4-(b)中各趋势面不难看出,在本例降水量空间地理分布规律分析中三维趋势面分析结果发生了畸变,出现较大的失真。因此,从本例分析结果来看,在有关降水量区域分布规律等民生地质数学方法分析中,应用二次趋势面进行拟合比较合理。

3　结论

本文应用MATLAB科学、工程软件,对数学地质方法在民生地质工程中的应用进行了一些有益的探索。结果表明，在民生地质工程的数学地质方法分析中,通过运用MATLAB强大的数值计算及其可视化功能,不仅可以提取到民生地质定量分析中的有用信息,而且有利于简化数学地质分析过程和总结地质变量变化规律。研究发现，在有关降水量的空间分布等具体民生地质工程问题中，利用MATLAB进行趋势面可视化分析时，较宜选取二次趋势面进行趋势分析,而三次趋势面分析结果则往往出现畸变和失真。

参考文献：

[1]尹大庆,林东维,林景晔.石油数学地质统计方法及推广应用[J].油气田地面工程,2013(11):8-9.

[2]汤朝阳,段其发，等.数学地质在鄂西南环境地质评价中的应用[J].资源环境与工程,2006，20(2):151-155.

[3]许开达,姚鸿远,张革胜.数学地质与计算机在地质方面的应用[J].钢铁设计,1993(1):10-17.

[4]徐建华.现代地理学中的数学方法[M].北京:高等教育出版社,2002:37-98.

[5]刘汝敏,曾少军，等.MATLAB在低维趋势面分析中的应用初探[J].工程地质计算机应用,2011(3):2-5.