第八讲测量程序设计在测绘中应用初步性

1、测量程序设计应用测量程序设计应用1.1.测量数据处理：进行平差计算（近似、测量数据处理：进行平差计算（近似、精确）；广义测量平差计算（滤波、配精确）；广义测量平差计算（滤波、配置、抗差估计等）；其它测量中有关计置、抗差估计等）；其它测量中有关计算算2.2.从事科学研究所必要的计算（读研）从事科学研究所必要的计算（读研）等等测量程序设计流程测量程序设计流程一、一、Matlab用于经典平差计算用于经典平差计算 1.条件平差条件平差函数模型函数模型11010rrnrnALA0WAV随机模型随机模型12020PQD minPVVT平差准则平差准则 条件平差就是在满足条件平差就是在满足r个条件方程式的条

2、件下，个条件方程式的条件下，求使函求使函 数数 最小的最小的 值。值。PVVTV1）、列平差值条件方程：）、列平差值条件方程：11010rrnrnALA000022110221102211rLrLrLrbLbLbLbaLaLaLannnnnn), 2 , 1(,nirbaiii条件方程系数条件方程系数000,rba常数项常数项2）、条件方程：）、条件方程：11010rrnrnALA 代入平差值条件方程中，得到代入平差值条件方程中，得到将将VLL11110rrnrWVA为条件方程闭合差为条件方程闭合差1rW0AALW3）、改正数方程：）、改正数方程：按求函数条件极值的方法引入联系系数向量按求函数

3、条件极值的方法引入联系系数向量TrbarkkkK),(1组成拉氏函数，并计算其最小值：组成拉氏函数，并计算其最小值：)(2WAVKPVVTT 将将对对V求一阶导数并令其为零求一阶导数并令其为零 022AKPVVTT 则：则：KAPVTKQAKAPVTT14）、法方程：）、法方程： 由条件方程由条件方程 AV+W=0和改正数方程和改正数方程 ，则，则得：得： 0WKAQAT记作：记作：1110rrrrraaWKNNaa为满秩方阵，为满秩方阵， rARAQARNRTaa)()()(由于由于)()(011AALAQAWNKTaaKQAVTVLLKQAVT5）、检核和精度计算）、检核和精度计算检查检查

4、AV+W=0 ?或或?00ALAAQNQAQQrPVVaaTLT120,例例1 1 教材误差理论与测量平差基础教材误差理论与测量平差基础P74P74 设对下图中的三个内角作同精度观测，得观测值：设对下图中的三个内角作同精度观测，得观测值：L L1 1=42=42o o12122020，L L2 2=78=78o o09090909，L L3 3=59=59o o38384040，试按条件平差求三个内角得平差值。试按条件平差求三个内角得平差值。clcclcdisp(disp(条件平差示例条件平差示例1)1)disp(disp(三角形内角观测值三角形内角观测值) ) L1 = 42 12 20L1

5、 = 42 12 20; ; % %注意格式注意格式L2 = 78 9 9L2 = 78 9 9; ;L3 = 59 38 40L3 = 59 38 40; ;L = L1; L2; L3L = L1; L2; L3; ;disp(disp(度分秒转换为弧度度分秒转换为弧度)LL = dms2rad(mat2dms(L)LL = dms2rad(mat2dms(L)A = 1 1 1A = 1 1 1; ;w = sum(LL(:) w = sum(LL(:) pi pi ; ; % %闭合差闭合差w = dms2mat(rad2dms(ww = dms2mat(rad2dms(w); ; %

6、 %弧度化为度分秒弧度化为度分秒P = eye(3);P = eye(3);Naa = ANaa = A* *inv(P)inv(P)* *AA; ;K = -inv(Naa)K = -inv(Naa)* *w w; ; V =inv(PV =inv(P) )* * A A* *K K; ; LP = L + VLP = L + V LP = LP = dms2rad(mat2dms(LP)dms2rad(mat2dms(LP) sumLP = sum(LP) sumLP = sum(LP) if(sum(LP) = pi) if(sum(LP) = pi) disp( disp(检核正确检核

7、正确) else else disp( disp(检核错误检核错误) end end例例2 2 教材误差理论与测量平差基础教材误差理论与测量平差基础P75P75在下图中，、为已知水准点，其高程为在下图中，、为已知水准点，其高程为HA=12.013m, HA=12.013m, HB = 10.013mHB = 10.013m。为了确定点及点的高程，共观测了。为了确定点及点的高程，共观测了四个高差，高差观测值及相应的水准路线的距离为：四个高差，高差观测值及相应的水准路线的距离为：h1 = -1.004m, S1 = 2km;h1 = -1.004m, S1 = 2km;h2 = 1.516m, S

8、2 = 1km;h2 = 1.516m, S2 = 1km;h3 = 2.512m, S3 = 2km;h3 = 2.512m, S3 = 2km;h4 = 1.520m, S4 = 1.5kmh4 = 1.520m, S4 = 1.5km试求和点高程的平差值。试求和点高程的平差值。clcclcclearclearh1 = -1.004;h1 = -1.004;h2 = 1.516;h2 = 1.516;h3 = 2.512;h3 = 2.512;h4 = 1.520;h4 = 1.520;HA = 12.013HA = 12.013HB = 10.013HB = 10.013L = h1 h

9、2 h3 h4L = h1 h2 h3 h4s1 = 2; s1 = 2; s2 = 1;s2 = 1;s3 = 2;s3 = 2;s4 = 1.5;s4 = 1.5;s=s1 s2 s3 s4s=s1 s2 s3 s4A = 1 1 -1 0; 0 1 0 -1A = 1 1 -1 0; 0 1 0 -1w1 = h1 + h2 - h3 + HA - HB;w1 = h1 + h2 - h3 + HA - HB;w2 = h2 - h4;w2 = h2 - h4;w = w1; w2w = w1; w2; ;P = diag(1./s)P = diag(1./s); ;Naa = ANaa

10、 = A* *inv(P)inv(P)* *AA; ;K = -inv(Naa)K = -inv(Naa)* *w w; ;V = inv(P)V = inv(P)* *AA* *K K; ;LP= L + V;LP= L + V;if if LP(1,1)+LP(2,1)-LP(3,1)+HA-HB=0 & LP(1,1)+LP(2,1)-LP(3,1)+HA-HB=0 & LP(2,1)-LP(4,1)=0LP(2,1)-LP(4,1)=0 disp( disp(检核正确检核正确)elseelse disp( disp(检核错误检核错误)endend disp( disp(平差后的高程值平

11、差后的高程值) HC = HA + LP(1,1) HC = HA + LP(1,1) HD = HA + LP(1,1) + LP(4,1) HD = HA + LP(1,1) + LP(4,1)例例3 如图如图3.1所示水准网，所示水准网，A、B两点为高程已知，各两点为高程已知，各观测高差及路线长度如表观测高差及路线长度如表3.1所列。用所列。用条件平差法条件平差法计计算求知点的高程平差值及算求知点的高程平差值及p2和和p3之间平差后高差值之间平差后高差值的中误差。的中误差。表表3.1解解 n=7,t=3,r=n-t=4.平差值条件方程为：平差值条件方程为：000031643765521B

12、AHhhHhhhhhhhhh改正数条件方程为：改正数条件方程为：0403070731643765521vvvvvvvvvvv令令C=1，观测值的权倒数为：，观测值的权倒数为：21122111PKAPVWNKAAPNTT111,-)(mmVT2.0 1.2- 3.8- 0.3- 4.4- 2.8 0.43-得：得：7h76543216570110000hhhhhhhhhh0.7416AQfNQAfQffQfTTTffT1,01100005.1mm7ffThQrPVV下面求平差后下面求平差后的中误差：的中误差：中误差为中误差为利用利用Matlab编写上题的代码为：编写上题的代码为：clear A=

13、-1 1 0 0 -1 0 0;0 0 0 0 1 1 -1;0 0 1 -1 0 -1 0;1 0 -1 0 0 0 0;q=1 1 2 2 1 1 2;Q=diag(q); %生成对角矩阵生成对角矩阵QP=inv(Q); %inv()为矩阵求逆运算为矩阵求逆运算W=-7;7;3;-4;N=A*Q*A; K=-inv(N)*W;V=Q*A*KL=1.359; 2.009; 0.363; -0.640; 0.657; 1.000; 1.650;LP=L+(V/1000)f=0;0;0;0;1;1;0;sigma=sqrt(V*P*V)/4) Qff=f*Q*f-f*Q*A*inv(N)*A*Q

14、*f;Dff=sigma*sqrt(Qff)vfunction Lv,V,d0=szwtjpc(L,A,W,P) vnl=length(L); vma, na=size(A);vnw=length(W);vif nl=nav disp( 观测值个数与改正数个数不相等观测值个数与改正数个数不相等! )v returnvendvif ma=nwv warning(所列条件方程有误所列条件方程有误) v returnvendv mp,np=size(P);v if mp=npv warning(权阵应为方阵权阵应为方阵) v returnvendvQ=inv(P);vV=-Q*A*inv(A*Q*A

15、)*WvLv=L+V/1000vd0=sqrt(V*P*V/nw)v%运行时必须把函数放在运行时必须把函数放在G盘中，更改路径盘中，更改路径2、间接平差、间接平差 在一个控制网中，设有在一个控制网中，设有t t个独立参数，将每一个个独立参数，将每一个观测值都表达成所选参数的函数，以此为基础进行平观测值都表达成所选参数的函数，以此为基础进行平差，求得参数的估计值。差，求得参数的估计值。选择参数应做到选择参数应做到足数足数（参数的个数等于必要观测（参数的个数等于必要观测数）和数）和独立独立（参数间不存在函数关系）。利用参（参数间不存在函数关系）。利用参数将观测值表示为数将观测值表示为dXAL其中其

16、中L L为观测值，为观测值，为误差，或者表示为为误差，或者表示为lXA其中其中l lL Ld.d.设误差设误差和参数和参数X X的估计值分别为的估计值分别为V V 和和X则有则有lXAVxXX0若则误差方程表示为则误差方程表示为lxAV其中常数项为其中常数项为)(0dAXLl由最小二乘准则，所求参数的改正数应该满足由最小二乘准则，所求参数的改正数应该满足minPVVT目标函数对目标函数对x x求一阶导数，并令其为零求一阶导数，并令其为零022)(PAVxddVPVxdPVVdTTT转置后得到转置后得到0PVAT把误差方程代入上式后得到把误差方程代入上式后得到0PlAxPAATT设设PAANTA

17、PlAWT则法方程为则法方程为0WxNA由此求得参数改正数的唯一解为由此求得参数改正数的唯一解为WNxA1将其代入误差方程，可求得改正数将其代入误差方程，可求得改正数V V，最后得到观测，最后得到观测值得平差值及参数的平差值分别为值得平差值及参数的平差值分别为VLLxXX0精度评定精度评定单位权中误差的估计公式为单位权中误差的估计公式为tnPVVT0由协因数传播定律得参数的权逆阵为由协因数传播定律得参数的权逆阵为1 AxxXXNQQ参数的中误差为参数的中误差为jjjXXXQ0若设参数的函数为若设参数的函数为),(21tXXXfF以以xXX0代入，并按泰勒级数展开，取一次项代入，并按泰勒级数展开

18、，取一次项ttxfxffF110记为记为xfxfxfTttF11由此得到参数函数的权逆阵为由此得到参数函数的权逆阵为fNffQfQATXXTFF1中误差为中误差为FFFQ0改正数为改正数为lIPAANlxAVTA)(1TAVVAANQQ1平差值及其协因数阵分别为平差值及其协因数阵分别为VLLTAVVVLLVLLLLAANQQQQQ1例例1 1 教材教材P124P124 A A和和B B是已知高程的水是已知高程的水准点，准点，C C、D D、E E是待定点，是待定点，A A、B B高程、观测高差和高程、观测高差和相应的水准路线长度见相应的水准路线长度见下表，试按间接平差求下表，试按间接平差求各待

19、定点的高程平差值。各待定点的高程平差值。程序代码如下：程序代码如下：disp(水准网间接平差示例水准网间接平差示例)disp(已知高程已知高程)Ha = 5.015 % 已知点高程已知点高程Hb = 6.016 % 已知点高程已知点高程disp(观测高差，单位观测高差，单位m)L = 1.359;2.009; 0.363; 1.012; 0.657; -0.357 disp(系数矩阵系数矩阵B)B = 1 0; 0 1; 1 0; 0 1; -1 1; -1 0l = 0; 0; 4; 3; 7; 2;S = 1.1, 1.7, 2.3, 2.7, 2.4, 4.0;P = C./S ; P

20、= diag(P); % 定义权阵定义权阵disp(参数的解参数的解)x = inv(B*P*B)*B*P*ldisp(误差误差V(mm), 各待定点的高程平差值各待定点的高程平差值L1（m）)V = B*x - l ; % 误差方程误差方程(mm)L1 = L + V/1000 % 观测值的平差值，观测值的平差值，disp(精度评定精度评定)n = 6; % 观测值的个数观测值的个数t = 2; % 必要观测数必要观测数delta = sqrt(V*P*V/(n t)disp(C是单位权观测高差的线路公里数，是单位权观测高差的线路公里数，S是线路长度是线路长度)C = ones(1,6);d

21、isp(水准网间接平差示例水准网间接平差示例)disp(已知高程已知高程)Ha = 5.015 % 已知点高程已知点高程Hb = 6.016 % 已知点高程已知点高程disp(观测高差观测高差)L = 1.359;2.009; 0.363; 1.012; 0.657; -0.357 function jiesuan2( )fid=fopen(E:szjj.txt,r) %注意路径注意路径temp=fscanf(fid,%f%f%f%f,4,6)fclose(fid); temp=temp; B=temp(:,1,2); l=temp(:,3); S=temp(:,4); S=S; disp(C

22、是单位线路公里数，是单位线路公里数，S是线路长度是线路长度) C = ones(1,6)P = C./S ; P = diag(P) % 定义权阵定义权阵disp(参数的解参数的解)x=inv(B*P*B)*B*P*l disp(误差误差V(mm), 各待定点的高程平差值各待定点的高程平差值L1（m）) V=B*x-l % 误差方程误差方程(mm)L1 = L + V/1000 % 观测值的平差值，观测值的平差值，disp(精度评定精度评定)n = 6; % 观测值的个数观测值的个数t = 2; % 必要观测数必要观测数delta = sqrt(V*P*V/(n t)例例2. 图图4.3中中,

23、A,B,C是已知点是已知点, 为待定点为待定点,网中网中观测了观测了12个角度和个角度和6条边长。已知测角中误差为条边长。已知测角中误差为 ，边长测量中误差为边长测量中误差为 cm，起算数据及观测值分，起算数据及观测值分别列表于表别列表于表4.1和表和表4.2。21PP、5 . 1 0 . 2ABC 表表4.1表表4.2.2211TYXYXX )./(5625. 025 . 1222222cmpSS秒.11814418118lxBV.)(1PlBPBBxTT解：本题解：本题n =18,即有即有18个误差方程，其中有个误差方程，其中有12个角度误差个角度误差方程和方程和4个边长误差方程。必要观测

24、数个边长误差方程。必要观测数t=4.采用间接平差采用间接平差，现取待定点坐标的平差值为参数，即现取待定点坐标的平差值为参数，即以测角中误差为该边角网的单位权观测值中误差，则测边的权为：以测角中误差为该边角网的单位权观测值中误差，则测边的权为：解题的关键是列误差方程解题的关键是列误差方程在此基础上解出在此基础上解出完成精度评定。完成精度评定。解题具体过程如下。解题具体过程如下。21APAP和.,73312LLLABAPABAP.sin,cos11110101APAPAAPAPASYYSXX（一）。首先计算待定点的（一）。首先计算待定点的近似近似坐标坐标。利用观测值和已知坐标求出。利用观测值和已知

25、坐标求出的近似坐标方位角，具体计算公式为的近似坐标方位角，具体计算公式为进一步可由坐标增量公式进一步可由坐标增量公式求得待定的近似坐标为：求得待定的近似坐标为：)556.2475,928.5656(),099.2944,875.663(0201PP点近似坐标类似计算2PkjkjkkjkjkjjkjkjjkjkjkySXxSYySXxSY)()()()(200200200200（二）（二）.计算坐标方位角改正数方程的系数。计算坐标方位角改正数方程的系数。按按计算。计算计算。计算时时以米为单位，以米为单位，而而因数值较小，因数值较小，一般一般采用厘米为单位，此时采用厘米为单位，此时有关计算见表有关

26、计算见表-1-1。1 x1 y2 x2 y表表-1000jkjkjkYXS、kkjjyxyx,65.2062 （三）（三）.计算测边的误差方程的系数。即计算测边的误差方程的系数。即ikjkjkkjkjkjjkjkjjkjkilySYxSXySYxSXv00000000的系数。的系数。1 x1 y2 x2 y表表-2)(00jhjkiiLl0jkiiSLll（四）（四）.确定误差方程。在计算确定误差方程。在计算或或时，利用近似坐标求出每个边的近似坐标方位角和近似边长。时，利用近似坐标求出每个边的近似坐标方位角和近似边长。从而可确定从而可确定的系数。进一步可得出误差方程的系数项，见表的系数。进一步

27、可得出误差方程的系数项，见表-3。1 x1 y2 x2 ylpiiS表表-3可得法方程为可得法方程为（五）（五）.平差值计算与精度评定。由法方程可解得参数改正数为：平差值计算与精度评定。由法方程可解得参数改正数为：2965. 00249. 00481. 00074. 02709. 00775. 01583. 09888. 00137. 07970. 0,8997. 03647. 38065. 65297. 312211对称bbNyxyx02148. 19415. 81620. 67469. 24194. 32866. 01437. 00229. 02866. 03174. 43363. 086

28、08. 01437. 03363. 00459. 10861. 00229. 08608. 00861. 04270. 12211yxyxx lxBV21,PP由上式可解出由上式可解出，再根据，再根据得各个改正数为得各个改正数为V= 1.9 -1.6 1.4 -3.2 0.7 -4.8 1.3 4.3 -2.5 0 6.1 -2.8 -10.7 -8.2 -3.4 -1.8 2.3T进而可求出观测值的平差值。进而可求出观测值的平差值。单位权中误差单位权中误差的枞横坐标的协因数和的枞横坐标的协因数和点位中误差点位中误差分别为分别为2967. 30tnPVVT2965. 0,2709. 0,988

29、8. 0,7970. 022221111YYXXYYXXQQQQ2.4833(cm)4.1622(cm)222202111101YYXXPYYXXPQQQQ利用利用Matlab编写上题的代码为：编写上题的代码为：clearB=0 0 0.5532 -0.8055;0 0 0.2434 0.5528;0 0 -0.7966 0.2572;0 0 -0.2434. -0.5528;0 0 0.6298 0.6368;0 0 -0.3864 -0.084;-0.2244 -0.3379 0.7966. -0.2572;0.0384 0.4095 -0.835 -0.1523;0.186 0.0716

30、 0.0384 0.4095;.0.2998 0.1901 -0.0384 -0.4095;-0.0384 -0.4095 -0.348 0.3255;.-0.2614 0.2194 0.3864 0.084;0 0 0.307 0.952;0 0 -0.915 0.403;0 0 0.212. -0.977;-0.643 -0.766 0 0;-0.833 0.553 0 0;-0.996 0.093 0.996 -0.093;l=-2.6;-0.11;-0.02;0.48;0.30;0.42;0.02;0.05;-0.97;0.24;1.05;-5.09;0.1;12.1.;9.7;10.

31、7;0;0.1;X0=663.875,2944.099,5656.928,2475.556;p=1 1 1 1 1 1 1 1 1 1 1 1 0.5625 0.5625 0.5625 0.5625 0.5625 0.5625;P=diag(p);N=B*P*B;W=B*P*l;x=inv(N)*W;XP=X0+x/100invN=inv(N) v=B*x-l;sigma=sqrt(l*P*l-W*x)/14)sigmaP1=sigma*sqrt(invN(1,1)+invN(2,2)sigmaP2=sigma*sqrt(invN(3,3)+invN(4,4)v二、二、MATLABMATLAB

32、在数字测量技术中简单应用在数字测量技术中简单应用v1.1.坐标正、反算坐标正、反算v1 1）坐标正算）坐标正算v 坐标正算即已知坐标正算即已知A A点坐标点坐标 ,A ,A点至点至B B点水平点水平距离距离 以及坐标方位角以及坐标方位角 ，求，求B B点坐标点坐标 的过程。计算公式如下：的过程。计算公式如下：v ABSABABAABABABABAABABAByySyyxxSxxsincosvax=input(请输入请输入A的纵坐标的纵坐标:)vay=input(请输入请输入A点的横坐标点的横坐标:)vs=input(请输入请输入AB长度长度:)vad=input(请输入请输入AB的坐标方位角度

33、的坐标方位角度:)vam=input(请输入请输入AB的坐标方位角分的坐标方位角分:)vas=input(请输入请输入AB的坐标方位角秒的坐标方位角秒:)va=ad+am/60+as/3600;va=a*pi/180;vbx = ax + s * cos(a) vby = ay + s * sin(a) v2）坐标反算坐标反算v设设A A、B B两已知点的坐标分别为两已知点的坐标分别为 和和 v，则直线，则直线ABAB的坐标方位角的坐标方位角 和水平距离和水平距离 为：为：v v ),(AAyx),(BByx2222)()(arctanarctanABABABABABABABABABAByxy

34、yxxSxyxxyyvclearvclcvx=2507.69 2299.82 2192.43 2179.72 2166.72;vy=1215.63 1303.80 1383.98 1556.41 1757.29;vplot(x,y)vhold onvfor i=1:5v plot(x(i),y(i),o)vendvfor i=1:4 a(i)=x(i+1)-x(i); b(i)=y(i+1)-y(i); s(i)=sqrt(a(i)2+b(i)2);endvs %给出各边长给出各边长vfor i=1:4v if a(i)=0 m(i)=atand(abs(b(i)/a(i); v %结果是度，

35、得到是与坐标纵轴结果是度，得到是与坐标纵轴x轴构成的锐角轴构成的锐角 k(i)=m(i); if a(i)0 & b(i)0 c(i)=k(i); elseif a(i)0 c(i)=180-k(i); elseif a(i)0 & b(i)0v c(i)=90;v elsev c(i)=270;v endv endv d=fix(c)v e=fix(c-d)*60)v f=(c-d)*60-e)*60v g=d e fv v 2.可以进行地形图绘制（若有可以进行地形图绘制（若有11个碎部点坐标）个碎部点坐标）v 115.14v129.56v138.66v2112.55v2210.51v232

36、.70v3122.03v3222.97v3332.94v4255.58v4324.55vfid=fopen(E:dengzhixian.txt,r);vdata = fscanf(fid,%f %f %f,3,inf);v%按行取，按列存按行取，按列存vdata=data;vfclose(fid);vx =data(:,1);y =data(:,2);z =data(:,3);vxb=min(x); xe=max(x); %得到最小值最大值得到最小值最大值vyb=min(y); ye=max(y); %得到最小值最大值得到最小值最大值vXl=linspace(xb,xe,200); vYl=l

37、inspace(yb,ye,200); vX,Y=meshgrid(Xl,Yl); %坐标矩阵坐标矩阵vmethod=cubic; %插值方法插值方法vZ =griddata(x,y,z,X,Y,method);vfigure; %准备图形窗口准备图形窗口vhold on; vcontour(X,Y,Z);%画等高线画等高线vcolormap cool; vbox on; %加边框加边框vxlabel(x);ylabel(y);vaxis equal;vtitle(等值线图等值线图:(单位单位:scale);vhold offv3.导线网近似平差计算导线网近似平差计算v 规范规定，对于二级及以

38、下等级的图根导线规范规定，对于二级及以下等级的图根导线允许对单一导线、单结点导线网采用近似平差方允许对单一导线、单结点导线网采用近似平差方法进行计算。导线近似平差的基本思路是将角度法进行计算。导线近似平差的基本思路是将角度观测值和边长观测值观测值和边长观测值分别分别进行平差计算。进行平差计算。v附合导线附合导线坐标近似平差计算步骤如下：坐标近似平差计算步骤如下：v（1 1）方位角计算；）方位角计算；v（2 2）由转角公式计算其角度闭合差；）由转角公式计算其角度闭合差；v（3 3）将角度闭合差除以测边数分配到各观测角中；将角度闭合差除以测边数分配到各观测角中；v（4 4）计算调整后的导线边的方位

39、角；）计算调整后的导线边的方位角；v（5）计算坐标增量；）计算坐标增量；v（6）计算坐标增量闭合差并分配；）计算坐标增量闭合差并分配；v（7）用调整后的方位角计算各点坐标。）用调整后的方位角计算各点坐标。v对于闭合导线和支导线的情况可类似计算对于闭合导线和支导线的情况可类似计算v注意：在进行控制网精平时，一般采用间接注意：在进行控制网精平时，一般采用间接平差，误差方程形式统一、规律性强、便于平差，误差方程形式统一、规律性强、便于程序开发；而使用条件平差等，其方程形式程序开发；而使用条件平差等，其方程形式多样、规律性差，不利于程序开发。所以，多样、规律性差，不利于程序开发。所以，一般建议使用间接

40、平差方法进行程序编写。一般建议使用间接平差方法进行程序编写。 如下图所示，已知四控制点平面坐标分如下图所示，已知四控制点平面坐标分别为：别为：A（599.951,224.856）B（704.816,141.165）、）、C（747.166,572.726）、）、D(889.339,622.134)，各观测角数据如，各观测角数据如下：下：1 =741030，2 =2790512，3=675529，4 =2761011，5 =802346，其导，其导线各边长分别为：线各边长分别为：S1=143.825m ，S2=124.777m,S3=188.950m,S4=117.338m,vA Dv 2 4v s1v 1 s2 s3 s4 5v B Cv 3vn=input(请输入导线边数请输入导线边数); %观测数据输入观测数据输入vfor i=1:nv a=num2str(i);v string1=(char(请输入第请输入第),a,(条导线边长条导线边长); v bian(i)=input(string1);vendvfor k=1:(1+n)v s=num2str(k);v string2=(char(请输入第请输入第),s,(个转折角个转折角); v b=input(string2,s); %度分秒为字符串度分秒为字符串v du=str2num(b(