劉世金, 劉大利
(湖北國土資源職業學院,湖北 武漢 430090)
數學地質(mathematical geology)是由數學、地質學和計算機科學結合而成的一門新興邊緣地質科學[1]。它利用數學和計算機科學的理論和方法,定量研究和解決地質科學問題。自20世紀50年代末期開始的計量運動以來,數學地質方法已在礦產資源定量預測、地質多元統計、環境地質評價等民生地質工程領域得到了廣泛的應用[2]。
地質學的定量化離不開數學地質和電子計算[3]。數學地質中的統計分析方法,是現代地質學發展史上計量革命的主要成果之一,具體包括相關分析方法、回歸分析方法、方差分析方法、時間序列分析方法、主成分分析方法、聚類分析方法、趨勢面分析方法和馬爾可夫預測方法等。
“民生地質”是指依靠地質科學的技術、方法、手段來研究和開展地質環境保護、地質安全保障和地質資源利用,以保障人民生產生活,促進經濟社會發展。民生地質涉及地質災害防治、區域地質調查、城市地質、農業地質、水文地質、環境地質等具體地質工程。民生數學地質就是數學地質在民生地質工程中的應用,其研究包括礦產資源定量預測(如蒙特卡羅模擬、Weng旋回模型法等)、地質多元統計分析、地質過程數學模擬等數學地質內容,涉及到的數學方法除地質統計分析方法外,還有小波分析方法、模糊數學方法等等[4]。
數學地質方法是研究和解決民生地質工程的主要方法之一,然而其繁瑣的數學理論推導與復雜的計算機編程往往讓有些地質工作者望而卻步。MATLAB是當今國際公認的最優秀科技應用軟件,具有強大的數值、符號和矩陣計算能力,同時它將計算、可視化和編程等功能集于一易于掌握和使用的環境[5]。本文以數學地質方法中的趨勢面分析為例,調用MATLAB軟件提取民生地質工程中相關有用信息。


圖1 趨勢面示意圖Fig.1 Schematic diagram of trend surface
在數學地質中,用于趨勢面分析的常用數學模型是多項式和傅里葉級數(本文以低維多項式趨勢面分析為例)。
多項式趨勢面的一般形式為:
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,從而確定趨勢面方程。
根據最小二乘法原理,進行多項式趨勢面分析的難點就是求解正規方程組式(6)及其可視化分析。因為其中涉及到矩陣求逆等復雜的數學運算,當階數較高時計算量較大。而這一短板,恰恰是MATLAB矩陣計算和繪圖功能的強項所在。MATLAB是一款集數值分析、矩陣計算和圖形圖像處理等強大的可視化功能于一體的優秀科學和工程應用軟件,而且計算精確、使用簡捷。
表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)中各趨勢面不難看出,在本例降水量空間地理分布規律分析中三維趨勢面分析結果發生了畸變,出現較大的失真。因此,從本例分析結果來看,在有關降水量區域分布規律等民生地質數學方法分析中,應用二次趨勢面進行擬合比較合理。
本文應用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.