matlab在地質中的應用

Matlab在地質學中的應用西北大學地質學系2023/2/4主講者：吳鴻天目錄1二維三維圖形的繪制2字符數組與符號運算3簡單編程的數據處理Matlab在地質中的應用舉例1Matlab概述2矩陣的建立與計算3Matlab的簡單操作Matlab的簡要介紹12Matlab的簡要介紹Part11.1Matlab概述MATLAB代表

MATrix

LABoratory1.1Matlab概述Matlab的首創(chuàng)者是美國新墨西哥大學計算機系的系主任CleveMoler博士，他在教授線性代數課程發(fā)現(xiàn)其他語言很不方便，便構思開發(fā)了MATLAB。最初采用FORTRAN語言編寫，20世紀80年代后出現(xiàn)了MATLAB的第二版，全部采用C語言編寫。

1984年Moler博士和一批數學家及軟件專家創(chuàng)建了MathWorks公司，專門開發(fā)MATLAB。1993年出現(xiàn)了微機版，到現(xiàn)在最新版本是2013b。1.1Matlab概述

Matlab是美國MathWorks公司出品的商業(yè)數學軟件，主要包括Matlab和Simulink兩大部分，在數學類科技應用軟件中在數值計算方面首屈一指。

Matlab可以進行矩陣運算、繪制函數和數據、實現(xiàn)算法、創(chuàng)建用戶界面、連接其他編程語言的程序等，主要應用于工程計算、控制設計、信號處理與通訊、圖像處理、信號檢測、金融建模設計與分析等領域。1.1Matlab概述

主要功能：數值分析數值和符號計算工程與科學繪圖控制系統(tǒng)設計與仿真數字圖像處理數字信號處理通訊系統(tǒng)設計與仿真財務與金融工程1.1Matlab概述

顯著特點：

1Matlab具有強大的矩陣運算能力，MatrixLaboratory（矩陣實驗室），使得矩陣運算非常簡單。

2它是一種演算式語言，包含大量計算算法，封裝了大量的數學函數，極大地減少了用戶編寫復雜函數的難度，而且數學表達式和運算規(guī)則與通常的習慣相同。3Matlab對許多專門的領域都開發(fā)了功能強大的模塊集和工具箱用戶可以直接使用工具箱學習、應用和評估不同的方法而不需要自己編寫代碼。1.1Matlab概述Matlab界面介紹1.1Matlab概述命令窗口

命令窗口是命令行語句和命令文件執(zhí)行的主要窗口。

在命令窗口中直接輸入命令或MATLAB函數，系統(tǒng)自動反饋結果。1.1Matlab概述編輯窗口

用來創(chuàng)建和修改M-files(MATLAB腳本)路徑是File—New—M-file或Openfile...，編輯完成后按F5運行程序。

1.1Matlab概述圖形窗口

圖形窗的打開方式有三種：1在主窗口上依次選擇File→New→Figure或Openfile…操作時，選擇打開圖形文件。2在命令窗口輸入figure命令。3執(zhí)行結果為圖形的語句。1.1Matlab概述圖形窗口如執(zhí)行以下語句x=0:0.1:6;y=sin(x);plot(x,y);在Matlab中可得到右側圖形1.1Matlab概述Matlab語言的基本要素

Matlab語言主要是由幾個基本要素構建而成的，其中包括變量、常量、數值、字符串、運算符及標點等。變量的命名規(guī)則：變量名區(qū)分大小寫；變量名長度不超過31位，第31位之后的字符被忽略；變量名以字母開頭，變量名中可以包含字母、數字、下劃線，但不能使用標點。

1.2矩陣的建立與計算矩陣是Matlab語言的核心由m×n個數（i=1，2，...，m；j=1，2，...，n）排成的m行n列的數表稱為m行n列矩陣，簡稱m×n矩陣。為表示它是一個整體，總是加一個括弧。線性方程組的解取決于系數常數項1.2矩陣的建立與計算對線性方程組的研究可轉化為對這張表的研究。線性方程組的系數與常數項按原位置可排為1.2矩陣的建立與計算1.2矩陣的建立與計算

由空間基本定理知，有且只有一組實數(x,y,z)，使得a=向量OP=xi+yj+zk，因此把實數對(x,y,z)叫做向量a的坐標，記作a=(x,y,z)。

某廠向三個商店發(fā)送四種產品，其發(fā)送的數量可用矩陣來刻劃。若用aij

表示為工廠向第i店發(fā)送第j種產品數量，則矩陣

