基于MATLAB的测量平差数据处理讲解

1、基于MATLAB的测量平差数据处理摘要MATLAB是目前在研究机构广泛应用的一种数值计算及图形工具软件，它的特点是语法结构简明、数值计算高效、图形功能完备，特别适合非专业编程员完成数值计算、科学试验处理等任务。以往的测量数据处理方法需要编制特定的处理矩阵运算程序，而且程度复杂，难度大。本文介绍一种基于MATLAB的水准网和测边网的程序设计方法，与其它算法语言相比，具有编程简单，运算速度快的特点。文中分别阐述了水准网和测边网程序的理论基础、实现步骤和运行结果。通过实例的分析，总结出利用MATLAB对测量数据处理有很大的应用价值，它缩短了编程的时间，提高工作效率。关键词：MATLAB；平差；程序设

2、计ABSTRACTMATLABisonespeciesofnumerical-valuescalculationandgraphictoolssoftwarewhichiswidelyusedtoapplyatresearchinstitutionsatpresent.Theparticularitiesare:concisegrammar-structure、highlyefficientinnumericalvaluescalculating、completefunctionofgraphs、especiallyitisadaptedtoevildoingprofessionalprogr

3、ammertoaccomplishthetasksthatarenumerical-valuescalculatingandscientificexperimentstreating.Theancientmethodsofmeasureddata-processingneedestablishingspecialproceedingsoftreatingmatricesoperation,moreover,itiscomplexandgreatlydifficult.Thisarticleintroducesoneprogrammingmethoddealingwithlevelingandm

4、easuringedgenetworkbasedonMATLAB.Comparedwithotheralgorithmlanguage,ithasparticularitieswhicharesimplyprogrammingandquicklyoperating.Thearticleseparatelyexpatiatethetheoriesbasics、realizingstepsandrunningresultsatlevelingandmeasuringedgenetwork.Withtheanalysisofexamples,ithasprodigiousapplicationval

5、ueinmeasureddata-processingbyuseofMATLAB.Moreover,itshortensprogrammingtimeandimprovesworkingeffectiveness.Keywords：MATLAB；programming绪论作为一名测量技术人员，如果不掌握一门PC机编程语言与便携计算工具，要想提高测量工作的效率几乎寸步难行。测量需求的多样性与复杂性，造就了测量计算鲜明的个性化特点，这就是在商业测量计算软件高度发达的今天，掌握一种实用的程序语言进行编程计算仍有广泛的市场需求的重要原因。当今较流行的计算机程序语言基本上都是基于Windows的，例如T

6、urboPascal,VisualBasic,VisualC,BorlandC+等，这些程序语言的优势是基于对象及可利用Windows丰富的系统资源，应用它们可以开发出界面非常丰富和友好的应用程序，其劣势主要有以下几点：Windows程序都非常庞大，学习并熟练掌握它们并非易事。虽然市场上已有的多种专用的测量平差软件都是采用C语言开发的，但这些软件价格都比较贵，而且都带有加密狗，一次只能供一个用户使用。出于商业目的，开发商不会公开程序源代码，这为修改程序功能以适应用户的特殊需求带来了不便。在测量生产中，经常需要根据工程的实际情况进行一些个性化的数值计算工作，这些数值计算工作无固定模式，这就需要求

7、测量技术人员最好能熟练掌握一种适用于数值计算的程序语言，以便提高测量计算的效率。C语言的数值计算语句不够丰富，例如，在测量平差计算中，经常需要进行的矩阵运算，尤其是解法方程的矩阵求逆不能直接使用语句实现，而必须应用计算机算法编程实现。如果不是基于商业软件开发，只为满足实际测量工作计算需要，则C语言的劣势就变成了MATLAB语言的优势。内蒙古科技大学毕业设计MATLAB软件简介1.MATLAB软件简介MATLAB是从Matrix（矩阵）和Laboratory（实验室）各取前3个字母组成的,意思是矩阵实验室,是美国MathWorks公司于20世纪80年代中期推出的一种交互式、面向对象的科技应用软件

