劉世金, 劉大利

(湖北國土資源職業(yè)學(xué)院,湖北 武漢　430090)

1　數(shù)學(xué)地質(zhì)與民生地質(zhì)工程

1.1　基本概念

數(shù)學(xué)地質(zhì)(mathematical geology)是由數(shù)學(xué)、地質(zhì)學(xué)和計算機(jī)科學(xué)結(jié)合而成的一門新興邊緣地質(zhì)科學(xué)[1]。它利用數(shù)學(xué)和計算機(jī)科學(xué)的理論和方法,定量研究和解決地質(zhì)科學(xué)問題。自20世紀(jì)50年代末期開始的計量運動以來,數(shù)學(xué)地質(zhì)方法已在礦產(chǎn)資源定量預(yù)測、地質(zhì)多元統(tǒng)計、環(huán)境地質(zhì)評價等民生地質(zhì)工程領(lǐng)域得到了廣泛的應(yīng)用[2]。

地質(zhì)學(xué)的定量化離不開數(shù)學(xué)地質(zhì)和電子計算[3]。數(shù)學(xué)地質(zhì)中的統(tǒng)計分析方法,是現(xiàn)代地質(zhì)學(xué)發(fā)展史上計量革命的主要成果之一,具體包括相關(guān)分析方法、回歸分析方法、方差分析方法、時間序列分析方法、主成分分析方法、聚類分析方法、趨勢面分析方法和馬爾可夫預(yù)測方法等。

1.2　民生數(shù)學(xué)地質(zhì)研究內(nèi)容與方法

“民生地質(zhì)”是指依靠地質(zhì)科學(xué)的技術(shù)、方法、手段來研究和開展地質(zhì)環(huán)境保護(hù)、地質(zhì)安全保障和地質(zhì)資源利用,以保障人民生產(chǎn)生活,促進(jìn)經(jīng)濟(jì)社會發(fā)展。民生地質(zhì)涉及地質(zhì)災(zāi)害防治、區(qū)域地質(zhì)調(diào)查、城市地質(zhì)、農(nóng)業(yè)地質(zhì)、水文地質(zhì)、環(huán)境地質(zhì)等具體地質(zhì)工程。民生數(shù)學(xué)地質(zhì)就是數(shù)學(xué)地質(zhì)在民生地質(zhì)工程中的應(yīng)用,其研究包括礦產(chǎn)資源定量預(yù)測(如蒙特卡羅模擬、Weng旋回模型法等)、地質(zhì)多元統(tǒng)計分析、地質(zhì)過程數(shù)學(xué)模擬等數(shù)學(xué)地質(zhì)內(nèi)容,涉及到的數(shù)學(xué)方法除地質(zhì)統(tǒng)計分析方法外,還有小波分析方法、模糊數(shù)學(xué)方法等等[4]。

數(shù)學(xué)地質(zhì)方法是研究和解決民生地質(zhì)工程的主要方法之一,然而其繁瑣的數(shù)學(xué)理論推導(dǎo)與復(fù)雜的計算機(jī)編程往往讓有些地質(zhì)工作者望而卻步。MATLAB是當(dāng)今國際公認(rèn)的最優(yōu)秀科技應(yīng)用軟件,具有強(qiáng)大的數(shù)值、符號和矩陣計算能力,同時它將計算、可視化和編程等功能集于一易于掌握和使用的環(huán)境[5]。本文以數(shù)學(xué)地質(zhì)方法中的趨勢面分析為例,調(diào)用MATLAB軟件提取民生地質(zhì)工程中相關(guān)有用信息。

2　基于MATLAB的趨勢面分析

2.1　概念

圖1　趨勢面示意圖Fig.1　Schematic diagram of trend surface

2.2　原理

在數(shù)學(xué)地質(zhì)中,用于趨勢面分析的常用數(shù)學(xué)模型是多項式和傅里葉級數(shù)(本文以低維多項式趨勢面分析為例)。

多項式趨勢面的一般形式為:

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

(1)

式中：z為地質(zhì)變量；x、y為觀測點的地理坐標(biāo)。確定上述多項式曲面，即根據(jù)有限個觀測值Mi(xi,yi,zi)(i=1,2,…,n)確定多項式中的各系數(shù)。設(shè)a1,a2,a3,…的估計值為b1,b2,b3,…,則可得近似趨勢面方程為：

(2)

將各觀測值Mi(xi,yi,zi)代入上式后,即可得關(guān)于估計值b1,b2,b3,…的線性方程組：

(3)

實現(xiàn)原理:通過回歸分析原理,運用最小二乘法擬合一個二維非線性函數(shù),模擬地質(zhì)變量在空間上的分布規(guī)律,展示其在地域空間上的變化趨勢。

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

(4)

而Q是關(guān)于b1,b2,b3,…的二次函數(shù),且Q?0,故有：

(5)

即解正規(guī)方程組：

A·B=C

(6)

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

解上述p階正規(guī)方程組,可得b1,b2,…bp,從而確定趨勢面方程。

2.3　方法