表示了工廠向三個商店發(fā)送四種產品的數量。1.2矩陣的建立與計算1.2矩陣的建立與計算矩陣的生成注意點：輸入矩陣時要以“[]”為其標識，即矩陣的元素應在“[]”內部。矩陣的同行元素之間可由空格或“，”分隔，行與行之間用“;”或回車符分隔。矩陣元素可為運算表達式。無任何元素的空矩陣也合法。1.2矩陣的建立與計算>>A=[1,2,3;4,5,6;7,8,9]A=123456789創(chuàng)建帶有運算表達式的矩陣：>>b=[sin(pi/3),cos(pi/4);log(9),tan(pi/4)]>>b=[sin(pi/3),cos(pi/4);log(9),tan(pi/4)]b=0.86600.70712.19721.00001.2矩陣的建立與計算矩陣的計算

在matlab語言中運算符包括算術運算符、關系運算符和邏輯運算符。其中算術運算符見右表：

其中點運算是指操作元素點對點的運算，也就是說矩陣內元素對元素之間的運算，點運算要求參與運算的變量在結構上必須是相似的。操作符定義+算術加-算術減*算術乘.*點乘^算術乘方.^點乘方\算術左除.\點左除/算術右除./點右除1.2矩陣的建立與計算矩陣的計算矩陣的計算包括矩陣的加法、減法、乘法、轉置等。A

=123

456780

B=147258360C=A+B返回:C=261061014101401.2矩陣的建立與計算矩陣的計算

矩陣乘法用“*”符號表示，只有當A矩陣列數與B矩陣的行數相等時，二者才可以進行乘法運算，否則就是錯誤的。計算方法和線性代數中所介紹的完全相同。

如：A=[1,2;3,4]；B=[5,6;7,8]；C=A*B，結果為C=12×56=1×5+2×71×6+2×8=1922

34783×5+4×73×6+4×843501.2矩陣的建立與計算矩陣的計算

數組的乘（矩陣的點乘）用符號.*表示，如果A與B矩陣具有相同階數，則A.*B表示A和B單個元素之間的對應相乘。1.2矩陣的建立與計算>>a=[1,2,3;4,5,6;7,8,9]a=123456789>>b=[123456789]b=123456789>>c=a*bc=303642668196102126150>>d=a.*bd=1491625364964811.2矩陣的建立與計算如何畫出sin(x)xcos(x)函數在0-pi區(qū)間上的曲線？x=0:0.1:pi;y=sin(x).*cos(x);plot(x,y)提示：通過點運算獲得。1.2矩陣的建立與計算

字符串：與矩陣基本等價字符串用單引號輸入或賦值；字符串的每個字符都是都是字符數組的一個元素；字符串的每個字符（包括空格）都是字符數組的一個元素。>>s='iloveyou's=iloveyou>>size(s)ans=110>>s(2)ans=>>s(10)ans=usize():獲取矩陣的行數和列數1.3Matlab的簡單操作標點：在matlab中一些標點符號也被賦以特殊的意義，或表示要進行一定的運算等。標點定義標點定義：冒號，具有多種應用功能…續(xù)行符；區(qū)分行，及取消結果顯示%注釋標記，區(qū)分列，及函數參數分隔符等!調用操作系統(tǒng)運算()指定運算過程中的先后次序等[]矩陣定義的標志等{}用于構成單元數組等‘字符串的標示符1.3Matlab的簡單操作下面對上表中各標點作進一步的介紹。

首先，用冒號來定義行向量?；靖袷剑簒＝x1:step:x2（初始值：步長：尾元素數值限）

x＝x1:x2（默認步長為1）其次，通過使用冒號，可以截取指定矩陣中的部分。>>a=[1:18]；a=reshape(a,3,6)%重組矩陣a=147101316258111417369121518>>a1=a(:,1,3)

%提取矩陣a的第1到3列元素組成的子矩陣。

a1=147258369>>a2=a(1:2,3:end)a2=710131681114171.3Matlab的簡單操作

分號在matlab語言中所起的作用是作為分隔行的標志，這里包括在矩陣定義中對行的區(qū)分以及在程序代碼中對行的劃分。在矩陣定義中，分號之間的數據被認為是矩陣的同一行。對于以分號結尾的行語句，matlab語言不會把其運算結果顯示在命令窗口中。

矩陣的定義需要使用中括號。此外，中括號可以表示空矩陣，并通過該空矩陣實現(xiàn)對矩陣的部分刪除。1.3Matlab的簡單操作1.3Matlab的簡單操作>>a=[1,2,3,4;5,6,7,8;4,3,2,1;8,7,6,5]a=1234567843218765>>a(:,1:2)=[]a=34782165>>b=a([13],:)b=34211.3Matlab的簡單操作編程語句結構Matlab語言最基本的賦值語句結構為：變量名列表=表達式注1：整個賦值語句以；結束，則不在屏幕上返回