8、，是一个为科学和工程计算而专门设计的高级交互式软件包。MATLAB集成了图示与精确的数值计算,是一个可以完成各种计算和数据可视化的强有力工具,其优秀的数值计算能力和卓越的数据可视化能力使其很快在数学软件中脱颖而出,成为以矩阵运算为主要工作方式的线性代数、概率论和数理统计、自动控制、数字信号处理、动态系统仿真等领域教学和科研工作者的有力武器。随着该软件自身的发展及市场的需求,其功能日趋完善,其最高版本7.0版已经推出,随着版本的不断升级,它的数值计算及符号计算功能得到了进一步完善。MATLAB是以矩阵作为数据操作的基本单位，矩阵的生成、运算、转置、求逆等非常简单。在MATLAB环境中，不需要对创

9、建的变量对象给出类型说明和维数，所有的变量都作为双精度数来分配内存空间,MATLAB将自动地为每一个变量分配内存。MATLAB语言起源于矩阵运算，并已经发展成为一种高度集成的计算机语言，它提供了强大的科学运算、灵活的程序设计流程、高质量的图形可视化与界面设计、便捷的与其他程序和语言接口的功能oMATLAB系统主要包含5部分的内容:MATLAB工作环境、Mablab数学函数库、MATLAB语言体系、句柄图形、MATLAB应用程序接口（API）。MATLAB系统主要功能包括:数值计算功能、符号计算功能、数据分析和可视化、文字处理功能、SIMULINK动态仿真功能。同时,MATLAB又是开放的，除了

10、内部函数之外，所有的MATLAB主包文件和各工具包文件都是可读可改的源文件，用户可以作为参考掌握其用法，并可对其修改以适应自己的需要，也可加入自己编写的文件构成新的工具包。例如，随着GPS的广泛应用,OrionDynamicsandCon21rolCorporation、ConstellInc.GPSSoftLLC、NavsysCorporation等多家公司都相应开发出了适于GPS数据处理的MATLAB工具箱。MATLAB是一个集数值计算、图形管理、程序开发于一体的功能十分强大的系统。将MATLAB应用于测量数据的处理是一件非常有意义的工作Mo2hamed等曾成功地在MATLAB系统中利用白

11、滤波技术研究动态解算GPS载波相位信号的模糊度问题。因为测量数据的处理特别是测量平差主要应用矩阵运算，而MATLAB又特别易于做矩阵运算，因此，研究开发基于MATLAB的测量平差方法具有极好的应用价值。毕业设计水准网平差程序2MATLAB在测量平差中的应用测量平差数据处理主要是基于矩阵的运算,常用的矩阵运算主要是矩阵的生成、转置、求逆和矩阵求广义逆等。在MATLAB环境中，不需要对创建的变量对象给出类型说明和维数，所有的变量都作为MATLAB中的M文件的语法与其他的高级语言类似，是一种程序化的编程语言，同时也是一种解释性的编程语言，即逐行解释运行程序，使程序容易调试，计算更为简捷,而且对于平差

12、原理理解和掌握变得更容易。另外,MATLAB语言与数学语言比较接近，更容易掌握和理解。2.1测量平差原理的概述测量平差的函数模型有条件方程和观测方程。以条件方程为函数的模型的最小二乘平差称为条件平差；在条件方程中，根据需要如果还设有一定数量的未知数，则称为附有参数的条件平差；以观测方程为函数模型的最小二乘平差称为间接平差；如果观测方程中的某些参数不独立，则这些不独立参数必然存在一些条件，称这种平差模型为附有条件的间接平差。本文的两个程序都采用间接平差模型。对于一个实际平差问题，根据所选参数的个数、选什么量为参数以及参数之间是否函数独立，经过仔细推敲可以发现附有条件的间接平差模型本身就是各种经典

13、平差模型的概括模型，其余的经典平差模型，如条件平差模型、间接平差模型、附有未知数的条件平差模型和附有限制条件的条件平差模型都是它的特例。间接平差的公式汇集：间接平差模型为V二BX-l21）（2-2）（2-3）VTPVTmin系数矩阵B满秩，即rank（B）=t法方程及解为：Nx-f二0,（N-i二BTPB,f二BtPI）bbebbex=N-ifbbe参数的平差值：X二X+x0（2-4）观测量的平差值：L二L+V（2-5）单位权的中误差：花占(2-6)平差参数的协方差阵：D=b2N-1(2-7)XX0bb平差函数的协方差阵：Q2FtN-iF(2-8)榊0bb2.2平差程序总体方案MATLAB号称

