实验二 MATLAB的矩阵运算在线视频

这一讲我们来学习实验部分的第2个内容，基于Matlab的平差计算

我们从三个方面来学习这个内容啊，第一 也就是这个矩阵的表达形式

第二 矩阵的运算，第三 我们结合平差计算实例来看一下具体的计算过程

那么首先呢，我们看一下矩阵表达

其实说到矩阵,也就是我们平差计算里面涉及到最多的一种数据表达形式

比如说协方差传播律、协因数传播律，那么这里面就是它的矩阵表达的公式，K、Dxx方差

这都是系数阵和协方差阵，好，那同样还有条件平差模型以及间接平差模型

还有我们的另外两种平差模型，在解算过程当中都涉及到，比如说系数阵A、常数向量A0

还有我们的这个常数项，也就是我们的闭合差项W，等等

Δ,这些都是矩阵,那矩阵的运算其实是很复杂的

我们在你学习的时候呢，应该是我们发现如果是二阶的矩阵，我们手算还是比较容易

但是三阶就稍稍要麻烦一点，那么更高阶的就更不好算了

那么在早期没有这样的计算工具的时候，关于线性方程组的解

也就是怎么样解函数模型，这些都是要花大量的时间，在我们平差课里面去练习的

但是由于现在出现了，像上节课我们学的Excel或者是现在我们这个matlab等等

这样一些计算软件和计算工具的话，这部分的内容呢

在我们上课的过程当中就简化了好多，也容易了好多

我们就看一下，怎么样基于MATLAB软件来进行矩阵的运算

首先，我们看一下matlab它的一个界面，我们知道MATLAB其实就是matrix和laboratory的组合

那其实我们的中文意思就是矩阵实验室，也就是说matlab它的基本的数据对象就是矩阵

那我们看这是它的一个进入系统的界面，面对这个窗口，还有当前目录窗口

还有我们下面这个命令历史记录窗口，但是这些窗口都是可以随机组合的

那上面是它的菜单，对于一些相关的工具，在这呢，我们不去详细介绍它各个部件的使用

我们会结合涉及到的矩阵的相关运算来说它的具体的应用

那么至于它的一些其他的详细的功能，我们可以结合MATLAB工具书去做详细的以及具体的学习

好,那我们首先看一下，在平差里面涉及到的矩阵表达，我们在这给大家总结了几种质证表达方式

包括直接输入、通过数据文件来创建矩阵，通过m文件创建矩阵、以及通过函数创建矩阵

那不管是哪一种方法呢，都要满足我们左这样5条原则：

第一，矩阵元素必须是在中括号里面输入的

这是它的输入的一个基本的要求，也就是形式，你不这样输就是错的

第二，矩阵的同行元素之间用空格或逗号分隔，那就是说，可能矩阵是3行4列

那么这4列4个元素之间用空格或逗号来分开

第三，矩阵的行与行之间要用分号或回车来分隔

第四，矩阵的尺寸不需要预定义，比如说3×4的矩阵，你没有必要先给它预留这样3×4的空间

这跟我们一些C语言等等这样的是不一样的，那些要事先定义矩阵的大小

MATLAB不需要,你有多大的矩阵，输多少个元素就可以了

第五，矩阵元素可以是数值变量，也可以用表达式或者是函数来输入

如果矩阵元素是表达式的话，那么系统就会自动计算出它的结果

然后以矩阵的形式、数据的形式来直接表达出来

那我们先用第1种方式：直接输入法，那么直接输入法又对应两种不同的方法

第一 方法一我们看在命令窗口，按规则输入方式创建矩阵

比如说我们例1,在命令窗口中创建简单的这样一个数值矩阵，那我们看这个矩阵是3×3的

那怎么输入呢？看图1,我直接在它的命令窗口输入A，把矩阵的名命为A就可以写成A等于

刚才按照规则必须是中括号里面，然后是三行用分号隔开，每行三个元素用空格分开

那么最后回车得到的结果就是这样一个3×3的矩阵

第2个例子我们看一下，在命令窗口创建带运算表达式的矩阵

那也就是我们刚才说它原则的时候,这个矩阵呢，除了你逐个向上面输入1 2 3这样的数值之外

还可以把每一个像1 2 3对应位置元素呢，直接写成函数

比如说我们这里面的sin(pi/3),pi是我们MATLAB里面的一个固定的一个常数

就是代表了3.1415926的圆周率，那我们看一下我们对应说出来的结果在右侧图2里面

就显示出来这样一个2×2的矩阵，输入的每一个元素都是函数

直接这样的带有数据的矩阵了，很简单