結果，否則立即返回結果。注2：多個語句可在同一行，用逗號分開。注3：表達是太長可以用續(xù)行符號…1.3Matlab的簡單操作編程語句結構循環(huán)語句for,while條件轉移ifend,ifelseifelseend開關語句switchcase注釋語句%中斷語句break暫停語句pause回顯語句echoon/off1.3Matlab的簡單操作for循環(huán)語句

for循環(huán)語句其功能是重復執(zhí)行一條或一組語句的固定次數(也可以是零次)。

表達式實際上是一個矩陣，因為在Matlab中基本元素就是矩陣。當矩陣的每一列元素依次賦給變量時，循環(huán)體語句就執(zhí)行一次。1.3Matlab的簡單操作for循環(huán)語句基本格式：for循環(huán)變量＝起始值：步長：終止值循環(huán)體end

步長缺省值為1，可以在正實數或負實數范圍內任意指定。對于正數，循環(huán)變量的值大于終止值時，循環(huán)結束；對于負數，循環(huán)變量的值小于終止值時，循環(huán)結束。循環(huán)結構可以嵌套使用。1.3Matlab的簡單操作for循環(huán)語句fori=1:5forj=1:5a(i,j)=1/(i+j-1);endenda=1.00000.50000.33330.25000.20000.50000.33330.25000.20000.16670.33330.25000.20000.16670.14290.25000.20000.16670.14290.12500.20000.16670.14290.12500.11111.3Matlab的簡單操作while循環(huán)語句基本格式：while表達式循環(huán)體end

若表達式為真，則執(zhí)行循環(huán)體的內容，執(zhí)行后再判斷表達式是否為真，若不為真，則跳出循環(huán)體，向下繼續(xù)執(zhí)行。

其功能是重復執(zhí)行一條或者一組語句，重復次數由某個邏輯條件控制，—般來說其執(zhí)行次數是不能預先確定的。1.3Matlab的簡單操作while循環(huán)語句x=0;sum=0;whilex<101sum=sum+x;x=x+2;endsum

While循環(huán)和for循環(huán)的區(qū)別在于，while循環(huán)結構的循環(huán)體被執(zhí)行的次數不是確定的，而for結構中循環(huán)體的執(zhí)行次數是確定的。1.3Matlab的簡單操作if分支結構（1）if邏輯表達式執(zhí)行語句

end

當邏輯表達式的值為真時，執(zhí)行該結構中的執(zhí)行語句，執(zhí)行完之后繼續(xù)向下進行；若為假，則跳過結構中的內容，向下執(zhí)行。（2）if邏輯表達式執(zhí)行語句1else執(zhí)行語句2end

if-else的執(zhí)行方式為：如果邏輯表達式的值為真，則執(zhí)行語句1，然后跳過語句2，向下執(zhí)行；如果為假，則執(zhí)行語句2，然后向下執(zhí)行。1.3Matlab的簡單操作if分支結構

（3）if邏輯表達式1執(zhí)行語句1elseif邏輯表達式2

執(zhí)行語句2

…end

if-elseif的執(zhí)行方式為：如果邏輯表達式1的值為真，則執(zhí)行語句1；如果為假，則判斷邏輯表達式2，如果為真，則執(zhí)行語句2，否則向下執(zhí)行。1.3Matlab的簡單操作if分支結構n=input('Pleaseenteranumber,n=');ifn<10n+10elsen-10end例：輸入一個數，如果小于10就輸出這個數加10；否則輸出這個數減10。1.3Matlab的簡單操作switch分支結構switch表達式（％可以是標量或字符串）case值1語句1case值2語句2….otherwise語句3end

執(zhí)行方式：表達式的值和哪種情況(case)的值相同，就執(zhí)行哪種情況中的語句，如果不同，則執(zhí)行otherwise中的語句。格式中也可以不包括otherwise，這時如果表達式的值與列出的各種情況都不相同，則繼續(xù)向下執(zhí)行。Matlab在地質中的應用舉例Part22.1二維三維圖形的繪制二維圖形及相關函數

繪制二維圖形最常用的函數是plot函數。其調用格式為：

plot(X,Y)

X為橫坐標，Y為縱坐標。

說明：x,y可以是向量或矩陣。

當x，y均為向量時，要求向量x與向量y的長度一致，則plot(x,y)繪制出以x為橫坐標，y為縱坐標的二維圖形。

2.1二維三維圖形的繪制二維圖形及相關函數例如：x，y為同維向量時