14、为全球工程师的共同语言，其语法和C语言相似，但它有强大的数值计算和绘图功能，这使之在工程应用方面的计算更出色，本文就基于这种程序设计语言环境设计一个控制网平差程序。该程序包含了一个高程控制网平差程序和测边网平差程序。本程序适用于各种等级的高程网和测边网，程序在设计过程中，始终考虑数据的储存量。因而本程序不储存误差方程的系数和常数项，对待定点数较多的平差网，组成法方程的系数矩阵是个稀疏矩阵，如待定点的编号恰当，法方程的系数会集中在主元系数的两侧形成带状。为减少法方程系数的储存量，只要按行储存下三角阵或按列储存上三角阵中第一个非零系数起的系数，就是通常叫做维变带宽储存方法。3.水准网平差程序3.1

15、程序的功能本程序适用于二、三、四等水准网平差计算，平差的水准网可以是独立的、也可以是附合网，其主要功能是完成水准网的平差计算和精度评定计算。平差计算采用间接平差法，以归心的观测值为高差，以未知点高程为未知参数。精度评定计算包括计算单位权中误差和每个待定点的高程中误差。3.2水准模型网的间接平差3.2.1“权”值的确定当在相同的条件下进行水准测量时，其精度是相同的，因而观测结果的可靠性也是同样的。但如果在不同的条件下进行水准测量时，高程的精度就有所不同，此时称为不等精度观测，所求出的未知量的值、高程的最或是值并对其精度进行评定时，就需要“权”了。由于观测的不等精度，因而观测值的可靠程度不同，求未

16、知量的最或是值时，这样的一个因素就必须考虑了，这个因素是:可靠性大的某观测值，其精度高，对测量的最后结果的影响也就越大。此时用“权”值来表示观测值的可靠程度那么，“权”值愈尢观测值的可靠程度就愈高。另外，在观测过程中,观测值的中误差愈小，观测结果愈可靠,它的“权”值就愈大。因而根据中误差来确定“权”值是非常适当的。设以Pi表示观测值Li的“权”,m为中误差，则“权”值的定义为：P.=吕(3-1)式中:A为任意的正常数，在一组观测值中为一个定数。在实际测量中，通常是观测值的中误差事先并不知道，因而必须先确定观测值的“权”,然后才能求出未知量的最或是值。此时可以利用距离(S)或测站数(N)来确定观

17、测值一高程的“权”。根据偶然误差传播定律,各观测点高程Hi的中误差mi由测站数Ni确定时，则有：(3-2)式中:m为一组观测值的中误差,为一个定数.由(3-2T)、(3-2-2)两式可得：P_C(3-3)式中上二亠.由于卫5邯为定数故C也为定数.7U同样可得出:(3-4)式中:C为定数，s为测距.i由(3-3)、(3-4)两式可以得出这样一个结论:当测站观测高差等精度时，观测总高差的“权”与测站数或距离成反比。3.2.2水准路线的平差计算1.附合路线的平差计算假定在图1示的A、B两水准点之间布设一条水准路线,A、B两水准点的高程为已知，分别设为H、H、n、n、C为中间水准点。假定观测了所有的点

18、的高TOC o 1-5 h zAB12程,现拟求C点的高程H的最或是值。CH可由水准路线AfC、BfC分别观测的高差Ah、Ah计算得出，由此而CACBC得到的观测高程分别设为Hcl、Hc2，其值为：Hc1=H+Ah;Hc2=H+AhAACBBC当Hcl、Hc2在不等精度条件下观测得出时,它们的“权”也不同，分别设为Pel、Pc2,这样C点的高程H的最或是值为：C(3-5)”PHPHTOC o 1-5 h zH=clclc21cP+Pclc2根据A点的高程H,AfC水准路线观测的高差Ah以及BfC水准路线观测AAC的高差Ah，可推算出B点的观测高程H为：BCBH=H+Ah-AhBAACBC水准路

