matlab在地质中的应用

Matlab在地质学中的应用西北大学地质学系2023/2/4主讲者：吴鸿天目录1二维三维图形的绘制2字符数组与符号运算3简单编程的数据处理Matlab在地质中的应用举例1Matlab概述2矩阵的建立与计算3Matlab的简单操作Matlab的简要介绍12Matlab的简要介绍Part11.1Matlab概述MATLAB代表

MATrix

LABoratory1.1Matlab概述Matlab的首创者是美国新墨西哥大学计算机系的系主任CleveMoler博士，他在教授线性代数课程发现其他语言很不方便，便构思开发了MATLAB。最初采用FORTRAN语言编写，20世纪80年代后出现了MATLAB的第二版，全部采用C语言编写。

1984年Moler博士和一批数学家及软件专家创建了MathWorks公司，专门开发MATLAB。1993年出现了微机版，到现在最新版本是2013b。1.1Matlab概述

Matlab是美国MathWorks公司出品的商业数学软件，主要包括Matlab和Simulink两大部分，在数学类科技应用软件中在数值计算方面首屈一指。

Matlab可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。1.1Matlab概述

主要功能：数值分析数值和符号计算工程与科学绘图控制系统设计与仿真数字图像处理数字信号处理通讯系统设计与仿真财务与金融工程1.1Matlab概述

显著特点：

1Matlab具有强大的矩阵运算能力，MatrixLaboratory（矩阵实验室），使得矩阵运算非常简单。

2它是一种演算式语言，包含大量计算算法，封装了大量的数学函数，极大地减少了用户编写复杂函数的难度，而且数学表达式和运算规则与通常的习惯相同。3Matlab对许多专门的领域都开发了功能强大的模块集和工具箱用户可以直接使用工具箱学习、应用和评估不同的方法而不需要自己编写代码。1.1Matlab概述Matlab界面介绍1.1Matlab概述命令窗口

命令窗口是命令行语句和命令文件执行的主要窗口。

在命令窗口中直接输入命令或MATLAB函数，系统自动反馈结果。1.1Matlab概述编辑窗口

用来创建和修改M-files(MATLAB脚本)路径是File—New—M-file或Openfile...，编辑完成后按F5运行程序。

1.1Matlab概述图形窗口

图形窗的打开方式有三种：1在主窗口上依次选择File→New→Figure或Openfile…操作时，选择打开图形文件。2在命令窗口输入figure命令。3执行结果为图形的语句。1.1Matlab概述图形窗口如执行以下语句x=0:0.1:6;y=sin(x);plot(x,y);在Matlab中可得到右侧图形1.1Matlab概述Matlab语言的基本要素

Matlab语言主要是由几个基本要素构建而成的，其中包括变量、常量、数值、字符串、运算符及标点等。变量的命名规则：变量名区分大小写；变量名长度不超过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)。

某厂向三个商店发送四种产品，其发送的数量可用矩阵来刻划。若用aij

表示为工厂向第i店发送第j种产品数量，则矩阵

表示了工厂向三个商店发送四种产品的数量。1.2矩阵的建立与计算1.2矩阵的建立与计算矩阵的生成注意点：输入矩阵时要以“[]”为其标识，即矩阵的元素应在“[]”内部。矩阵的同行元素之间可由空格或“，”分隔，行与行之间用“;”或回车符分隔。矩阵元素可为运算表达式。无任何元素的空矩阵也合法。1.2矩阵的建立与计算>>A=[1,2,3;4,5,6;7,8,9]A=123456789创建带有运算表达式的矩阵：>>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区间上的曲线？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中一些标点符号也被赋以特殊的意义，或表示要进行一定的运算等。标点定义标点定义：冒号，具有多种应用功能…续行符；区分行，及取消结果显示%注释标记，区分列，及函数参数分隔符等!调用操作系统运算()指定运算过程中的先后次序等[]矩阵定义的标志等{}用于构成单元数组等‘字符串的标示符1.3Matlab的简单操作下面对上表中各标点作进一步的介绍。

首先，用冒号来定义行向量。基本格式：x＝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语言中所起的作用是作为分隔行的标志，这里包括在矩阵定义中对行的区分以及在程序代码中对行的划分。在矩阵定义中，分号之间的数据被认为是矩阵的同一行。对于以分号结尾的行语句，matlab语言不会把其运算结果显示在命令窗口中。

矩阵的定义需要使用中括号。此外，中括号可以表示空矩阵，并通过该空矩阵实现对矩阵的部分删除。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：表达是太长可以用续行符号…1.3Matlab的简单操作编程语句结构循环语句for,while条件转移ifend,ifelseifelseend开关语句switchcase注释语句%中断语句break暂停语句pause回显语句echoon/off1.3Matlab的简单操作for循环语句

for循环语句其功能是重复执行一条或一组语句的固定次数(也可以是零次)。

表达式实际上是一个矩阵，因为在Matlab中基本元素就是矩阵。当矩阵的每一列元素依次赋给变量时，循环体语句就执行一次。1.3Matlab的简单操作for循环语句基本格式：for循环变量＝起始值：步长：终止值循环体end