那我们要注意的第1点，就是在矩阵较大的时候我们可以分行输入

就是比如说刚开始123 456 789，在3和6后面就可以直接按回车直接分行显示

这样也比较直观一些，第二任何矩阵元素内部啊，都不能有空格

比如说我这有一个矩阵元素是23，你在2和3之间出了空格，那它就代表了2和3，不是23了

这是要注意的细节问题，让我们看一下第2种方法，就是在它的工作空间窗口输入，workspace里面

我们看这样一个图，也就是说在work space窗口中，我们点击新建变量按钮

它会出现一个新变量，那我们在这个下面这个长方格新变量这输入它的变量名

然后双击这个变量就会出现右侧这样一个输入窗口

那在这里面你输入你对应的行列对应元素就可以得到矩阵了

这是另外一种方法，直接输入法的里面的第2种方式，我们再看

通过数据文件创建矩阵，比如说我在MATLAB中要创建如下的一个权阵

我们在平台里面是不是有权阵？按经验定权，比如说水准路线的根据路线长度定权是吧

有这样的一个权阵，一个对角阵，那这个怎么来输入呢？

首先,我们第1步说通过数据文件嘛，那我可以在记事本里面输入权阵的矩阵元素，像图5

我把它保存为一个txt文件，定名为quan.txt，这种拼写，那么这样的话我接下来就可以调用这个文件

然后直接输出这个权阵了，为什么用文件方式呢？

因为我在外业用全站仪观测数据的时候，这些数据是记录成一个电子文件的对不对？

那么我们希望把这些数据直接导入到这个软件里面进行计算，那我们平差计算也一样

我把我需要的一些数据也可以直接写成这样的数据文件，以免再平差过程计算过程中手动输入出错

所以有了这样一个记事本之后呢，把这个数据文件第1种方式就可以在命令窗口直接输入quan

直接我们用用load这样的一个命令读入这个权这个数据文件就可以了

那么它会给你显示权是什么，你看直接就出来这样一个权阵

那我们要注意的是呢，你这个记事本数据文件必须是在你当前的工作目录之下

如果不是当前工作目录，是找不到它的啊，这是第1个方法，

第2个方法呢，也就是在这个工作空间work space窗口中导入数据

首先我们点击这个导入数据图标，那么会出现导入数据的这个向导

然后呢，我们按提示操作，当然我们这里面提示你这个记事本里面这个数据

是按空格还是按什么进行分隔的，你选好对应就可以导入这个权这样一个矩阵了，好，这是第2种方法

第3种表达形式就是通过m文件来创建矩阵，那比如说我要创建这样两个系数正阵A和B

那怎么样通过m文件创建呢？那首先呢，我们要将矩阵，按创建原则就是前面我们提到了5条原则

写入到一个M文件里面，那m文件呢，其实相当于是Matlab的一个编程文件

这里面可以写一些脚本，写一些代码

然后把它保存再运行，就可以对应执行它里面的命令了，那在这里面我们输入A B这样两个矩阵

输入方法、原则和我们前面提到的5条是一样的，然后把它保存命名

接下来就可以在命令窗口或程序中直接运行这个m文件

当然找到它的名字直接运行就可以了(文件名)，那这个矩阵就可以在工作空间中产生了A和B

这个很简单，应用起来也很方便

第4种方式就是通过函数创建矩阵，比如说我要创建等精度观测值的权阵，也就是一个单位阵

那我们怎么创建呢？我们的方法就是用这样一个函数我不直接输入数据

我直接用一个eye就是eye(4)，这样一个函数就可以产生P权阵了

那这就是我们MATLAB里面它固有的函数，不需要我们自己去编辑了

它里面已经内嵌了这样的一个函数了，我们直接用就好了，

像这样的函数MATLAB里面有很多，我们在这呢给大家列举了几个常用的

比如说我们这里面产生一个对角阵，对吧？

其实我们这里面也可以用对角阵这个函数来生成P也是可以的，对角阵元素就都是1

以上是我们关于Matlab中矩阵的表达，那现在我们看一下MATLAB中的矩阵运算

在平差里面涉及到大量的矩阵运算，比如说条件平差里面

在这些公式里面涉及到矩阵的加减、矩阵相乘，比如说矩阵的球逆、以及矩阵的转置

那现在我们就依次来看一下，在MATLAB里面如何实现这些项运算

首先，我们看矩阵的加减，比如说有系数阵a和b，我想要让它们俩求和也就是相加，那怎么做呢？

第1步我们可以在命令行，按前述矩阵表达方式输入A和B

然后直接在命令行输入C=A+B，就可以得到C矩阵了，那我们看它就是A、B对应元素相加