19、线AfB的高程闭合差为：f=H-H=HH（3-6）hBBclc2由（3-6）式得到：H=H-fc2clh HYPERLINK l bookmark82 o Current Document CC由（3-3）式得到：P=亠P=亠（N、N分别表示水准路线AfC、BclN、c2NACBCfC的测站数，水准路线AfB的测站数N二N+N）ABACBC将上述表达式代入（3-2-5）式中,得到：(3-7)如果以水准路线AfC的距离S、BfC的距离S、AfB的距离SBCBCAB(S二S+S)ABACBC来确定高程观测值的“权”值时，同样可以得到:H.图3-1水准路线图2闭合路线的平差计算闭合路线的平差计算原理

20、与附合路线相同，因而(3-7)、(3-8)两式的结论适用于闭合路线的平差计算。(3)具有一个结点的水准网的平差计算工PHAiAiXpAii=1(3-9)如图2所示为具有一个结点的水准网,B,C,D,为已知高程水准点,B-A,C-A,D-A,为水准路线，则接点A的高程最或是值为：口PH+PH+PH+H二一A2A2A3A3AP+P+P+A1A2A3式中H,H,H分别为水准路线B-A,C-A,D-A,计算A的观测高程，各高A1A2A3程相应的“权”值为P,P,PA1A2A3设H,H,H的算术平均值为H0，各高程观测值与H0的差值分别为SA1,A1A2A3ASA2,5A3,则有：41=耳-找Jtl=戌

21、+越1狂=验-戌或毘2=&2&2=屁_出H石=止+(3-10)将(3-10)式代入(3-9)式得到:(3-11)当以测站数和距离来确定“权”值时，（3-11）分别可以转化为：耳=疋+于当，（3-12）耳=此十占；（3-13）上述结论也可应用于小三角水准网平差计算。3.2.3精度评定单位权中误差：bVTPV0Vnt(3-12)平差参数的协方差阵：D-b2N-1XX0bb(3-13)平差函数的协方差阵：Q-b2FtN-1F林0bb(3-14)3.3水准网间接平差程序信息设计1.数据文件的组织下面给出一个水准网输入数据文件的例子：336（已知点个数、未知点个数、观测值个数）1011021031041

22、05106（点号）34.78835.25937.825（已知点高程）1011.6254.5（起点点号、终点点号、高差观测值、距离观测值）101102-0.4183.11020.7143.41031.2433.8106103-0.5774.3101-0.7862.5（其中编号数组未知点在前，已知点在后）2.水准网平差变量约定表3-1变量约定表变量名说明ed已知点个数dd未知点个数sd总点数gd观测值个数pn点号h0已知点高程be起点点号en终占占号k、八、八、Jh1咼差观测值s距离观测值3.4水准网程序与使用说明3.4.1水准网程序流程图图3-2水准网流程图程序的全部代码见附录一。3.4.2水准

23、网程序的使用本程序使用MATLAB的矩阵功能计算法方程，在运行程序前首先要有其始数据。其始数据是一文件的形式保存在磁盘中，文件的格式在上文已经说明过,编好文件后,以后缀名为.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.&s5=4.3,

24、s6=2.5.图3-3水准网图首先编数据文件，命名为data1.txt.数据的格式如下：33610210310410510634.78835.25937.8251011.6524.5102-0.4183.11020.7143.41031.2433.8103-0.5774.3101-0.7862.5进入MATLAB界面，在命令窗口直接输入level3运行程序。弹出如下窗口Inputfilename查我范團Q：口work_T十匡呼datal.如七data3.txt0新建文本文裆1.3Jout0新建丈本文挡。毗U1打开|*.txtJ取消文件名逼：文件类型(I):图3-4数据读入文件选择datal.t

25、xt即可运行出如下结果：FileEdi+TeKtCellToolsDebugDesktopWindowHelp站a#4A4*VSStack：Base田待定点高程值及中误差:点号近似高程咼差高程平差値高程平差値中误差pnh0(m)dh(m)h(m)vwh10136.4400-0.008036.4320C.016210236.0220-0.029535.9925C.015910337.2650-0.034137.2309C.016310434.78800.000034.7880C.000010535.25900.000035.25900.000010637.82500.000037.82500.0

