东华理工大学毕业设计摘要基于MATLAB 的水准网和测边网平差程序设计摘要MATLAB 是目前在研究机构广泛应用的一种数值计算及图形工具软件，它的特点是语法结构简明、数值计算高效、图形功能完备，特别适合非专业编程员完成数值计算、科学试验处理等任务。以往的测量数据处理方法需要编制特定的处理矩阵运算程序，而且程度复杂，难度大。本文介绍一种基于MATLAB 的水准网和测边网的程序设计方法，与其它算法语言相比，具有编程简单，运算速度快的特点。文中分别阐述了水准网和测边网程序的理论基础、实现步骤和运行结果。通过实例的分析，总结出利用MATLAB 对测量数据处理有很大的应用价值，它缩短了编程的时间，提高工作效率。关键词：MATLAB；水准网；测边网；程序设计： 1 东华理工大学毕业设计摘要ABSTRACT MATLAB is one species of numerical-values calculation and graphic tools software which is widely used to apply at research institutions at present. The particularities are: concise grammar-structure、highly efficient in numerical values calculating、complete function of graphs、especially it is adapted to evildoing professional programmer to accomplish the tasks that are numerical-values calculating and scientific experiments treating. The ancient methods of measured data-processing need establishing special proceedings of treating matrices operation, moreover, it is complex and greatly difficult. This article introduces one programming method dealing with leveling and measuring edge network based on MATLAB. Compared with other algorithm language, it has particularities which are simply programming and quickly operating. The article separately expatiate the theories basics、realizing steps and running results at leveling and measuring edge network. With the analysis of examples, it has prodigious application value in measured data-processing by use of MATLAB. Moreover, it shortens programming time and improves working effectiveness. Key words：MATLAB；leveling network；measuring edge network；programming ： 2 东华理工大学毕业设计目录目录................................................................ .................................. 绪论.................................................................. 4 ................................................... 1. MATLAB 软件简介................................................... 5 2．MATLAB 在测量平差中的应用.......................................... 6 .......................................... 2.1 测量平差原理的概述................................................ 6 2.2 平差程序总体方案.. (7)
3.1 程序的功能........................................................ 8 3.2 水准模型网的间接平差.............................................. 8 3.2.1 “权”值的确定................................................ 8 3.2.2 水准路线的平差计算............................................ 9 3.2.3 精度评定..................................................... 11 3.3 水准网间接平差程序信息设计....................................... 11 3.4 水准网程序与使用说明........................................... 12 3.
4.1 水准网程序流程图............................................ 12 3.4.2 水准网程序的使用............................................ 12 3.5 案例............................................................
13 ................................................ 4. 测边网平差程序设计................................................ 15 4.1 数学模型........................................................ 15 4.1.1 误差方程和法方程的组成....................................... 15 4.1.2 边长观测的权................................................ 15 4.1.3 解算法方程.................................................. 16 4.1.4 精度评定....................................................
19 4.2 测边网平差信息设计............................................ 20 4.2.1 主要的技术要求.............................................. 21 4．3 利用MATLAB 的绘图语句绘制网图................................... 21 4.4 测边网程序和使用说明............................................ 22 4.5 程序代码说明：................................................. 23 4.6 程序的使用算例.................................................. 25 结论................................................................
29 ............................................................... 致谢...............................................................
30 ............................................................. 参考文献.............................................................
31 ............................................................... 附录一...............................................................
32 ............................................................... 附录二 (36)
附录三............................................................... 46 ............................................................... 3 东华理工大学毕业设计绪论绪论作为一名测量技术人员，如果不掌握一门PC 机编程语言与便携计算工具，要想提高测量工作的效率几乎寸步难行。测量需求的多样性与复杂性，造就了测量计算鲜明的个性化特点，这就是在商业测量计算软件高度发达的今天，掌握一种实用的程序语言进行编程计算仍有广泛的市场需求的重要原因。当今较流行的计算机程序语言基本上都是基于Windows 的，例如Turbo 这些程序语言的优势是基于对象及可利用Pascal,Visual Basic,Visual C,Borland C++等，Windows 丰富的系统资源，应用它们可以开发出界面非常丰富和友好的应用程序，其劣势主要有以下几点：1.Windows 程序都非常庞大，学习并熟练掌握它们并非易事。2.虽然市场上已有的多种专用的测量平差软件都是采用 C 语言开发的，但这些软件价格都比较贵，而且都带有加密狗，一次只能供一个用户使用。出于商业目的，开发商不会公开程序源代码，这为修改程序功能以适应用户的特殊需求带来了不便。 3.在测量生产中，经常需要根据工程的实际情况进行一些个性化的数值计算工作，这些数值计算工作无固定模式，这就需要求测量技术人员最好能熟练掌握一种适用于数值计算的程序语言，以便提高测量计算的效率。4.C 语言的数值计算语句不够丰富，例如，在测量平差计算中，经常需要进行的矩阵运算，尤其是解法方程的矩阵求逆不能直接使用语句实现，而必须应用计算机算法编程实现。如果不是基于商业软件开发，只为满足实际测量工作计算需要，则C 语言的劣势就变成了MATLAB 语言的优势。4 东华理大学毕业设计MATLAB 软件简介 1. MATLAB软件简介MATLAB是从Matrix(矩阵)和Laboratory(实验室)各取前3个字母组成的,意思是矩阵实验室,是美国MathWorks公司于20世纪80年代中期推出的一种交互式、面向对象的科技应用软件,是一个为科学和工程计算而专门设计的高级交互式软件包。MATLAB 集成了图示与精确的数值计算,是一个可以完成各种计算和数据可视化的强有力工具, 其优秀的数值计算能力和卓越的数据可视化能力使其很快在数学软件中脱颖而出,成为以矩阵运算为主要工作方式的线性代数、概率论和数理统计、自动控制、数字信号处理、动态系统仿真等领域教学和科研工作者的有力武器。随着该软件自身的发展及市场的需求,其功能日趋完善,前其最高版本7.0版已经推出,随着版本的不断升级,它的数值计算及符号计算功能得到了进一步完善。MATLAB是以矩阵作为数据操作的基本单位,矩阵的生成、运算、转置、求逆等非常简单。在MATLAB环境中,不需要对创建的变量对象给出类型说明和维数,所有的变量都作为双精度数来分配内存空间,MATLAB将自动地为每一个变量分配内存。MATLAB语言起源于矩阵运算,并已经发展成为一种高度集成的计算机语言,它提供了强大的科学运算、灵活的程序设计流程、高质量的图形可视化与界面设计、便捷的与其他程序和语言接口的功能。MATLAB系统主要包含5 部分的内容:MATLAB 工作环境、Mablab 数学函数库、MATLAB语言体系、句柄图形、MATLAB应用程序接口(API)。MATLAB系统主要功能包括:数值计算功能、符号计算功能、数据分析和可视化、文字处理功能、SIMULINK动态仿真功能。同时,MATLAB又是开放的,除了内部函数之外,所有的MATLAB 主包文件和各工具包文件都是可读可改的源文件,用户可以作为参考掌握其用法,并可对其修改以适应自己的需要,也可加入自己编写的文件构成新的工具包。例如,随着GPS 的广泛应用,Orion Dynamics and Con2t rol Corporation 、Constell Inc. GPSSoft LLC、NavsysCorporation等多家公司都相应开发出了适于GPS数据处理的MATLAB 工具箱。MATLAB 是一个集数值计算、图形管理、程序开发于一体的功能十分强大的系统。将MATLAB 应用于测量数据的处理是一件非常有意义的工作。
Mo2hamed 等曾成功地在MATLAB 系统中利用白滤波技术研究动态解算GPS 载波相位信号的模糊度问题。因为测量数据的处理特别是测量平差主要应用矩阵运算,而MATLAB 又特别易于做矩阵运算,因此,研究开发基于MATLAB 的测量平差方法具有极好的应用价值。
5 东华理大学毕业设计MATLAB 在测量平差中的应用2．MATLAB 在测量平差中的应用测量平差数据处理主要是基于矩阵的运算,常用的矩阵运算主要是矩阵的生成、转置、求逆和矩阵求广义逆等。在MATLAB环境中,不需要对创建的变量对象给出类型说明和维数,所有的变量都作为MATLAB中的M文件的语法与其他的高级语言类似,是一种程序化的编程语言,同时也是一种解释性的编程语言,即逐行解释运行程序,使程序容易调试,计算更为简捷,而且对于平差原理理解和掌握变得更容易。另外,MATLAB语言与数学语言比较接近,更容易掌握和理解。2.1测量平差原理的概述测量平差的函数模型有条件方程和观测方程。以条件方程为函数的模型的最小二乘平差称为条件平差；在条件方程中，根据需要如果还设有一定数量的未知数，则称为附有参数的条件平差；以观测方程为函数模型的最小二乘平差称为间接平差；如果观测方程中的某些参数不独立，则这些不独立参数必然存在一些条件，称这种平差模型为附有条件的间接平差。本文的两个程序都采用间接平差模型。对于一个实际平差问题,根据所选参数的个数、选什么量为参数以及参数之间是否函数独立,经过仔细推敲可以发现附有条件的间接平差模型本身就是各种经典平差模型的概括模型,其余的经典平差模型,如条件平差模型、间接平差模型、附有未知数的条件平差模型和附有限制条件的条件平差模型都是它的特例。间接平差的公式汇集：间接平差模型为? ? ? T V PV → min ? ? V = BX ?l 系数矩阵B 满秩，即rank(B)=t 法方程及解为：（2－1）$ N bb x ? f e = 0, ( N ?1bb = B T PB , f e = B T Pl ) （2-2）（2-3）（2-4）（2-5）（2-6）（2-7）$ x = N ?1bb f e $ 参数的平差值：X = X 0 + x 观测量的平差值：$ = L + V L 单位权中误差：V T PV σ0 = n?t D 平差参数的协方差阵：X X = σ 2 0 N ?1bb
6 东华理大学毕业设计MATLAB 在测量平差中的应用平差函数的协方差阵：φφ = σ 2 0 F T N ?1bb F Q$ $ （2-8）2.2 平差程序总体方案MATLAB 号称为全球工程师的共同语言，其语法和C 语言相似，但它有强大的数值计算和绘图功能，这使之在工程应用方面的计算更出色，本文就基于这种程序设计语言环境设计一个控制网平差程序。该程序包含了一个高程控制网平差程序和测边网平差程序。本程序适用于各种等级的高程网和测边网，程序在设计过程中，始终考虑数据的储存量。因而本程序不储存误差方程的系数和常数项，对待定点数较多的平差网，组成法方程的系数矩阵是个稀疏矩阵，如待定点的编号恰当，法方程的系数会集中在主元系数的两侧形成带状。为减少法方程系数的储存量，只要按行储存下三角阵或按列储存上三角阵中第一个非零系数起的系数，就是通常叫做维变带宽储存方法。
7 东华理大学毕业设计水准网平差程序3. 水准网平差程序3.1 程序的功能本程序适用于二、三、四等水准网平差计算，平差的水准网可以是独立的、也可以是附合网，其主要功能是完成水准网的平差计算和精度评定计算。平差计算采用间接平差法，以归心的观测值为高差，以未知点高程为未知参数。精度评定计算包括计算单位权中误差和每个待定点的高程中误差。3.2 水准模型网的间接平差 3.2.1 “权”值的确定当在相同的条件下进行水准测量时,其精度是相同的,因而观测结果的可靠性也是同样的。但如果在不同的条件下进行水准测量时,高程的精度就有所不同,此时称为不等精度观测,所求出的未知量的值、高程的最或是值并对其精度进行评定时,就需要“权”了。由于观测的不等精度,因而观测值的可靠程度不同,求未知量的最或是值时,这样的一个因素就必须考虑了,这个因素是:可靠性大的某观测值,其精度高,对测量的最后结果的影响也就越大。此时用“权”值来表示观测值的可靠程度,那么,“权”值愈大,观测值的可靠程度就愈高。另外,在观测过程中,观测值的中误差愈小,观测结果愈可靠,它的“权”值就愈大。因而,根据中误差来确定“权”值是非常适当的。设以Pi 表示观测值Li 的“权”,m 为中误差,则“权”值的定义为: (3-1) 式中:A 为任意的正常数,在一组观测值中为一个定数。在实际测量中,
通常是观测值的中误差事先并不知道,因而必须先确定观测值的“权” , 然后才能求出未知量的最或是值。此时可以利用距离(S)或测站数(N)来确定观测值———高程的“权”。根据偶然误差传播定律,各观测点高程Hi的中误差mi由测站数Ni确定时,则有: mi = m N i 式中:m 为一组观测值的中误差,为一个定数. 由(3-2-1)、(3-2-2)两式可得: (3-3) (3-2) 8 东华理大学毕业设计水准网平差程序同样可得出: 式中:C为定数, si 为测距. （3-4）由(3-3)、（3-4）两式可以得出这样一个结论:当测站观测高差等精度时,观测总高差的“权”与测站数或距离成反比。3.2.2 水准路线的平差计算1．附合路线的平差计算假定在图1 示的A 、B 两水准点之间布设一条水准路线,A 、B 两水准点的高程为已知,分别设为H A 、H B 、n1 、n2 ?、C为中间水准点。假定观测了所有的点的高程,现拟求C点的高程H C 的最或是值。
H C 可由水准路线A →C、B →C分别观测的高差Δ hAC 、Δ hBC 计算得出,由此而得到的观测高程分别设为Hc1、Hc2,其值为: Hc1= H A +Δ hAC ;Hc2= H B +Δ hBC 当Hc1、Hc2在不等精度条件下观测得出时,它们的“权”也不同,分别设为Pc1、Pc2,这样C点的高程H C 的最或是值为: Hc = Pc1 H c1 ? Pc 2 H c 2 Pc1 + Pc 2 (3-5) 根据A 点的高程H A ,A →C水准路线观测的高差Δ hAC 以及B →C水准路线观测的高差Δ hBC ,可推算出B点的观测高程H B 为: H ′B = H A +Δ hAC -Δ hBC 水准路线A →B 的高程闭合差为: f h = H ′B - H B = H c1 ? H c 2 (3-6) 由（3-6)式得到: H c 2 = H c1 - f h 由（3-3)式得到: Pc1 = C C 、Pc 2 = N （N AC 、N BC 分别表示水准路线A→C、B N AC BC →C 的测站数,水准路线A→B的测站数N AB = N AC + N BC ) 9 东华理大学毕业设计水准网平差程序将上述表达式代入(3-2-5)式中,得到: (3-7) 如果以水准路线A→C的距离S BC 、→C的距离S BC 、→B的距离S AB ( S AB = S AC + S BC ) B A 来确定高程观测值的“权”值时,同样可以得到: (3-8) 图3-1 水准路线图2．闭合路线的平差计算闭合路线的平差计算原理与附合路线相同,因而(3-7)、(3-8)两式的结论适用于闭合路线的平差计算。(3) 具有一个结点的水准网的平差计算如图2所示为具有一个结点的水准网,B,C,D,?为已知高程水准点,B→A,C→A,D →A,?为水准路线,则接点A的高程最或是值为: P H + P H + PA3 H A3 + L H A = A1 A1 A 2 A2 = PA1 + PA2 + PA3 + L ∑P i =1 n i =1 n Ai H Ai (3-9) Ai ∑P 式中H A1 , H A2 , H A3 L 分别为水准路线B→A,C→A,D→A,?计算A的观测高程,各高程相应的“权”值为PA1 , PA 2 , PA3 L 设H A1 , H A2 , H A3 L 的算术平均值为H 0 A ,各高程观测值与H 0 B 的差值分别为δA1, δA2,δA3,?,则有: 10 东华理大学毕业设计水准网平差程序(3-10) 将(3-10)式代入(3-9)式得到: (3-11) 当以测站数和距离来确定“权”值时,(3-11)分别可以转化为: (3-12) (3-13) 上述结论也可应用于小三角水准网平差计算。3.2.3 精度评定单位权中误差：σ 0 = V T PV n?t （3-12）（3-13）（3-14）平差参数的协方差阵：DX X = σ 2 0 N ?1bb Q$ $ 平差函数的协方差阵：φφ = σ 2 0 F T N ?1bb F 3.3 水准网间接平差程序信息设计1．数据文件的组织下面给出一个水准网输入数据文件的例子：3 3 6（已知点个数、未知点个数、观测值个数）101 102 103 104 105 106 （点号）34.788 35.259 37.825 （已知点高程）104 101 1.625 4.5（起点点号、终点点号、高差观测值、距离观测值）101 102 -0.418 3.1 105 102 0.714 3.4 102 103 1.243 3.8 106 103 -0.577 4.3 11 东华理大学毕业设计水准网平差程序103 101 -0.786 2.5 （其中编号数组未知点在前，已知点在后）2．水准网平差变量约定表3-1 变量约定表变量名ed dd sd gd pn h0 be en h1 s 说明已知点个数未知点个数总点数观测值个数点号已知点高程起点点号终点点号高差观测值距离观测值3.4 3.4.1 水准网程序与使用说明水准网程序流程图图3-2 水准网流程图程序的全部代码见附录一。 3.4.2 水准网程序的使用本程序使用MATLAB 的矩阵功能计算法方程，在运行程序前首先要有其始数据。其始数据是一文件的形式保存在磁盘中，文件的格式在上文已经说明过，编好文件后，12 东华理大学毕业设计水准网平差程序以后缀名为.TXT 的形式保存。执行时在MATLAB 命令窗口直接键入文件名即可。3.5 案例如下图水准网，104、105、106 为已知点，101、102、103 为待定点，已知点的高程分别为34.788，35.259 ，37.825 。观测高
差和观测路线长度分别为：h1=1.652, h2=-0.418 ,h3=0.714 ,h4=1.243, h5=-0.577, h6=-0.786. s1=4.5, s2=3.1, s3=3.4, s4=3.8, s5=4.3, s6=2.5. 图3-3 水准网图首先编数据文件，命名为data1.txt. 数据的格式如下：3 3 6 101 102 103 104 105 106 34.788 35.259 37.825 104 101 1.652 4.5 101 102 -0.418 3.1 105 102 0.714 3.4 102 103 1.243 3.8 106 103 -0.577 4.3 103 101 -0.786 2.5 进入MATLAB界面，在命令窗口直接输入level3运行程序。弹出如下窗口13 东华理大学毕业设计水准网平差程序图3-4 数据读入文件选择data1.txt即可运行出如下结果：图3-5 计算结果在图3-5 中，分别输出了高程的平差值及精度。结果是一文本的形式保存，用户可对它进行编辑。14 东华理工大学毕业设计测边网平差程序设计4. 测边网平差程序设计4.1 数学模型4.1.1 误差方程和法方程的组成控制网中的观测值为边长，误差方程非零项最多为4 个，所以误差方程系数矩阵采用压缩格式进行储存。可采用以下的方法：A(m,n)=>A(m,9) 其中，m 为观测值个数，n 为未知点个数的两倍。改进后的A 阵格式为Ai =（编号1，系数1，编号2，系数2，L 编号4，系数4，常数项）共9 列。即只存储误差方程的4 个非零参数系数。法方程系数阵N A 为对称阵，在存储时，只需要存其上三角部分就可以了。其占用的空间为：sum = n(n + 1) 2 现有 A 阵：A=（编号1，系数1，编号2，系数2，L 编号4，系数4，常数项）其中偶数项为系数，加上最后的A9 为常数项，在组成法方程时，从A2 开始分别与剩下的偶数项以及常数项相乘，然后再用A4 与剩余的项相乘，一直到A8 为止，这样就完成了N A = AT PA 的过程。需要注意的是：若A1，A3，A5，A7 小于零，则表示该点已知点，不参与法方程的组成。4.1.2 边长观测的权边长观测的精度一般与其长度有关，定权公式为psi = σ 02 σs 2 i (i = 1, 2,L , n) 式中σ si 2 为所测边长si 的方差，σ 0 2 为任意选定的单位权方差。为了定权psi 必须已知测边的先验方差σ si 2 ，但精确的已知是十分困难的，一般采用厂方给定的测距仪精度，即σ s = a + bSi i 式中，a 为固定误差（单位mm），b 为比例误差（单位：ppm），Si 为边长（单15 东华理工大学毕业设计测边网平差程序设计位km）。 4.1.3 解算法方程由于法方程是对称正定阵，因此，可采用改进的平方根法进行解算。平方根法是对称正定矩阵非常有效的三角分解方法，设A 为n 阶方阵，如果其所有顺序主子式均不为零，则其存在唯一的分解式：A=LDR ?1 ? ?1 ? d1 K 0 ? ? ? ? ? ，D= ? M O M ? ，R = ? L = ? l2 ? ? ?M ? ? O ?0 L d ? ? ? ? n? ? ?l L l ? ? n , n ?1 1 ? ? ? n1 r12 K r1n ? ? M ? O rn ?1, n ? ? ? 1 ? 其中由于此住A 对称性，得L = RT ，又根据A 阵正定的性质，可证明D 均为正数。现在设? d1 ? ? d1 ? ? ? ? ? O D= ? O ? ? =? ? ? ? ? dn ? ? d d1 ? ? 1 ? d1 ? ? ? O ? ? ? ? d d1 ? ? 1 即则为方便，记为D = D2 D2 A = LDLT = LD 2 D 2 LT = ( LD 2 )( D 2 LT ) = LL A = LLT 1 1 1 1 T 称为Cholesky 分解，即正定对称矩阵的平方根分解法。解AX = b 等阶于求解两个三角方程组：LY = b 和LT X = Y 在用平方根分解法计算时，需要进行n 次开方运算。为了避免开方，可以直接采用对称正定的A = LDLT 分解式对平方根法进行改进。从而解方程组AX = b 可以按如下步骤进行：把A 分解成A = LDLT ，则AX = b 变成( LDLT ) X = b ，即等价于? LY = b ? T ?1 ?L X = D Y 由此可以解出X 和Y。这称为改进的平方根法，在计算中避免了开方运算。平方根法和改进的平方根法的计算量和存储量比消去法节约近一半，而且不需要选主元，能得到比较精确的数值解。16 东华理工大学毕业设计测边网平差程序设计法方程用改进平方根法解算的过程如下：(1)分解：C = S T D ?1 S 其中? S11 K S1n ? ? S11 ? ? ? ? ? S =? O M ?, D = ? O ? ? ? ? ? Smn ? Smn ? ? ? ? ? ? 1 K c1n ? ? s12 ? ? O M ? = ? s11 ? L cnn ? ? M ? ? s1n ?s ?
11 ? ? ? ? ? S11 K S1n ? ?? ? O M ? ?? ? ?? S mn ? ?? ? 1? ? ? c11 ? ? M ?c ? n1 O sn ?1, n sn ?1, n ?1 L s1 j = c1 j c2 j = c3 j = 纯量计算公式为i =1 ? cij ? i ?1 s s s1 j = ? ki kj cij ? ∑ i > 1, j ≥ 2 ? k =1 skk ? s12 s1 j s11 s13 s1 j s11 s2 j , + s2 j = c2 j ? s12 s1 j s11 , j≥2 s13 s1 j s11 ? s23 s2 j s22 , j≥3 s23 s2 j s22 + s3 j , s3 j = c3 j ? （2）求逆R = S ?1 ? r11 K r1n ? ? ? R=? O M ? ? ? rnn ? ? 由RS=I 得? r11 r12 L r1n ? r22 L r2 n ? ? O M ? rnn ? ? ? s11 ? ? ??? ? ? ? ? ? ? s12 L s22 s1n ? ? ? L s2
n ? ? =? ? ? O M ? ? snn ? ? 1 0 L 0 1 L M M O 0 0 L 0 0 M 1 ? ? ? ? ? ? 17 东华理工大学毕业设计测边网平差程序设计纯量计算公式：rii = 1 sii r11 s12 s 22 r12 = ? r13 = r1n = ? ( r11 s13 + r12 s 23 ) s33 ? ( r11 s1 n + L + r1( n ?1) s ( n ?1) n s nn r22 s12 s33 r23 = ? r24 = r2 n = 通式为? ( r22 s 24 + r23 s34 ) s 44 ? ( r22 s 2 n + L + r2 ( n ?1) s( n ?1) n s nn 1 ? ?rii = s ii ? ? j ?1 ? ∑ rik skj ? , ?rij = ? k =1 sij ? ? （3）求积j >i Q = ( S T D ?1S ) ?1 = S ?1 D( S ?1 )T = RDRT ? r11 K r1n ?? s11 ? ? r11 ? ? ?? ?? ? Q = RDR = ? O M ?? O ?? M O ? ? ?? ? ? rnn ?? snn ? ? r1n L rnn ? ? ? T ? s11 r11 s12 r12 L ? s22 r22 L =? ? O ? ? s1n r1n ? ? ?? s2 n r2 n ? ? M ?? ?? snn rnn ? ? r11 rnn M r1n rnn M O rnn rnn L ? ? ? ? ? ? 18 东华理工大学毕业设计测边网平差程序设计4.1.4 精度评定（1）坐标改正数以及单位权中误差的计算m0 使用上三角一维数组形式存储坐标改正数的公式为：δ xi = ?∑ q ji w j ? ∑ qij w j , i = 1,L , n j =1 j =1 i ?1 n $ 其中，n=2 × dd, xi 的单位是cm. 平差值：$ X = X0 +x 写成分量的形式，为X i = X 0i + δ xi 如果近似坐标的误差较大，或网形较大，平差的结果不会精确，这时，就需要进行迭代平差，直到两次平差间互差在允许值内。由测量平差理论：σ = 同样可得到单位权中误差：V T PV n?t m0 = 其中，m-n 观测个数减去未知点个数；[ PVV ] m?n m = m1 + m2 + m3 , n = 2 × dd + ST , ST ? 方向观测的测站数[ PVV ] = [ pll ] + [δ x ω ] （2）点位误差椭圆误差椭圆表示了网中点或点与点之间的误差分布情况如图。在测量工作中，常用的误差椭圆对布网方案作精度分析。绘制误差椭圆只需要三个数据：椭圆长半轴a，短半轴b 和主轴方向φ,其求法为1 ? (σ x 2 + σ y 2 + (σ x 2 ? σ y 2 ) 2 + 4σ xy 2 ) ? ? 2 ?, 1 2 2 2 2 2 2 2 ? b = (σ x + σ y ? (σ x ? σ y ) + 4σ xy ) ? 2 ? a2 = tan 2? = 2σ xy σ x2 ? σ y2 19 东华理工大学毕业设计测边网平差程序设计图4-1 误差椭圆的表达顾及方差与权倒数的关系，得? (Qxx + Qyy + (Qxx ? Qyy ) 2 + 4Qxy 2 ) ? ? 2 ?, 2 σ0 2 2 2 ?
b = (Qxx + Qyy ? (Qxx ? Qyy ) + 4Qxy ) ? 2 ? a = 2 σ 02 tan 2? = 2Qxy Qxx ? Qyy 根据上述的理论，我们实际要求的是mxi 、m yi 、m xiyi 。只要得到了这些元素，就能依照上面的公式来求得椭圆的元素了。4.2 测边网平差信息设计外业测量的数据首先应进行预处理，包括测站平差、归心计算、观测值归化到椭球面的改正、椭球面归化到高斯平面的改正等，然后将预处理后的数据输入到以后缀名为TXT 的文本文件中，该数据文件的组织格式如下所列：表4-1 数据组织格式次序1 2 3 4 5 6 7 8 9 内容已知点个数ed,未知点个数dd,[控制参数] 点号pn。先输入已知点编号，各点输入顺序无要求已知点坐标,x0,y0,x1,y1,x2,y2… 测量边的个数m1 边长的固定误差ms,边长的比例误差pp。单位分别是：和×0.000001 cm 边长的起始点号e,终点点号d,边长sid。每一条边一行，依次列出推算近似坐标的路线经过的边数推算近似坐标的起算已知点坐标（按顺时针）推算近似坐标的路线经过的边的边号20 东华理工大学毕业设计测边网平差程序设计 4.2.1 主要的技术要求表4-2 测边网的技术要求等级二等三等四等一级小三角二级小三角平均边长（km）测距中误差（mm）9 5 2 1 0.5 ±30 ±30 ±16 ±16 ±16 测距相对中误差1∕30 万1∕16 万1∕12 万1∕6 万1∕3 万4．3 利用MATLAB 的绘图语句绘制网图1. 网形的绘制由于网形图与误差椭圆绘制在同一个图形上，因此必须对误差椭圆进行放大，在本文的程序中，使用了inputdlg 对话框输入语句，其中，确省的放大倍数为100。在程序中，使用了ed dd pn m1 x y e d sid ai bi fi 等变量，其意义与前面的变量相同。对绘制的网图有效放大和缩小功能，即点击鼠标左键放大图形，点击右键缩小图形，利用MATLAB 菜单本身的功能，可以将该图形输出为各种图形文件格式。2 误差椭圆的绘制无论多么复杂的图形，其基本单元还是点和线。换句话说，只需要利用基本元素的点或线，通过各种组合，也能画出复杂的图形。MATLAB 中没有提供直接绘制椭圆的命令，因此可以直接利用连线来画椭圆。测量中描述误差椭圆用长半轴A，短半轴 B 和方位角FI 三个量。在如图4-1 的x′Oy′中直角坐标系中，椭圆的标准方程为x′2 y′2 + =1 A2 B 2 如果以角度i 为变量，则椭圆的标准参数方程为x′ = Αcosi ? ? y′ = Βsini ? 设在测量坐标系xOy 中椭圆长半轴的方位角为?0 ，则有（4-2）（4-1）? x ? ?
cos ?0 ? ? = ? sin ? ? y? ? 0 用参数方程代入，得到? sin ?0 ? ? x′ ? ? ? cos ?0 ? ? y′ ? ? （4-3）21 东华理工大学毕业设计测边网平差程序设计x = A cos ?0 cos i ? B sin ?0 sin i ? ? y = A sin ?0 cos i + B cos ?0 sin i ? 图时，需要调换x，y 坐标。（4-4）测量坐标系与MATLAB 或AUTOCAD 绘制的数学坐标系的x,y 坐标轴不同，绘在上式中，i 取不同的值，就有一组（x，y），只需要将这些点连接起来，就可以绘出一个椭圆。图4-2 误差椭圆的参数方程4.4 测边网程序和使用说明使用本程序的全部数据都按规定的格式编辑成数据文件储存在磁盘上。数据文件的编辑取决于平差网型和观测值的编号，为此，先绘制平差网的略图，在图上标明各项数据的信息。以下就测边网的点号，边编号，输入数据，输出成果，运行程序等问题作简明说明。1 点号和观测边的编号已知点和待定点的编号为pn（取三为数），已知点在前，未知点在后，其顺序无要求。但为了减小法方程系数的带宽，应使相临的待定点编号的差数尽可能小。平差网的编号见图4-3。 2.推算近似坐标的路线近似坐标的路线是用户在测边网略图上指定出来的。如图4-3 的箭头，就是表示推算路线。路线的两个起算点必须为已知点，从两个已知推算出的第一个未知点开始，选择观测边，由观测边和已求得的近似坐标或已知坐标推算出观测边所对的未知点。本程序是按推算的三个点A、B、P 顺序为顺时针。22 东华理工大学毕业设计测边网平差程序设计图4-3 测边网略图 2 数据的输入为了在程序运行中数据的传递，定义了一些全局变量,参见4-1 表。（1）简单变量（2）数据文件外业测量的数据首先应进行预处理，包括测站平差、归心计算、观测值归化到椭球面的改正、椭球面归化到高斯平面的改正等，然后将预处理后的数据输入到以后缀名为TXT 的文本文件中，该数据文件的组织格式如表4-1。3 编辑。4.5 程序代码说明：程序的总体框架输出成果本程序的计算成果是以文件的形式输出到文本文件中，用户可以在文本中查看和图4-4 测边网总体流程图23 东华理工大学毕业设计测边网平差程序设计1．数据读入块本模块的功能是打开一个*.txt 的数据文件，同时生成一个*out.txt 的文本文件，记录用户数据和输出成果用。程序调用fscanf 函数把文件中的数据赋值给相应的变量，这些变量是后面计算的数据依据。2．近似坐标的计算测边网的观测数据是边长，所以在近似坐标的计算时只能用测边交会计算。以及A、B 两点坐标及A、B 到P 点的距离b、a，c 为A、B 两点的距离，A、B、P 三点按顺时针排列。则P 点的坐标计算公式如下：b2 + c 2 ? a 2 , f = b2 ? e2 2c X P = X A + e cos α AB ? f sin α AB ? ? YP = X A + e sin α AB + f cos α AB ? e= 计算近似坐标的流程图入图4-5。图4-5 计算近似坐标流程图3．误差方程和法方程的形成程序中，用数组 a 来存储误差方程的编号和系数，a(i,9)存储常数项，w 和 c 分别存储法方程系数和法方程常数。4．解算法方程。函数流程图如图4-6。函数的代码见附录三。5．精度评定本模块包括坐标改正数、单位权中误差和误差椭圆的计算。程序中定义了dxy 坐（标改正数）、pvv（即存储[ PVV ] = [ pll ] + [δ x ω ] ）uw0(单位权中误差)等。同时计算、出误差椭圆的三个参数长半轴ai，短半轴bi 和主轴方向fi。24 东华理工大学毕业设计测边网平差程序设计图4-6 平方根法求逆程序框架图6．控制网的绘制程序中调用inputdlg 函数来打开一个对话框来输入误差椭圆的参数；text 函数用来对点号标注；用ploth 函数及控制参数绘制线和误差椭圆。程序的完全代码见附录二。4.6 程序的使用算例有测边网如图4-3 所示。网中A、B、C、C 点已知，其余为未知点，现用某测距仪观测了13 条边长，测距精度σ s = (3mm + 1× 10?6 S ) 。起算数据和观测数据如下：表4-3 已知点坐标点名 A B X 53743.136 47943.002 Y 61003.826 66225.854 点名C D X 40049.229 36924.728 Y 53782.79 61027.086 表4-4 观测值编号1 2 3 4 5 边观测值5760.706 5187.342 7838.88 5483.158 5731.788 编号6 7 8 9 10 边观测值8720.162 5598.57 7494.881 7493.323 5438.382 编号11 12 13 边观测值5487.073 8884.587 7228.367 25 东华理工大学毕业设计测边网平差程序设计1.编辑数据文本文件如下图图4-7 数据文件2.在MATLAB 命令窗口键入nnbb 执行程序，运行中会弹出一个对话匡提
示用户输入误差椭圆的放大比例，默认为100，本例选择500。如下图图4-8a 对话框图4-8b 对话框3.计算出的结果如下：在图4-9a 中，第1、2 行分别是已知点和未知的个数；第4 行是点的编号；第5 至第8 行是已知点的坐标；第10 行是观测值的个数；第12 行是测距的固定误差和比例误差；至25 行是观测边的起点号、13 终点号和观测边长；至39 行27 （接到图4-9b）是点号转换为计算机顺序后的观测边起点号、终点号和观测边长。接图4-9a，在图4-9b 中，第42 至49 行是推算的近似坐标；第52 至64 行是计算的误差方程系数和常数。在图4-9 中，第67 至70 行是法方程的系数（上三角一维存储）；第73 至76 行是求逆后的法方程系数；第79 至86 是坐标改正数和坐标平差值；第89 至92 行是误差椭圆的参数；第95 行是单位权中误差。26 东华理工大学毕业设计测边网平差程序设计图4-9a 计算结果图4-9b 计算结果27 东华理工大学毕业设计测边网平差程序设计图4-9c 计算结果 4.输出的控制网图和误差椭圆图如下: 图4-10 控制网图28 东华理工大学毕业设计结论结论由本论文的两个程序例子解算过程不难发现，程序的平原理和其它高级语言编写的程序大致相同的。但MATLAB 语言本身就有独特的数值计算功能和图形绘制功能，这使它不用编写专门的函数而直接处理和计算测量的数据，这就大大缩短了编程的时间。MATLAB 软件将人们从乏味的Fortran、C 编程中解放出来，使他们真正的把精力放在科学研究的核心问题上。但用MATLAB 语言编程本身也有它的不足之处：1.由于MATLAB 是一个解释器，会逐行对程序代码进行解释后执行，当要处理的数据量很大时，计算机的运行速度明显变慢了好多。2.用MATLAB 进行界面开发时没有像其他面向对象语言方便。针对上面的问题，提出一些见解：1.为了提高整体程序的执行效率，应尽量多使用向量化的运算，而避免或少用for 循环、while 循环。事实上，如果程序是纯粹的数值运算、而没有使用for 循环或while 循环，那么其执行速度将会接近于纯粹用C 语言来写的程序代码。2.虽然MATLAB 是一个完整的科学计算与可视化环境，但在多数情况下，MATLAB 如编程的代码执行效率不佳、你希望隐藏你的程序代码以保护产权、或您想要进行任何只有在某些高级语言程序才能做到的事，这时就需要用MATLAB 所提供的应用程序接口。29 东华理工大学毕业设计致谢致谢本论文是在导师吴良才教授的热忱关怀和细心指导下完成的。从论文选题到最后的成文，无不凝结着导师的心血。导师开拓创新的精神、严谨踏实的治学态度以及宽厚的为人都使我深受熏陶，并将使我终身受益。在此谨向导师致以衷心的感谢和诚挚的敬意。此外，也感谢测绘专业的全体老师，正因为有了他们的虚心教导，使我的知识得到一点一滴的积累，才从容的完成今天的毕业设计。30 东华理工大学毕业设计参考文献参考文献[1] 崔明理. 控制测量手册. 山西: 山西科学技术出版社, 1999.1. [2] 吴俊昶. 控制网测量平差 2 版. 北京: 北京测绘出版社, 1998.1. [3] 谭辉. 测量程序与新型全站仪的应用. 北京: 机械工业出版社, 2006. [4] 孔祥元, 梅是义. 控制测量学下册2 版. 武汉: 武汉大学出版社, 2002.2. [5] 武汉大学测绘学院测量平差学科组. 误差理论与测量平差基础. 武汉: 武汉大学出版社, 2003.1. [6] 张智星. MATLAB 程序设计与应用. 北京: 清华大学出版社, 2002.4. [7] 郭九训. 控制网平差程序设计. 北京: 原子能出版社, 2004.8. [8] 姚连壁, 周小平. 基于MATLAB 的控制网平差程序设计. 上海: 同济大学出版社, 2006.6. [9] 清源计算机工作室. MATLAB 基础及其应用[M]. 北京: 机械工业出版社, 2000. [10] 程卫国,冯峰, 姚东, 徐昕, 等. MATLAB 5.3 应用指南[M ]. 北京: 人民邮电出版社, 1999. [11] CHAI Yan2jun , et al . Quasi2accurate Detection Methodsof Gross Errors with MATLAB Language [J ] . Journal of Shandong University of Science and Technology (NatureScience) , 2000 , (3) [12] NIE Gui2gen. Application of MATLAB to Surveying Da2ta Processing [ J ] . Bulletin of Surveying and Mapping ,2001 , (2) 31 东华理工大学毕业设计附录一附录一水准网平差程序代码：function level3 [ed,dd,sd,gd,pn,h0,k1,k2,h1,s]=readlevelnetdata; global pathname net_name s_datafile
a1_datafile; global ed dd sd pn gd h0 k1 k2 h1 s dh; [dh,h,V,L,uw0,uwh,uw1]=calculatelevelnet(ed,dd,sd,pn,gd,h0,k1,k2,h1,s);
writelevelnetdata(pn,k1,k2,h1,V,L,h0,dh,h,uwh,uw0); return function[dh,h,V,L,uw0,uwh,uw1]=calculatelevelnet(ed,dd,sd,pn,h0,be,en,hd, distance); % 水准平差网A=sparse(zeros(sd,gd)); b=(0:(gd-1)*sd); A(be'+b)=-1; A(en'+b)=1; A=A'; A=A(:,1:dd); l=zeros(gd,1); %求解常数项%权阵%高程改正数%待定点高程近似值%待定点高程平差值%高差观测值改正数%高差观测值平差值l=h0(be)-h0(en)+hd; p=diag(1./distance); dh=inv(A'*p*A)*A'*p*1; h00=h0(dd+1:sd); h0=h0(1:dd); h=h0+dh; V=A*dh-1; L=hd+V; %精度评定uw0=sqrt(V'*p*V/(gd-dd)); Qxx=inv(A'*p*A); uwh=uw0*sqrt(diag(Qxx)); uwh(dd+1:ed+dd)=0.0; 32 %输出水准网计算结果%求解系数阵%单位权中误差%待定点高程平差值中误差东华理工大学毕业设计附录一Qff=A*Qxx*A'; uwh=uw0*sqrt(diag(Qff)); h=[h;h00]; h0=[h0;h00]; dh=[dh;zeros(ed,1)]; return function [ed,dd,sd,gd,pn,h0,k1,k2,h1,s]=readlevelnetdata; global pathname net_name s_datafile b_datafile; global ed dd sd pn gd h0 k1 k2 h1 s k11 k12; k1=[];k2=[];h=[];s=[]; if(isempty(pathname)|isempty(net_name)) [filename,pathname]=uigetfile('*.txt','Input filename'); i=find('.'==filename); net_name=filename(1:i-1); end fid1=fopen(strcat(pathname,net_name,s_datafile),'rt'); if(fid1==-1) msgbox('Input File or Path is not correct','Warning','warn'); return; end %open afile to read ed=fscanf(fid1,'%f',1); dd=fscanf(fid1,'%f',1); sd=ed+dd; gd=fscanf(fid1,'%f',1); pn=fscanf(fid1,'%f',sd); %known data h0=fscanf(fid1,'%f',ed); h0(dd+1:ed+dd)=h0(1:ed) heightdiff=fscanf(fid1,'%f',[4,gd]); heightdiff=heightdiff'; k1=heightdiff(:,1); k2=heightdiff(:,2); %起点%终点%已知点高程%已知点个数%未知点个数%总点数%观测点个数%点号%open afile to read %高差平差值中误差%所有点高程33 东华理工大学毕业设计附录一k11=heightdiff(:,1); k12=heightdiff(:,2); h1=heightdiff(:,3); s=heightdiff(:,4); fclose('all'); %点号转换[k1,k01]=chkdat(sd,pn,k1); [k2,k02]=chkdat(sd,pn,k2); h0(1:dd)=20000.; ie=0; while(1)%计算近似高程for k=1:gd i=k1(k); j=k2(k); %起点%终点%高差%距离if(h0(i) [a]=fscanf(fit1,'%f',2*ed); for i=1:ed x0(i)=a(2*i-1); y0(i)=a(2*i); fprintf(fit2,'%15.3f %15.3f\n',x0(i),y0(i)); end fprintf(fit2,'\n'); m1=fscanf(fit1,'%d',1); fprintf(fit2,'%5d\n',m1); isid=0; if(m1>0) fprintf(fit2,'side\n'); [a]=fscanf(fit1,'%f',2); ms=a(1); pp=a(2); fprintf(fit2,'%6.2f %6.2f\n',ms,pp); [a]=fscanf(fit1,'%d %d %f',3*m1); for i=1:m1 e(i)=a(3*i-2);d(i)=a(3*i-1);sid(i)=a(3*i); fprintf(fit2,'%5d %5d %15.3f\n',e(i),d(i),sid(i)); end [e,i1]=chkdat(sd,pn,e); [d,i2]=chkdat(sd,pn,d); i3=0; for i=1:m1; if(e(i)==d(i)) i3=1; fprintf(fit2,'* * %5d %5d %5d * *\n',i,e(i),d(i)); end fprintf(fit2,'%5d %5d %15.3f\n',e(i),d(i),sid(i)); end 37 东华理工大学毕业设计附录二isid=i1+i2+i3; end idir=0; kk=isid if(kk>0) msgbox('Error by function rddat1','Warning','warn'); return; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%近似坐标计算xyknow(1:ed)=pn1(1:ed);xyunknow(1:dd)=pn1(ed+1:sd); [aa]=fscanf(fit1,'%d',2); A=aa(1);B=aa(2); mn=fscanf(fit1,'%d',1); [bb]=fscanf(fit1,'%d',mn); m=bb; for i=1:length(m) if i==1 if e(m(i))==A|e(m(i))==B P=d(m(i)); else P=e(m(i)); end else if e(m(i-1))==A|d(m(i-1))==A B=P; else A=P; end if e(m(i))==A|e(m(i))==B P=d(m(i)); else P=e(m(i)); end 38 东华理工大学毕业设计附录二end for j=1:m1 if (e(j)==P|d(j)==P)&(d(j)==A|e(j)==A) s1=sid(j); end if (e(j)==P|d(j)==P)&(d(j)==B|e(j)==B) s2=sid(j); end end deltx=x0(B)-x0(A); delty=y0(B)-y0(A); alfaAB=alfa(deltx,delty) ; ss=sqrt(deltx*deltx+delty*delty); ee=(s1*s1+ss*ss-s2*s2)/(2*ss);ff=sqrt(s1*s1-ee*ee); x0(P)=x0(A)+ee*cos(alfaAB)-ff*sin(alfaAB); y0(P)=y0(A)+ee*sin(alfaAB)+ff*cos(alfaAB); end fprintf(fit2,'\n'); fprintf(fit2,'计算的近似坐标\n'); for i=1:sd fprintf(fit2,'%15.3f %15.3f \n',x0(i),y0(i)); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% k=1; while(k) lo=2062.648062470964; n=2*dd; sum=n*(n+1)/2.0; sd=ed+dd; a=[]; 39 东华理工大学毕业设计附录二a(1:m1,1:9)=0.0; c(1:sum)=0.0; w(1:n)=0.0; for i=1:m1 %边长观测误差方程dx=x0(d(i))-x0(e(i)); dy=y0(d(i))-y0(e(i)); ss=sqrt(dx*dx+dy*dy); cosa=dx/ss; sina=dy/ss; a(i,1)=2*e(i)-1-2*ed+1.0e-9; a(i,2)=-cosa; a(i,3)=a(i,1)+1; a(i,4)=-sina; a(i,5)=2*d(i)-1-2*ed+1.0e-9; a(i,6)=cosa; a(i,7)=a(i,5)+1; a(i,8)=sina; a(i,9)=10.0*(ss-sid(i)); %ql(i)=(ms^2+(ss*pp*0.0001)^2); ql(i)=((ms+ss*pp*0.001)); ql(i)=100/ql(i)^2; fprintf(fit2,'%15.3f %15.3f %15.3f %15.3f %15.3f %15.3f %15.3f %15.3f %15.3f \n',a(i,:)); end % fprintf(fit2,'%15.3f for i=1:m1 for j=1:4 jj=fix(a(i,2*j-1)); if(jj<=0) continue; end w(jj)=w(jj)+a(i,2*j)*a(i,9)*ql(i); di=(jj-1)*(n-jj/2.0); for k=1:4 kk=fix(a(i,2*k-1)); \n',ql); %形成法方程40 东华理工大学毕业设计附录二if(kk<=0|jj>kk) continue; end c(di+kk)=c(di+kk)+a(i,2*k)*a(i,2*j)*ql(i); end end end fprintf(fit2,' \n'); fprintf(fit2,'%15.3f \n',w); fprintf(fit2,'%15.3f %15.3f %15.3f %15.3f %15.3f %15.3f %15.3f %15.3f %15.3f \n',c) %%%%%%%%%%%%%%坐标改正数和单位权中误差的计算sd=ed+dd; n=2*dd; k=0; for i=1:n dxy(i)=0.0; di=(i-1)*(n-i/2.0); for j=1:n dj=(j-1)*(n-j/2.0); if(j n=2.0*dd; maxmm=0.0; smm=0.0; for i=1:dd ii=ed+i; di=(2*i-2)*(n-(2*i-1)/2.0); dj=(2*i-1)*(n-i); ql=c(di+2*i-1); q2=c(dj+2*i); q3=c(di+2*i); dl=sqrt(abs(ql+q2-q3)); d11=uw0*dl; xx=ql-q2; yy=2*q3; zz=ql+q2; mx1=sqrt(ql); my1=sqrt(q2); mm1=sqrt(zz); if(abs(xx)<1e-10) fi(i)=sign(90.0,q3); else fi(i)=atan(yy/xx)*57.2958; end if(xx>=0.0&yy>=0.0) fi(i)=fi(i)/2.0; elseif(xx>=0.0&yy<=0.0) fi(i)=(fi(i)+360.0)/2.0; fi\n'); ai bi fi parameter of error ellipse\n'); pn(i) mx my mm a %fit2=fopen('D:\大牛文档\work\siliout.txt','wt') 43 东华理工大学毕业设计附录二elseif(xx<0.0) fi(i)=(fi(i)+180.0)/2.0; end ww=sqrt(xx*xx+yy*yy); a1=sqrt((zz+ww)/2.0); b1=sqrt((zz-ww)/2.0); ab1=a1-b1; mx=uw0*mx1; my=uw0*my1; mm=uw0*mm1; if(mm>maxmm) i1=ii; end smm=smm+mm; ai(i)=uw0*a1; bi(i)=uw0*b1; ab=ai(i)-bi(i); fprintf(fit2,'%10d my,mm,ai(i),bi(i),fi(i)); end smm=smm/dd; fprintf(fit2,' %fprintf(fit2,' pp= fprintf(fit2,' mse of unit weight= %9.6f\n',uw0); the maximum station error mm= %8.3f(cm) the average station error mm= %8.3f\n',smm); %10.3f %10.3f %10.3f %10.3f %10.3f %10.3f\n',pn(ii),mx, %4d\n',maxmm,p n(i1)); fprintf(fit2,' %14.4f\n',uw0); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%% %%%%%%%%%%%%%%绘制控制网图kk=0; %tile1={'误差椭圆比例'}; %prompt011={'请输入误差椭圆的比例'}; lines=1; d11={'100'}; c12=inputdlg('plese input the scale','scale',lines,d11);%输入误差比例椭圆对话框44 东华理工大学毕业设计附录二scale=str2num(str2mat(c12)); %scale=500; hold off title('控制网网图'); xlabel('Y 坐标'); ylabel('X 坐标'); plot(y(1:ed),x(1:ed),'k') hold on plot(y((ed+1):end),x((ed+1):end),'k.') for i=1:ed+dd %写点名text(y(i)+3,x(i)+3,num2str(pn(i))); end for i=1:m1%画边长连线plot([y(e(i)) y(d(i))],[x(e(i)) x(d(i))],'r'); end for i=1:dd %画误差椭圆fe(i)=fi(i)*pi/180.; plot(y(i+ed)+(ai(i)*cos(0:0.01:2*pi)*sin(fe(i))+bi(i)*sin(0:0.01:2*pi)*cos(fe(i)))*scale,... x(i+ed)+(ai(i)*cos(0:0.01:2*pi)*cos(fe(i))-bi(i)*sin(0:0.01:2*pi)*sin(fe(i)))*scale,'r'); end axis equal; zoom on fclose('all'); open(strcat(pathname,net_name,b_datafile)); return %%%%%%%%%%%%%%%%%%%%55%%%%%%%%%%%%%%%%%% function a=alfa(dx,dy)%%%%%计算方位角函数a11=atan(dy/dx); if (dy>0&dx>0) a=a11; elseif (dx<0) a=pi+a11; else a=2*pi+a11; 45 东华理工大学毕业设计附录二end 46 东华理工大学毕业设计附录三附录三法方程求逆函数function c=invsqr(c,n) ss=0.0; for i=1:n di=(i-1)*(n-i/2.0); for j=i:n ss=c(di+j); for k=1:i-1 dk=(k-1)*(n-k/2.0); ss=ss-c(dk+j)/c(dk+k); end if(i==j) c(di+j)=1/ss; else c(di+j)=ss*c(di+i); end end end for i=1:n-1 di=(i-1)*(n-i/2.0); for j=i+1:n ss=-c(di+j); for k=i+1:j-1 dk=(k-1)*(n-k/2.0); ss=ss-c(di+k)*c(dk+j); end c(di+j)=ss; end end for i=1:n-1 di=(i-1)*(n-i/2.0); for j=i:n dj=(j-1)*(n-j/2.0); 47 东华理工大学毕业设计附录三if(j==i) ss=c(di+j); else ss=c(di+j)*c(dj+i); end for k=j+1:n dk=(k-1)*(n-k/2.0); ss=ss+c(di+k)*c(dj+k)*c(dk+k); end c(di+j)=ss; end end return 48
相关主题
文本预览