很简单，也很容易实现，但我们要注意一点，就是相加减的这个两个矩阵必须是同阶的

就是A是两行三列，B也必须是两行三列，那如果是不同阶它会出现错误

这是加减，那减的话和加一样，我在这就不再说了，我们可以自己去练习一下

第二 矩阵相乘运算，A、B的乘积，你看A和B两个矩阵，怎么算它们的乘积呢？

那如果我们手算的话，可能要写长长的一串计算表达式

那在这里面计算的话，我们同样首先在命令行呢，按前述方式输入矩阵A和B

我们看图12,在这里面呢，这个AB的输入方式我们已经很熟悉了，那有了A和B之后，接下来呢？

我们就可以算它们的乘积，接下来直接输入C=A*B，注意这个乘号运算符是什么呀？

是星号，那这样我们就出现了C这样一个结果，那这个结果对不对呢？

我们看A是2×3，B是3×2，那么乘下来应该是2×2的一个矩阵，也就是C，我们对照一下就是它

那在这我们同样要注意矩阵相乘的要求，在线性代数里面我们提到了相乘的两个矩阵

它必须内接是一致的，内阶数，它的内阶数，如果是A乘B，那就应该是A的列应该等于B的行

那我们在这里面是不是A是三列,B是三行，这样才能相乘，否则是不能相乘的，这个要特别要注意

第三个 矩阵的转置运算，如果有系数阵A,我们在条件平台里面经常要算A的转置，对吧？

A的转置其实就是行列怎么样呀？行列数互换，然后元素位置调换位置

那在这我们在MATLAB里面怎么来实现呢？首先我输入这样一个矩阵A

在工作空间里有这样一个A矩阵变量，然后接下来第2步找它的转置

那这个转置符号也很简单，就是一个单引号，B=A'，就得到B的转置了，很容易实现啊

第四个 矩阵的求逆，这个运算符是inv，如果有权阵，那我们计算观测向量的协因数阵Q

那我们在学权的时候大家应该知道，权阵和协因数阵是互为逆阵的，那我求P的逆阵就可以了

所以首先我们还是输入权阵，在这输入权阵，我们可以按照这样一种方法来输入

diag[], 括号里面按矩阵方式写谁呢？中括号里边是 2 1 3 4，当然要用空格来间隔

这样得到的是权阵的矩阵形式，那这里面和前面的差异就是我们用了一个对角阵函数，对吧

因为权阵是一个对角阵，这也很简单，就不用输入那么多零了，很容易实现

有了权阵之后，接下来我们求逆呀，那求逆我们说是inv这样一个函数

那直接Q等于inv(P)得到就是一个它的逆阵，这样很容易吧

要不然的话我们的4×4的一个矩阵求逆起来是很麻烦的

当然你可能会说这个对角阵求逆就是直接求对角线元素的倒数

我手算也很简单，没错很简单，但是如果它不是这样的一个对角阵呢，对吧？

那那样就复杂了是不是？但在MATLAB里边如果是那样的矩阵，也同样的方法很容易也能出来

这就是它的优势，计算的方便所在。但我们注意它必须是什么样的方阵才有逆

满秩方阵，如果是非满秩的方阵，这个矩阵没有逆，虽然MATLAB算逆很容易

但它也求不出来，因为你的矩阵不满足要求，年原理上不行是求不出来的啊

这是关于我们矩阵的4个运算以及对应的运算符，我们记住以及它的运算的一个规则

记住之后就很容易实现。下面我们再看一下两个计算实例

首先，第1个就是有这样一个观测向量L，它的协方差是DLL, 有观测向量函数Y1和Y2

让我们求Y1、Y2的方差，那你可能看这个题很熟悉，没错

我们在Excel里面对它用Excel软件进行了计算,那在这我们看一下这样一个计算过程

怎么在MATLAB里面实现呢？你可以比一下Excel软件和MATLAB软件之间哪个更方便

看你更喜欢哪一个。那首先呢，我们这里面有两个方法

第1个就是在命令窗口完成，其实关键就是找到系数阵，还有方差阵，还有系数阵的转置

那首先呢，如果在命令窗口实现的话，我在命令行按照矩阵输入原则，输入系数阵K

这就是矩阵输入的方法是吧？再同样输入方差阵DLL，那这个方差阵它是一个对角阵

你可以用对角阵函数，当然你也可以用其他的方法都行，然后就是矩阵的相乘

我们要把这三个矩阵相乘，那么也就是求出来结果我们记为Dyy

它应该是系数阵K乘以方差证DLL，再乘以系数阵的转置，那这个转置你在前面可以不求