根據(jù)最小二乘法原理,進(jìn)行多項式趨勢面分析的難點就是求解正規(guī)方程組式(6)及其可視化分析。因為其中涉及到矩陣求逆等復(fù)雜的數(shù)學(xué)運算,當(dāng)階數(shù)較高時計算量較大。而這一短板,恰恰是MATLAB矩陣計算和繪圖功能的強(qiáng)項所在。MATLAB是一款集數(shù)值分析、矩陣計算和圖形圖像處理等強(qiáng)大的可視化功能于一體的優(yōu)秀科學(xué)和工程應(yīng)用軟件,而且計算精確、使用簡捷。

2.4　應(yīng)用實例

表1為某地區(qū)某月份12個氣象站的平均降水量信息和觀測站地理坐標(biāo)位置數(shù)據(jù)。為揭示降水量的空間地理分布規(guī)律,下面以降水量為因變量z,地理位置的橫、縱坐標(biāo)分別為自變量x、y,將降水量在地理空間位置上的分布情況做趨勢面分析。

表1　某地區(qū)降水量及觀測點地理位置數(shù)據(jù)Table 1　Precipitation and geographical locating data of observation points

根據(jù)最小二乘法原理,將表1中數(shù)據(jù)代入正規(guī)方程組式(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)

應(yīng)用MATLAB進(jìn)行低維趨勢面分析時,可根據(jù)上面二、三次趨勢面方程調(diào)用MATLAB庫函數(shù)surf和contour,分別繪制出二、三次趨勢面及其對應(yīng)的等值線圖，如圖2、圖3所示。

圖2　二次趨勢面及其等值線圖Fig.2　Quadratic trend surface and contour map

圖3　三次趨勢面及其等值線圖Fig.3　Cubic trend surface and contour map

對比二、三次趨勢面等值線圖可知,本例中二次趨勢面的回歸方程較三次趨勢面顯著。

為了進(jìn)一步對比本例中的二次趨勢面與三次趨勢面擬合程度,下面分別調(diào)用plot3和griddata繪制出各觀測值及其插值后的空間分布，分別如圖4-(a)和圖4-(b)所示。

圖4　各觀測值及其插值后的空間分布圖Fig.4　Observation values and spatial distribution map

由圖4-(a)可見,各觀測值在空間的分布表現(xiàn)為12個離散點,從中無法分析降水量的變化趨勢情況。為此,選用MATLAB庫函數(shù)griddata中的′V4′插補(bǔ)方法對觀測值進(jìn)行插值處理,插值后數(shù)據(jù)變化趨勢明顯增強(qiáng),如圖4-(b)所示。同時,對比圖2、圖3和圖4-(b)中各趨勢面不難看出,在本例降水量空間地理分布規(guī)律分析中三維趨勢面分析結(jié)果發(fā)生了畸變,出現(xiàn)較大的失真。因此,從本例分析結(jié)果來看,在有關(guān)降水量區(qū)域分布規(guī)律等民生地質(zhì)數(shù)學(xué)方法分析中,應(yīng)用二次趨勢面進(jìn)行擬合比較合理。

3　結(jié)論

本文應(yīng)用MATLAB科學(xué)、工程軟件,對數(shù)學(xué)地質(zhì)方法在民生地質(zhì)工程中的應(yīng)用進(jìn)行了一些有益的探索。結(jié)果表明，在民生地質(zhì)工程的數(shù)學(xué)地質(zhì)方法分析中,通過運用MATLAB強(qiáng)大的數(shù)值計算及其可視化功能,不僅可以提取到民生地質(zhì)定量分析中的有用信息,而且有利于簡化數(shù)學(xué)地質(zhì)分析過程和總結(jié)地質(zhì)變量變化規(guī)律。研究發(fā)現(xiàn)，在有關(guān)降水量的空間分布等具體民生地質(zhì)工程問題中，利用MATLAB進(jìn)行趨勢面可視化分析時，較宜選取二次趨勢面進(jìn)行趨勢分析,而三次趨勢面分析結(jié)果則往往出現(xiàn)畸變和失真。

參考文獻(xiàn)：

[1]尹大慶,林東維,林景曄.石油數(shù)學(xué)地質(zhì)統(tǒng)計方法及推廣應(yīng)用[J].油氣田地面工程,2013(11):8-9.

[2]湯朝陽,段其發(fā)，等.數(shù)學(xué)地質(zhì)在鄂西南環(huán)境地質(zhì)評價中的應(yīng)用[J].資源環(huán)境與工程,2006，20(2):151-155.

[3]許開達(dá),姚鴻遠(yuǎn),張革勝.數(shù)學(xué)地質(zhì)與計算機(jī)在地質(zhì)方面的應(yīng)用[J].鋼鐵設(shè)計,1993(1):10-17.

[4]徐建華.現(xiàn)代地理學(xué)中的數(shù)學(xué)方法[M].北京:高等教育出版社,2002:37-98.

[5]劉汝敏,曾少軍，等.MATLAB在低維趨勢面分析中的應(yīng)用初探[J].工程地質(zhì)計算機(jī)應(yīng)用,2011(3):2-5.