26、000高差观测値平差値:plainfileLn1Col1DV1pnpnh0(m)V(m)H(m)1041011.6520-0.00801.6440101102-0.4180-0.0215-0.43951051020.71400.01950.73351021031.2430-0.00461.238106103-0.5770-0.0171-0.5941103101-0.7860-0.0129-0.7989D12345678910111213141516171819单位权中误差：0.0118m图3-5计算结果在图3-5中，分别输出了高程的平差值及精度。结果是一文本的形式保存，用户可对它进行编辑。水准

27、网平差程序代码functionlevel3ed,dd,sd,gd,pn,h0,k1,k2,h1,s=readlevelnetdata;globalpathnamenet_names_datafilea1_datafile;globaledddsdpngdh0k1k2h1sdh;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);%输出水准网计算结果returnfunctiondh,h,V,L,uw0,uwh,

28、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+

29、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;Qff=A*Qxx*A;uwh=uw0*sqrt(diag(Qff);%高差平差值中误差h=h;h00;%所有点高程h0=h0;h00;dh=dh;zeros(ed,1);returnfunctioned,dd,sd,gd,pn,h0,k1,k2,h1,s=readlevelnetdata;globalpathnamenet_names_datafi

30、leb_datafile;globaledddsdpngdh0k1k2h1sk11k12;k1=;k2=;h=;s=;if(isempty(pathname)|isempty(net_name)filename,pathname=uigetfile(*.txt,Inputfilename);i=find(.=filename);net_name=filename(1:i-1);endfid1=fopen(strcat(pathname,net_name,s_datafile),rt);if(fid1=-1)msgbox(InputFileorPathisnotcorrect,Warning,w

31、arn);return;end%openafiletoread%openafiletoreaded=fscanf(fid1,%f,1);%已知点个数dd=fscanf(fid1,%f,1);%未知点个数sd=ed+dd;%总点数gd=fscanf(fid1,%f,1);%观测点个数pn=fscanf(fid1,%f,sd);%点号%knowndatah0=fscanf(fid1,%f,ed);%已知点高程h0(dd+1:ed+dd)=h0(1:ed)heightdiff=fscanf(fid1,%f,4,gd);heightdiff=heightdiff;k1=heightdiff(:,1);

32、%起点k2=heightdiff(:,2);%终点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(l)%计算近似高程fork=1:gdi=k1(k);j=k2(k);if(h0(i)le4)h0(j)=h0(i)+h1(k);ie=ie+1;endif(h0(i)le4&h0(j)A（

33、m,9）其中，m为观测值个数，n为未知点个数的两倍。改进后的A阵格式为A=（编号1,系数1,编号2,系数2,编号4,系数4,常数项）i共9列。即只存储误差方程的4个非零参数系数。法方程系数阵N为对称阵,在存储时,只需要存其上三角部分就可以了。其占A用的空间为：n（n+1）sum=2现有A阵：A=（编号1,系数1,编号2,系数2,编号4,系数4,常数项）其中偶数项为系数，加上最后的A9为常数项，在组成法方程时，从A2开始分别与剩下的偶数项以及常数项相乘,然后再用A4与剩余的项相乘,一直到A8为止,这样就完成了N=AtPA的过程。需要注意的是：若A1,A3,A5,A7小于零，则表A示该点已知点,不

34、参与法方程的组成。4.1.2边长观测的权边长观测的精度一般与其长度有关,定权公式为C2P二0-siQ2si（i二12，n）式中O2为所测边长s的方差，a2为任意选定的单位权方差。si0i为了定权P必须已知测边的先验方差a2sisi但精确的已知是十分困难的,一般采用厂方给定的测距仪精度,即a=a+bSsii式中，a为固定误差（单位mm）,b为比例误差（单位：ppm）,S为边长（单i位km)。4.1.3解算法方程由于法方程是对称正定阵，因此，可采用改进的平方根法进行解算。平方根法是对称正定矩阵非常有效的三角分解方法，设A为n阶方阵，如果其所有顺序主子式均不为零，则其存在唯一的分解式：A=LDR其中

35、L=ri1.2，D=rdi0、，R=rirr12.in1ni11n,n-1丿10d丿nrn-i,n1丿由于此住A对称性，得L=RT，又根据A阵正定的性质，可证明D均为正数。现在设(ddiD=d即111丄，A=LDLt=LD2D2Lt=(LD2)(D2Lt)=LLT则为方便A=LLt称为Cholesky分解，即正定对称矩阵的平方根分解法。解AX=b等阶于求解两个三角方程组：LY=b和LtX=Y在用平方根分解法计算时，需要进行n次开方运算。为了避免开方，可以直接采用对称正定的A=LDL分解式对平方根法进行改进。从而解方程组AX=b可以按如下步骤进行：把A分解成A=LDLt，则AX=b变成(LDLT

36、)X二b，即等价于LY=bLtX=D-iY由此可以解出X和Y。这称为改进的平方根法，在计算中避免了开方运算。平方根法和改进的平方根法的计算量和存储量比消去法节约近一半，而且不需要选主元，能得到比较精确的数值解。法方程用改进平方根法解算的过程如下(1)分解：C=StD-is其中1S11S1n,D=1S11kS丿mnkS丿mn1c11c)1ns-2-s111S11S)1ncn1C丿nnS丿mns-ns11n1,72sn-1,n-1s1j=c1jssc2j12_Ljsc3js11ss131j2js2jc2js11s11s+232j+ss223j,s3j=c3js131js1122cijs1jcij-

37、1k=1ssskki1,j2纯量计算公式为i=12)求逆R=S-11r11r、1nR=r丿nn由RS=I得1r11r12r22r1nr2nsiis12s22rnn227纯量计算公式：r=一iisiirsr=ii1212s22(rs+rs)r二1113122313s33(rs+.+rsr=111n1(n1)(n1)n1nsnnrsr=221223s33(rs+rs)r=22_24233424s44(rs+rsr=222n2(n1)(n1)n2nsnn通式为1r=iisii玄rsikkjr二一i,jiijsij3）求积Q=(StD-1S)-1=S-1D(S-1)t=RDRtQ=RDRt厂sr111

38、1s12sr11r1n(s11rnn八snnr11r1nr12r22sr1n1nsr2n2nr11rnnr丿nn、rnnsrnnnnr1nrnnrnn4.1.4精度评定（1）坐标改正数以及单位权中误差的计算m0使用上三角一维数组形式存储坐标改正数的公式为：5x=一可qw-工qw,i=1,nijijjjj=1j=1其中，n=2xdd,x的单位是cm.i平差值：X=Xo+x写成分量的形式，为X=X0+5xiii如果近似坐标的误差较大，或网形较大，平差的结果不会精确，这时，就需要进行迭代平差，直到两次平差间互差在允许值内。由测量平差理论：=Vtpv.n-1同样可得到单位权中误差：fiPVVTm=v-

39、0m一n其中，m-n观测个数减去未知点个数；m=m+m+m,n=2xdd+ST,ST一方向观测的测站数123PVV=pll+5x（2）点位误差椭圆误差椭圆表示了网中点或点与点之间的误差分布情况如图。在测量工作中，常用的误差椭圆对布网方案作精度分析。绘制误差椭圆只需要三个数据：椭圆长半轴a，短半轴b和主轴方向，其求法为2otan2p=xyo2o2xya2=(O2+O2+2o2)2+402)2xyxyxyb2=丄(02+o2一(o2-o2)2+4o2)2xy耳xyxy1793.00266225.85440049.22953782-79036924.72861027.S8613760.7061055

40、187.342065483.1587838.S80061051017493.3237494.8810304108058884.5：871QS081021087228.36721012图4-7数据文件2.在MATLAB命令窗口键入nnbb执行程序，运行中会弹出一个对话匡提示用户输入误差椭圆的放大比例，默认为100,本例选择500。如下图图4-8a对话框图4-8b对话框3计算出的结果如下：在图4-9a中，第1、2行分别是已知点和未知的个数；第4行是点的编号；第5至第8行是已知点的坐标；第10行是观测值的个数；第12行是测距的固定误差和比例误差；13至25行是观测边的起点号、终点号和观测边长；27至

41、39行（接到图4-9b）是点号转换为计算机顺序后的观测边起点号、终点号和观测边长。接图4-9a，在图4-9b中，第42至49行是推算的近似坐标；第52至64行是计算的误差方程系数和常数。在图4-9中，第67至70行是法方程的系数（上三角一维存储）；第73至76行是求逆后的法方程系数；第79至86是坐标改正数和坐标平差值；第89至92行是误差椭圆的参数；第95行是单位权中误差。图4-9a计算结果-lai如IjjslsM|豪詆耳53me|jH1出画1程）51netp田m曰日百32百3斷述IGZ3331?wfit.&ra3447TW.EE35S7T4SL32336T5開迟3B33?35翻1&5&2E

42、wsF.oraS$34,5&r72S,33FtoUii鼻的腔貝止标125370.1360JDDXB2(3a丁齟乳D02S6225.S2.WO15；淞如.理Q閃亞.DM16isseo.2?0和5汕.5師4：.390SSO30.2Mta43762225T?6S.592戲E師.BM3D51讶農方程弔加和糟功52-&.-D.11L-a.DDDD.SS4I.DOQ乩ILL2-HH-CL&Mdcen53-r.DDDD.QQE-C.ODDD.D9T1.000-O.fliK2-QXi-CLOS?CLCOD54-r.ODDD.NE-fi.ODDD.Tfi4n.oao-O.SJC帚cm-

43、clr&i-CLCOD551.-D.ma2.DDD1.DDdn.oao-DLS帚cm-L-iXO-CLCOD5G3.DDDD.B574.DDD-D.ELEE.aao-D.3E?fi.SBX53EQ.CGD573.ODD.99Da.DDDD.42-.oaa-D.9CO-2.C833-OlL程-H.知50-3.ODD-D.GB4-2.non-n.6-0000=6&16-0000=MS3LS3F59-l.DOO-D.0130.000L刖85-0110flB9L35-0000=40890,000-D.39flb.noaD.92L乩QilD乩临6-000-CLSeinmoSJ5.000f除6.DOO-n

44、.Jflfl1.000a.ass2rijjOL-lif-5O.mn02-LD00-d.rwn.noa-n.7(iaTb0iOa.TL-iaoooa-LS.ldS3LDOO0.OFL2.0Q0-G.XTMH-fcrfH：gx，allipc*稠jbi.23T2.521.2.5902.20?Ilfl.QLL90106人6323.025i.1853.2022.刘J2F.33T9L1QT人觸*舊岡I:：.袱:.0TT2.2LQ62104.57&3.355l.呦3.53d2.SSL11-4.QPi93IrSftOlL1J1114Jdith=lEearn94thjivtA2statifiisityxflir

45、ru=3.93Z旳LM匸3.523796-电1Jn|U1?OS3图4-9c计算结4.输出的控制网图和误差椭圆图如下:卜Figure1FileEditViewInsertToolsDesktopWindowHelp勺q包背可口eT测边网平差程序代码：globaledddsdddlpnx0y0mlm2m3msppedsidmdgfaabbccrtrrttqlwglobalpathnamenet_names_datafileb_datafilefilename;xO=;yO=;e=;d=;sid=;g=;f=;pn=;%数据的读入if(isempty(pathname)lisempty(net_na

46、me)filename,pathname=uigetfile(*.txt,Tnputfilename);i=find(.=filename);net_name=filename(1:i-1);endfit1=fopen(strcat(pathname,filename),rt);if(fit1=-1)msgbox(InputFileorPathisnotcorrect,Waming,warn);%openafilereturn;endtoreadb_datafile=out.txtfit2=fopen(strcat(pathname,net_name,b_datafile),wt);if(fi

47、t1=-1)msgbox(InputFileorPathisnotcorrect,Warning,warn);return;enda=fscanf(fit1,%d,2);ed=a(1);dd=a(2);sd=ed+dd;fprintf(fit2,%5dn,ed,dd,dd1);fprintf(fit2,n);pn=fscanf(fit1,%d,sd);pn1=pn;pn1,i2=chkdat(sd,pn,pn1);fprintf(fit2,%5d%5d%5d%5d%5d%5d%5d%5d%5d%5d%5dn,pn);fprintf(fit2,n);a=fscanf(fit1,%f,2*ed);

48、fori=1:edx0(i)=a(2*i-1);y0(i)=a(2*i);fprintf(fit2,%15.3f%15.3fn,x0(i),y0(i);endfprintf(fit2,n);m1=fscanf(fit1,%d,1);fprintf(fit2,%5dn,m1);isid=0;if(m10)fprintf(fit2,siden);a=fscanf(fit1,%f,2);ms=a(1);pp=a(2);fprintf(fit2,%6.2f%6.2fn,ms,pp);a=fscanf(fit1,%d%d%f,3*m1);fori=1:m1e(i)=a(3*i-2);d(i)=a(3*i

49、-1);sid(i)=a(3*i);fprintf(fit2,%5d%5d%15.3fn,e(i),d(i),sid(i);ende,i1=chkdat(sd,pn,e);d,i2=chkdat(sd,pn,d);i3=0;fori=1:m1;if(e(i)=d(i)i3=1;fprintf(fit2,*%5d%5d%5d*n,i,e(i),d(i);endfprintf(fit2,%5d%5d%15.3fn,e(i),d(i),sid(i);endisid=i1+i2+i3;endidir=0;kk=isidif(kk0)msgbox(Errorbyfunctionrddat1,Warnin

50、g,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;fori=1:length(m)ifi=1ife(m(i)=A|e(m(i)=BP=d(m(i);elseP=e(m(i);endelseife(m(i-1)=A|d(m(i-1)=AB=P;elseA=P;endife(m(i)=A|e(m(i)=BP=d(m(i);else

51、P=e(m(i);endendforj=1:m1if(e(j)=P|d(j)=P)&(d(j)=A|e(j)=A)s1=sid(j);endif(e(j)=P|d(j)=P)&(d(j)=B|e(j)=B)s2=sid(j);endenddeltx=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(al

52、faAB);y0(P)=y0(A)+ee*sin(alfaAB)+ff*cos(alfaAB);endfprintf(fit2,n);fprintf(fit2,计算的近似坐标n);fori=1:sdfprintf(fit2,%15.3f%15.3fn,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=;a(1:m1,1:9)=0.0;c(1:sum)=0.0;w(1:n)=0.0;fprintf(fit2,误差方程系数和常数项n);fori=1:m1%边长观测误差方程dx=x0(

53、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)=(ms2+(ss*pp*0.0001)人2);ql(i)=(ms+ss*pp*0.001);ql(i)=100/q

54、l(i)A2;fprintf(fit2,%15.3f%15.3f%15.3f%15.3f%15.3f%15.3f%15.3f%15.3f%15.3fn,a(i,:);end%fprintf(fit2,%15.3fn,ql);fori=1:m1%形成法方程forj=1:4jj=fix(a(i,2*j-1);if(jj=0)continue;endw(jj)=w(jj)+a(i,2*j)*a(i,9)*ql(i);di=(jj-1)*(n-jj/2.0);fork=1:4kk=fix(a(i,2*k-1);if(kkkk)continue;endc(di+kk)=c(di+kk)+a(i,2*k)

55、*a(i,2*j)*ql(i);endendendfprintf(fit2,n);fprintf(fit2,%15.3fn,w);fprintf(fit2，法方程系数n);fprintf(fit2,%15.3f%15.3f%15.3f%15.3f%15.3f%15.3f%15.3f%15.3f%15.3fn，c)%坐标改正数和单位权中误差的计算sd=ed+dd;n=2*dd;k=0;fori=1:ndxy(i)=0.0;di=(i-1)*(n-i/2.0);forj=1:ndj=(j-1)*(n-j/2.0);if(j1.0)%k=1;%enddxy(i)=-dxy(i)/100.0;end%

56、end%while(k11)while(k11)x(1:ed)=x0(1:ed);y(1:ed)=y0(1:ed);fori=1:ddx(ed+i)=x0(ed+i)+dxy(2*i-1);y(ed+i)=y0(ed+i)+dxy(2*i);endx0(1:sd)=x(1:sd);y0(1:sd)=y(1:sd);fori=1:sdif(i=ed)vx(i)=0.0;vy(i)=0.0;elsevx(i)=dxy(2*(i-ed)-1);vy(i)=dxy(2*(i-ed);endendvy%l4.4fnfprintf(fit2,坐标改正数和坐标平差值n);fprintf(fit2,pnvxxyn);fori=l:sdfprintf(fit2,%6d%8.4f%14.4f%8.4f,pn(i),vx(i),x(i),vy(i),y(i);endpvv=0.0;fori=l:np