在这直接用一个表达式包含了相乘运算和转置运算同时完成，也就是这样一个表达式

Dyy=K*DLL*K',转置是一个单引号，直接得到Dyy

那这里面61和96就是Y1和Y2的方差，那么68呢，就是它们的协方差，是不是很简单？

那方法二 可以利用m文件实现，那也就是说，你可以在m文件里面

建立这样一个函数function，开始协方差计算，然后这里面有这样 总体4行

那么这4行就能实现协方差计算了，第1行输入系数阵，第2行输入方差阵

第3行就是我们刚才的那个计算函数，对吧，K乘以DLL乘以K的转置

那我把它保存，保存之后呢，我们看一下在命令窗口里面，调入XFC，也就刚才我们存的一个文件名

那它会提示你什么呢？输入什么呀？系数阵K，[1 2 ,000 --> 3；2 0 ,000 --> 4]对吧，按照这个输入方式

然后,再输入方差阵D，然后直接它就计算出了结果，那这个方法和前面方法相较而言有什么差异呢？

其实,计算原理是相同的，但是这个呀，如果你发生的话是不是可以函数仍然可以直接用，对不对？

就不需要再怎么样？就是说它可以对应计算不同的同样(类似)的问题，是不是？

数据不同同样的问题都可以计算，就更具有通用性了

接下来我们看一下第2个问题，这样一个水准网，A、B是已知点

然后观测了这样7段高差，让我们算P1、P2、P3点的高程，这就是一个平差解算的问题了

那如果用条件平差法计算，那怎么来计算呢？然后我们第2个问题要评定它的精度是吧

各点高程的中误差，以及平差后高差的中误差，那如果说真正让我们解算的话

其实，我们要通过列条件方程，然后定权，然后依据最小二乘原理解算，对不对？

涉及到了这么多公式，那其实如果你把方程列出来，已知数据有，这个观测数据也有

我们不考虑它的循环，以及这个几次计算，我们只做一次计算的话

如果你有方程、有系数阵、有常数项，那么我们就可以把这些怎么样放到Matlab里面去计算

那就不用我们去手算，首算我们求这个K的时候、Naa的逆、A乘以P的逆再乘以A的转置，是不是很麻烦？

这么多阶，所以呢现在我们怎么办？我们可以首先把它们输入到各自的数据文件里面去

A、W、P、L、f输入到各自的数据文件里面去，然后第2步呢，我们建立如下m函数文件

sz取名，输入后保存，那在这里面我们看首先我们清空工作空间，就把原来里面调用这个函数之后

把工作空间里面其它变量能清除掉，先不要，然后导入谁呢？A、W、P、L、f

就把我们需要的刚才这些矩阵导进来，接着我们看，我们导进来的权

刚刚看了它是一个行向量的是不是？那我们把这个行向量按照生成对角阵的形式构成权阵

那接下来再执行我们的法方程系数阵，A乘P的逆再乘A的转置计算

这里还有求逆运算、转置运算，计算K的时候还是求逆运算、转置运算

以及求V求L^，最后求它的单位权中误差，以及函数的方差、中误差是吧

这些函数都集成到这个m文件里面了，把它保存之后接下来我们在工作空间调入sz

当然，必须把这个sz.m这个m文件保存到当前工作目录下，否则会找不到的

那么我们直接调用，就出现V、LL，也就是观测值的改正数以及观测值的平差值

还有就是单位权证误差以及这个观测值函数的中差，就出来了，这就很容易了

这是条件平差，刚才我们所解决的其实是在方程已经建立起来了

那么我们知道了各个系数阵、常数向量等等，然后把矩阵导进来了

但在很多计算的时候呢，我们希望我们导进来的就是观测数据，就是已知数据

然后我就能直接输出结果，那其实这样的话里面就得包含什么呢

包含你列方程的过程，是不是？包含你定权的过程等等

那这样的话我们就要有一定的编程基础，比如说我用间接平差来解决刚才这个问题

这是间接平差的流程，那在这儿呢，我们就要考虑怎么样从头到尾只给观测数据和这个已知数据

来进行最终结果的求解呢？我们给这样一个计算流程

当然这是一次计算的流程，没有涉及到循环

那你课下呢，你可以试着按照这个流程，如果你有一定的编程基础

那么你可以按照分步把这些实现，其实也是比较容易的

那这，我们就作为一个大作业课下来实现

以上呢，就是我们基于Matlab的平台解算的一个简要介绍

那至于一些涉及到的其它的细节问题呢

MATLAB的专业工具书和你一些编程书来进行详细的应用

我们这呢只起到一个抛砖引玉的作用

谢谢！