x=0:0.01*pi:pi;y=sin(x).*cos(x);plot(x,y)2.1二維三維圖形的繪制二維圖形及相關函數

當x為向量、y為矩陣時，plot(x,y)用不同顏色的曲線繪制出y行或列對于x的圖形。y矩陣的行或列的選擇取決于x，y的維數，若y為方陣或y矩陣的列向量長度與x向量的長度一致，則繪制出y矩陣的各個列向量相對于x的一組二維圖形；若y矩陣的行向量長度與x向量的長度一致，則繪制出y矩陣的各個行向量相對于x的一組二維圖形。2.1二維三維圖形的繪制二維圖形及相關函數clearallx=0:pi/50:3*pi;y(1,:)=cos(x);y(2,:)=0.3*cos(x);y(3,:)=0.6*cos(x);plot(x,y)2.1二維三維圖形的繪制二維圖形及相關函數

若x,y是同維的矩陣，則plot(x,y)繪制出y的列向量相對于x的列向量的一組二維圖形。例：x=[465;132];y=[798;121110];plot(x,y)2.1二維三維圖形的繪制二維圖形及相關函數plot(x,y,'參數')說明：x,y可以是向量或矩陣，與前面相同。參數選項為一個字符串，它決定了二維圖形的顏色、線型以及數據點的圖標等。2.1二維三維圖形的繪制選項說明選項說明選項說明-實線g綠色xx符號：點線b藍色++號-.點劃線w白色*星號--虛線k黑色s方形y黃色.點d菱形r紅色o圓p正五邊形常用曲線顏色和線型格式2.1二維三維圖形的繪制二維圖形及相關函數例如：x=1:0.1*pi:2*pi;y=sin(x);z=cos(x);plot(x,y,'--k',x,z,'-.rd')%結果怎樣？2.1二維三維圖形的繪制極坐標系繪圖函數polar

polar(theta,tho,s)

其中theta為弧度表示的角度向量；tho為相應的幅向量；s為圖形屬性設置選項，可省略。例如：x=0:0.01*pi:4*pi;y=sin(x/2)+x;polar(x,y,'-')2.1二維三維圖形的繪制二維圖形及相關函數：繪制玫瑰花圖

玫瑰花圖實際上是極坐標中的直方圖，它根據數據值的范圍進行分組，并顯示數據的分布特征。每個組顯示為一瓣。

用rose函數繪制玫瑰花圖，該函數最常用的語法格式為rose(theta)

用theta數據繪制玫瑰花圖，玫瑰花的瓣數一般為20，或更少。向量theta中的數據用弧度表示，確定每個花瓣的角度?；ò甑拈L度反映了theta中落在該范圍內的元素個數。2.1二維三維圖形的繪制二維圖形及相關函數：繪制玫瑰花圖例：用玫瑰花圖顯示50個隨機數的分布特征。clearalltheta=2*pi*rand(1,50);%rand為生成均勻分布的隨機數函數rose(theta)2.1二維三維圖形的繪制三維圖形繪制1plot3函數

最基本的三維圖形函數為plot3，它是將二維函數plot的有關功能擴展到三維空間，用來繪制三維圖形。

函數格式：plot3(x1,y1,z1,c1,x2,y2,z2,c2,…)

其中x1,y1,z1…表示三維坐標向量，c1，c2…表示線形或顏色。

函數功能：以向量x，y，z為坐標，繪制三維曲線。

2.1二維三維圖形的繪制三維圖形繪制例如：x=0:pi/50:10*pi;y=sin(x);z=cos(x);plot3(x,y,z);2.1二維三維圖形的繪制三維圖形繪制

plot3也可以以矩陣作為輸入參數，這時要求3個輸入矩陣必須結構相同。例如：[x,y]=meshgrid(-2:0.1:2,-2:0.1:2);z=x.*exp(-x.^2-y.^2);plot3(x,y,z)meshgird函數：按采樣點生成繪制3D圖

形所需的網格數據2.1二維三維圖形的繪制三維圖形繪制

2

mesh函數

mesh函數用于繪制三維網格圖。在不需要繪制特別精細的三維曲面結構圖時，可以通過繪制三維網格圖來表示三維曲面。三維曲面的網格圖最突出的優(yōu)點是：它較好地解決了實驗數據在三維空間的可視化問題。

函數格式：mesh(x,y,z,c)