步长缺省值为1，可以在正实数或负实数范围内任意指定。对于正数，循环变量的值大于终止值时，循环结束；对于负数，循环变量的值小于终止值时，循环结束。循环结构可以嵌套使用。1.3Matlab的简单操作for循环语句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循环语句基本格式：while表达式循环体end

若表达式为真，则执行循环体的内容，执行后再判断表达式是否为真，若不为真，则跳出循环体，向下继续执行。

其功能是重复执行一条或者一组语句，重复次数由某个逻辑条件控制，—般来说其执行次数是不能预先确定的。1.3Matlab的简单操作while循环语句x=0;sum=0;whilex<101sum=sum+x;x=x+2;endsum

While循环和for循环的区别在于，while循环结构的循环体被执行的次数不是确定的，而for结构中循环体的执行次数是确定的。1.3Matlab的简单操作if分支结构（1）if逻辑表达式执行语句

end

当逻辑表达式的值为真时，执行该结构中的执行语句，执行完之后继续向下进行；若为假，则跳过结构中的内容，向下执行。（2）if逻辑表达式执行语句1else执行语句2end

if-else的执行方式为：如果逻辑表达式的值为真，则执行语句1，然后跳过语句2，向下执行；如果为假，则执行语句2，然后向下执行。1.3Matlab的简单操作if分支结构

（3）if逻辑表达式1执行语句1elseif逻辑表达式2

执行语句2

…end

if-elseif的执行方式为：如果逻辑表达式1的值为真，则执行语句1；如果为假，则判断逻辑表达式2，如果为真，则执行语句2，否则向下执行。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

执行方式：表达式的值和哪种情况(case)的值相同，就执行哪种情况中的语句，如果不同，则执行otherwise中的语句。格式中也可以不包括otherwise，这时如果表达式的值与列出的各种情况都不相同，则继续向下执行。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中的数据用弧度表示，确定每个花瓣的角度。花瓣的长度反映了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函数用于绘制三维网格图。在不需要绘制特别精细的三维曲面结构图时，可以通过绘制三维网格图来表示三维曲面。三维曲面的网格图最突出的优点是：它较好地解决了实验数据在三维空间的可视化问题。

函数格式：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：凌河地形地质图来进行试验

在构造地质学课程的学习中，我们要接触大量的地形地质图，地形地质图上最基本的一项信息就是等高线地形图，读懂等高线地形图，通过空间想象将三维的立体地形在脑海中构建起来，这对于接下来正确获取其它地学专业信息有极大的帮助。但是，在刚刚接触地质图的时候，我们通常很难建立起三维的立体地形，或者是建立起来了也没有一个客观的对照标准，不知道自己建立的影像到底对不对，老师在讲解的时候，也只是简单的语言描述，我们根本没法建立清晰可见的三维立体地形对照影像，不利于我们后续课程的学习。2.1二维三维图形的绘制2.1二维三维图形的绘制

为了获取这一地区的点的高度，在原图上，以左下角为坐标原点，图的下边界为X轴，图的左边界为Y轴，分别在两轴上每隔1cm取点，根据图的比例尺1:20000，图上的1cm在实际中表示200米，根据取的点作方格，如图2，再根据地质图上原有的等高线的信息，就可以得到一系列的（X,Y,H）。其中（X,Y）表示在这张地形地质图上某点的位置坐标，H表示该点的高程（通过图上的等高线可以得出来）。下表就是得到的数据：Matlab在构造地质学中的应用2.1二维三维图形的绘制在编程中为了导入的方便，将相同的数据改用txt格式进行导入(h.txt)。2.1二维三维图形的绘制

用原始数据作该地区的地貌图。x=0:200:3000;y=0:200:2400;loadh.txt;meshz(x,y,h);xlabel('X'),ylabel('Y'),zlabel('H')运行结果如图3。明显山区表面很粗糙，效果不好。2.1二维三维图形的绘制

每间隔50加密网格，用最近邻点插值方法作该地区的地貌图，可以使地图的显示效果更好。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）创建一个包含下面代码的M文件：fori=360:-5:0view(i,0);drawnowend（2）保存名称为spin.m。（3）用格线创建球体显示效果：axesm('globe','Grid','on','Gcolor',[.7,.8,.9],'GlineStyle','-')（4）隐藏图形包围盒的边线，进行透视投影：set(gca,'Box','off','Projection','perspective')（5）创建海平面数据网格，使球体不透明：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字符数组与符号运算在数值计算中，最主要的部分就是解线性方程组。由矩阵的定义可知，线性方程组的运算可以转化为矩阵的运算。另外，在地质学中应用比较广泛的还有插值、曲线的拟合等等。趋势面分析：为研究某一地质特征的区域性变化规律和局部范围的异常变化,常用趋势面分析方法来处理地质问题。用matlab处理的过程就是进行多个变量的函数的曲线拟合过程。2.2字符数组与符号运算具体过程为：输入原始数据y建立趋势方程y计算趋势值和偏差值y作分析图和剩余图。在matlab中：一次趋势方程为：W=