其中x，y控制X和Y軸坐標，矩陣z是由(x，y)求得Z軸坐標，(x,y,z組成了三維空間的網格點；c用于控制網格點顏色。

meshc(x,y,z)繪制將網格與等高線結合的三維面圖。

meshz(x,y,z)繪制以邊界面屏蔽的網格圖

surf(x,y,z)繪制三維著色面圖2.1二維三維圖形的繪制三維圖形繪制例如：[X,Y]=meshgrid([-4:0.5:4]);z=sqrt(X.^2+Y.^2);subplot(2,2,1)mesh(z)subplot(2,2,2)meshc(z)subplot(2,2,3)meshz(z)subplot(2,2,4)surf(z)sqrt()求平方根，subplot（）是將多個圖畫到一個平面上的工具。2.1二維三維圖形的繪制Matlab在構造地質學中的應用選擇構造地質學課本的附圖1：凌河地形地質圖來進行試驗

在構造地質學課程的學習中，我們要接觸大量的地形地質圖，地形地質圖上最基本的一項信息就是等高線地形圖，讀懂等高線地形圖，通過空間想象將三維的立體地形在腦海中構建起來，這對于接下來正確獲取其它地學專業(yè)信息有極大的幫助。但是，在剛剛接觸地質圖的時候，我們通常很難建立起三維的立體地形，或者是建立起來了也沒有一個客觀的對照標準，不知道自己建立的影像到底對不對，老師在講解的時候，也只是簡單的語言描述，我們根本沒法建立清晰可見的三維立體地形對照影像，不利于我們后續(xù)課程的學習。2.1二維三維圖形的繪制2.1二維三維圖形的繪制

為了獲取這一地區(qū)的點的高度，在原圖上，以左下角為坐標原點，圖的下邊界為X軸，圖的左邊界為Y軸，分別在兩軸上每隔1cm取點，根據圖的比例尺1:20000，圖上的1cm在實際中表示200米，根據取的點作方格，如圖2，再根據地質圖上原有的等高線的信息，就可以得到一系列的（X,Y,H）。其中（X,Y）表示在這張地形地質圖上某點的位置坐標，H表示該點的高程（通過圖上的等高線可以得出來）。下表就是得到的數據：Matlab在構造地質學中的應用2.1二維三維圖形的繪制在編程中為了導入的方便，將相同的數據改用txt格式進行導入(h.txt)。2.1二維三維圖形的繪制

用原始數據作該地區(qū)的地貌圖。x=0:200:3000;y=0:200:2400;loadh.txt;meshz(x,y,h);xlabel('X'),ylabel('Y'),zlabel('H')運行結果如圖3。明顯山區(qū)表面很粗糙，效果不好。2.1二維三維圖形的繪制

每間隔50加密網格，用最近鄰點插值方法作該地區(qū)的地貌圖，可以使地圖的顯示效果更好。x=0:200:3000;y=0:200:2400;loadh.txt;xi=0:50:3000;yi=0:50:2400;hi=interp2(x,y,h,xi,yi','cubic');surfc(xi,yi,hi);xlabel('X'),ylabel('Y'),zlabel('H')運行結果如圖，效果很好。2.1二維三維圖形的繪制

顯示一個旋轉的地球（1）創(chuàng)建一個包含下面代碼的M文件：fori=360:-5:0view(i,0);drawnowend（2）保存名稱為spin.m。（3）用格線創(chuàng)建球體顯示效果：axesm('globe','Grid','on','Gcolor',[.7,.8,.9],'GlineStyle','-')（4）隱藏圖形包圍盒的邊線，進行透視投影：set(gca,'Box','off','Projection','perspective')（5）創(chuàng)建海平面數據網格，使球體不透明：base=zeros(180,360);baseref=[1900];hs=meshm(base,baseref,size(base));colormapcopper（6）提升網格，使它浮在球面上方2.5%半徑

距離的地方：setm(gca,'Galtitude',0.025)（7）準備用地形地貌數據替換目前的球面數據：clmo(hs)loadtopo2.1二維三維圖形的繪制

顯示一個旋轉的地球（8）將高程放大顯示，繪制表面：topo=topo/(almanac('earth','radius')*20);hs=meshm(topo,topolegend,size(topo),topo)；demcmap(topo)（9）旋轉地球spin2.2字符數組與符號運算在數值計算中，最主要的部分就是解線性方程組。由矩陣的定義可知，線性方程組的運算可以轉化為矩陣的運算。另外，在地質學中應用比較廣泛的還有插值、曲線的擬合等等。趨勢面分析：為研究某一地質特征的區(qū)域性變化規(guī)律和局部范圍的異常變化,常用趨勢面分析方法來處理地質問題。用matlab處理的過程就是進行多個變量的函數的曲線擬合過程。2.2字符數組與符號運算具體過程為：輸入原始數據y建立趨勢方程y計算趨勢值和偏差值y作分析圖和剩余圖。在matlab中：一次趨勢方程為：W=