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

1、第八讲测量程序设计在测绘中应用 初步性 测量程序设计应用测量程序设计应用 1.1.测量数据处理：进行平差计算（近似、测量数据处理：进行平差计算（近似、 精确）；广义测量平差计算（滤波、配精确）；广义测量平差计算（滤波、配 置、抗差估计等）；其它测量中有关计置、抗差估计等）；其它测量中有关计 算算 2.2.从事科学研究所必要的计算（读研）从事科学研究所必要的计算（读研） 等等 第八讲测量程序设计在测绘中应用 初步性 测量程序设计流程测量程序设计流程 第八讲测量程序设计在测绘中应用 初步性 一、一、Matlab用于经典平差计算用于经典平差计算 1.条件平差条件平差 函数模型函数模型 1 1 0 1

2、 0 r r nrn ALA 0WAV 随机模型随机模型 12 0 2 0 PQD minPVV T平差准则平差准则 条件平差就是在满足条件平差就是在满足r个条件方程式的条件下，个条件方程式的条件下， 求使函求使函 数数 最小的最小的 值。值。 PVV T V 第八讲测量程序设计在测绘中应用 初步性 1）、列平差值条件方程：）、列平差值条件方程： 1 1 0 1 0 r r nrn ALA 0 0 0 02211 02211 02211 rLrLrLr bLbLbLb aLaLaLa nn nn nn ), 2 , 1(,nirba iii 条件方程系数条件方程系数 000 ,rba常数项常数

3、项 第八讲测量程序设计在测绘中应用 初步性 2）、条件方程：）、条件方程： 1 1 0 1 0 r r nrn ALA 代入平差值条件方程中，得到代入平差值条件方程中，得到 将将 VLL 1111 0 rrnr WVA 为条件方程闭合差为条件方程闭合差 1r W 0 AALW 第八讲测量程序设计在测绘中应用 初步性 3）、改正数方程：）、改正数方程： 按求函数条件极值的方法引入联系系数向量按求函数条件极值的方法引入联系系数向量 T rba r kkkK),( 1 组成拉氏函数，并计算其最小值：组成拉氏函数，并计算其最小值： )(2WAVKPVV TT 将将对对V求一阶导数并令其为零求一阶导数并

4、令其为零 022 AKPV V TT 则：则： KAPV T KQAKAPV TT 1 第八讲测量程序设计在测绘中应用 初步性 4）、法方程：）、法方程： 由条件方程由条件方程 AV+W=0和改正数方程和改正数方程 ，则得：，则得： 0WKAQAT 记作：记作： 111 0 rrr rr aa WKN Naa为满秩方阵，为满秩方阵， rARAQARNR T aa )()()( 由于由于 )()( 0 11 AALAQAWNK T aa KQAV T VLL KQAV T 第八讲测量程序设计在测绘中应用 初步性 5）、检核和精度计算）、检核和精度计算 检查检查AV+W=0 ?或或?0 0 ALA

5、 AQNQAQQ r PVV aa T L T 12 0 , 第八讲测量程序设计在测绘中应用 初步性 例例1 1 教材误差理论与测量平差基础教材误差理论与测量平差基础P74P74 设对下图中的三个内角作同精度观测，得观测值：设对下图中的三个内角作同精度观测，得观测值： L L1 1=42=42o o12122020，L L2 2=78=78o o09090909，L L3 3=59=59o o38384040， 试按条件平差求三个内角得平差值。试按条件平差求三个内角得平差值。 第八讲测量程序设计在测绘中应用 初步性 clcclc disp(disp(条件平差示例条件平差示例1)1) disp(

6、disp(三角形内角观测值三角形内角观测值) ) L1 = 42 12 20L1 = 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(:)

7、pi pi ; ; % %闭合差闭合差 w = dms2mat(rad2dms(ww = dms2mat(rad2dms(w); ; % %弧度化为度分秒弧度化为度分秒 P = eye(3);P = eye(3); Naa = ANaa = A* *inv(P)inv(P)* *A A; ; 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

8、(LP) sumLP = sum(LP) if(sum(LP) = pi) if(sum(LP) = pi) disp( disp(检核正确检核正确) else else disp( disp(检核错误检核错误) end end 第八讲测量程序设计在测绘中应用 初步性 例例2 2 教材误差理论与测量平差基础教材误差理论与测量平差基础P75P75 在下图中，、为已知水准点，其高程为在下图中，、为已知水准点，其高程为HA=12.013m, HA=12.013m, HB = 10.013mHB = 10.013m。为了确定点及点的高程，共观测了。为了确定点及点的高程，共观测了 四个高差，高差观测值及

9、相应的水准路线的距离为：四个高差，高差观测值及相应的水准路线的距离为： h1 = -1.004m, S1 = 2km;h1 = -1.004m, S1 = 2km; h2 = 1.516m, S2 = 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 试求和点高程的平差值。试求和点高程的平差值。 第八讲测量程序设计在测绘中应用 初步性 clcclc clearclear h1 = -1.004;h1 = -1

10、.004; h2 = 1.516;h2 = 1.516; h3 = 2.512;h3 = 2.512; h4 = 1.520;h4 = 1.520; HA = 12.013HA = 12.013 HB = 10.013HB = 10.013 L = h1 h2 h3 h4L = h1 h2 h3 h4 s1 = 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 s4 A = 1 1 -1 0; 0 1 0 -1A = 1 1 -1 0; 0 1 0 -1 w1 = h1 + h

11、2 - 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 = A* *inv(P)inv(P)* *A A; ; K = -inv(Naa)K = -inv(Naa)* *w w; ; V = inv(P)V = inv(P)* *A A* *K K; ; LP= L + V;LP= L + V; 第八讲测量程序设计在测绘中应用 初步性 if if LP(1,1)+LP(2,

12、1)-LP(3,1)+HA-HB=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); %生成对角矩阵生成对角矩阵Q P=inv(Q); %inv()为矩阵求逆运算为矩阵求逆运算 W=-7;7;3;-4; N=A*Q*A; K=-inv(N)*W; V=Q*A*K L=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

13、*A*inv(N)*A*Q*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=na v disp( 观测值个数与改正数个数不相等观测值个数与改正数个数不相等! ) v return vend vif ma=nw v warning(所列条件方程有误所列条件方程有误) v return vend 第八讲测量程序设计在测绘中应用 初步性 v mp,np=size(P); v if mp=

14、np v warning(权阵应为方阵权阵应为方阵) v return vend vQ=inv(P); vV=-Q*A*inv(A*Q*A)*W vLv=L+V/1000 vd0=sqrt(V*P*V/nw) v%运行时必须把函数放在运行时必须把函数放在G盘中，更改路径盘中，更改路径 第八讲测量程序设计在测绘中应用 初步性 2、间接平差、间接平差 在一个控制网中，设有在一个控制网中，设有t t个独立参数，将每一个个独立参数，将每一个 观测值都表达成所选参数的函数，以此为基础进行平观测值都表达成所选参数的函数，以此为基础进行平 差，求得参数的估计值。差，求得参数的估计值。 选择参数应做到足数（参

15、数的个数等于必要观测选择参数应做到足数（参数的个数等于必要观测 数）和独立（参数间不存在函数关系）。利用参数）和独立（参数间不存在函数关系）。利用参 数将观测值表示为数将观测值表示为 dXAL 其中其中L L为观测值，为观测值，为误差，或者表示为为误差，或者表示为 lXA 其中其中l lL Ld.d. 第八讲测量程序设计在测绘中应用 初步性 设误差设误差和参数和参数X X的估计值分别为的估计值分别为V V 和和 X 则有则有 lXAV xXX 0 若 则误差方程表示为则误差方程表示为 lxAV 其中常数项为其中常数项为 )( 0 dAXLl 第八讲测量程序设计在测绘中应用 初步性 由最小二乘准

16、则，所求参数的改正数应该满足由最小二乘准则，所求参数的改正数应该满足 minPVV T 目标函数对目标函数对x x求一阶导数，并令其为零求一阶导数，并令其为零 02 2 )( PAV xd dV PV xd PVVd TT T 转置后得到转置后得到 0PVAT 把误差方程代入上式后得到把误差方程代入上式后得到 0PlAxPAA TT 第八讲测量程序设计在测绘中应用 初步性 设设 PAAN T A PlAW T 则法方程为则法方程为0WxN A 由此求得参数改正数的唯一解为由此求得参数改正数的唯一解为 WNx A 1 将其代入误差方程，可求得改正数将其代入误差方程，可求得改正数V V，最后得到观

17、测，最后得到观测 值得平差值及参数的平差值分别为值得平差值及参数的平差值分别为 VLL xXX 0 第八讲测量程序设计在测绘中应用 初步性 精度评定精度评定 单位权中误差的估计公式为单位权中误差的估计公式为 tn PVV T 0 由协因数传播定律得参数的权逆阵为由协因数传播定律得参数的权逆阵为 1 Axx XX NQQ 参数的中误差为参数的中误差为 jj jXX X Q 0 第八讲测量程序设计在测绘中应用 初步性 若设参数的函数为若设参数的函数为 ) , , ( 21t XXXfF 以以 xXX 0 代入，并按泰勒级数展开，取一次项代入，并按泰勒级数展开，取一次项 tt xfxffF 110

18、记为记为 xfxfxf T tt F 1 1 由此得到参数函数的权逆阵为由此得到参数函数的权逆阵为 fNffQfQ A T XX T FF 1 第八讲测量程序设计在测绘中应用 初步性 中误差为中误差为 FFF Q 0 改正数为改正数为 lIPAANlxAV T A )( 1 T AVV AANQQ 1 平差值及其协因数阵分别为平差值及其协因数阵分别为 VLL T AVVVLLVLL LL AANQQQQQ 1 第八讲测量程序设计在测绘中应用 初步性 例例1 1 教材教材P124P124 A A和和B B是已知高程的水是已知高程的水 准点，准点，C C、D D、E E是待定点，是待定点， A A

19、、B B高程、观测高差和高程、观测高差和 相应的水准路线长度见相应的水准路线长度见 下表，试按间接平差求下表，试按间接平差求 各待定点的高程平差值。各待定点的高程平差值。 第八讲测量程序设计在测绘中应用 初步性 线路编号观测高差 （m） 线路长度 （km） 已知高程 （m） 1 2 3 4 5 6 7 +1.359 +2.009 +0.363 +1.012 +0.657 +0.238 -0.595 1.1 1.7 2.3 2.7 2.4 1.4 2.6 HA=5.016 HB=6.016 第八讲测量程序设计在测绘中应用 初步性 程序代码如下：程序代码如下： disp(水准网间接平差示例水准网间

20、接平差示例) 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 0 l = 0; 0; 4; 3; 7; 2; 第八讲测量程序设计在测绘中应用 初步性 S = 1.1, 1.7, 2.3, 2.7, 2.4, 4.0; P = C./S ; P = diag(P); % 定义权阵定义权阵 d

21、isp(参数的解参数的解) 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) disp(C是单位权观测高差的线路公里数，是单位权观测高差的线路公里数，S是线路长度是线路长度) C = ones(1,6); 第八讲测量程序设计在测绘

22、中应用 初步性 disp(水准网间接平差示例水准网间接平差示例) 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);

23、l=temp(:,3); S=temp(:,4); S=S; disp(C是单位线路公里数，是单位线路公里数，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; % 观测值的个

24、数观测值的个数 t = 2; % 必要观测数必要观测数 delta = sqrt(V*P*V/(n t) 第八讲测量程序设计在测绘中应用 初步性 例例2. 图图4.3中中,A,B,C是已知点是已知点, 为待定点为待定点,网中网中 观测了观测了12个角度和个角度和6条边长。已知测角中误差为条边长。已知测角中误差为 ， 边长测量中误差为边长测量中误差为 cm，起算数据及观测值分，起算数据及观测值分 别列表于表别列表于表4.1和表和表4.2。 21 PP、 5 . 1 0 . 2 第八讲测量程序设计在测绘中应用 初步性 A B C 点号点号 坐坐 标标 (m) 坐标方位角坐标方位角 边边 长长 (m

25、) xy 4899.846 8781.945 4548.795 130.812 1099.443 7572.622 14 00 35.77 123 10 57.97 4001.117 7734.443 表表4.1 第八讲测量程序设计在测绘中应用 初步性 编号编号观测角观测角编号编号观测角观测角编号编号观测边观测边 （m） 1 2 3 4 5 6 84 07 38.2 37 46 34.9 58 05 44.1 33 03 03.2 126 01 55.7 20 55 02.3 7 8 9 10 11 12 74 18 16.8 77 27 59.1 28 13 43.2 55 21 09.9

26、72 22 25.8 52 16 20.5 13 14 15 16 17 18 2463.94 3414.71 5216.23 6042.94 5085.08 5014.99 表表4.2 第八讲测量程序设计在测绘中应用 初步性 . 2211 T YXYXX )./(5625. 0 2 5 . 1 22 2 2 2 2 cmp S S 秒 . 11814418118 lxBV .)( 1 PlBPBBx TT 解：本题解：本题n =18,即有即有18个误差方程，其中有个误差方程，其中有12个角度误差个角度误差 方程和方程和4个边长误差方程。必要观测数个边长误差方程。必要观测数t=4.采用间接平差

27、，采用间接平差， 现取待定点坐标的平差值为参数，即现取待定点坐标的平差值为参数，即 以测角中误差为该边角网的单位权观测值中误差，则测边的权为：以测角中误差为该边角网的单位权观测值中误差，则测边的权为： 解题的关键是列误差方程解题的关键是列误差方程 在此基础上解出在此基础上解出 完成精度评定。完成精度评定。 解题具体过程如下。解题具体过程如下。 第八讲测量程序设计在测绘中应用 初步性 21 APAP和 ., 733 12 LLL ABAPABAP .sin,cos 1111 0 1 0 1APAPAAPAPA SYYSXX （一）。首先计算待定点的近似坐标。利用观测值和已知坐标求出（一）。首先计

28、算待定点的近似坐标。利用观测值和已知坐标求出 的近似坐标方位角，具体计算公式为的近似坐标方位角，具体计算公式为 进一步可由坐标增量公式进一步可由坐标增量公式 求得待定的近似坐标为：求得待定的近似坐标为： )556.2475,928.5656(),099.2944,875.663( 0 2 0 1 PP 点近似坐标类似计算 2 P 第八讲测量程序设计在测绘中应用 初步性 k jk jk k jk jk j jk jk j jk jk jk y S X x S Y y S X x S Y )( )( )( )( 20 0 20 0 20 0 20 0 （二）（二）.计算坐标方位角改正数方程的系数。

29、按计算坐标方位角改正数方程的系数。按 计算。计算时计算。计算时 以米为单位，而以米为单位，而 因数值较小，因数值较小，一般一般采用厘米为单位，此时采用厘米为单位，此时 有关计算见表有关计算见表-1-1。 1 x 1 y 2 x 2 y 方向方向 的系数（秒的系数（秒/cm） AP1 AP2 BP2 CP1 CP2 P1 P2 -0.2244 0 0 0.2614 0 -0.0384 -0.3379 0 0 -0.2194 0 -0.4095 0 -0.7966 -0.2434 0 0.3864 0.0384 0 0.2572 -0.5528 0 0.084 0.4095 表表-1 000 jk

30、jkjk YXS、 kkjj yxyx, 65.2062 第八讲测量程序设计在测绘中应用 初步性 （三）（三）.计算测边的误差方程的系数。即计算测边的误差方程的系数。即 ik jk jk k jk jk j jk jk j jk jk i ly S Y x S X y S Y x S X v 0 0 0 0 0 0 0 0 的系数。的系数。 1 x 1 y 2 x 2 y 方向方向 边长误差方程系数边长误差方程系数 AP1 AP2 BP2 CP1 CP2 P1 P2 -0.833 0 -0.643 0 -0.996 -0.996 0.553 0 0 -0.766 0 0.093 0 0.307

31、 -0.915 0 0.212 0.996 0 0.952 0.403 0 -0.977 -0.093 表表-2 第八讲测量程序设计在测绘中应用 初步性 )( 00 jhjkii Ll 0 jkii SLl l （四）（四）.确定误差方程。在计算确定误差方程。在计算 或或 时，利用近似坐标求出每个边的近似坐标方位角和近似边长。时，利用近似坐标求出每个边的近似坐标方位角和近似边长。 从而可确定从而可确定 的系数。进一步可得出误差方程的系数项，见表的系数。进一步可得出误差方程的系数项，见表-3。 1 x 1 y 2 x 2 y lp i i S 表表-3 编号编号 1 1 2 2 3 3 4 4

32、5 5 6 6 7 7 8 8 9 9 1010 1111 1212 0 0 0 0 0 0 0 0 0 0 0 0 -0.2244-0.2244 0.03840.0384 0.1860.186 0.29980.2998 -0.0384-0.0384 -0.2614-0.2614 0 0 0 0 0 0 0 0 0 0 0 0 -0.3379-0.3379 0.40950.4095 0.07160.0716 0.19010.1901 -0.4095-0.4095 0.21940.2194 0.55320.5532 0.24340.2434 -0.7966-0.7966 -0.2434-0.24

33、34 0.62980.6298 -0.3864-0.3864 0.79660.7966 -0.835-0.835 0.03840.0384 -0.0384-0.0384 -0.348-0.348 0.3864 0.3864 -0.8055-0.8055 0.55280.5528 0.25720.2572 -0.5528-0.5528 0.63680.6368 -0.084-0.084 -0.2572-0.2572 -0.1523-0.1523 0.40950.4095 -0.4095-0.4095 0.32550.3255 0.0840.084 -2.6-2.6 -0.11-0.11 -0.0

34、2-0.02 0.480.48 0.300.30 0.420.42 0.020.02 0.050.05 -0.97-0.97 0.240.24 1.051.05 -5.09-5.09 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1313 1414 1515 1616 1717 1818 0 0 0 0 0 0 -0.643-0.643 -0.833-0.833 -0.996-0.996 0 0 0 0 0 0 -0.766-0.766 0.5530.553 0.0930.093 0.3070.307 -0.915-0.915 0.2120.2

35、12 0 0 0 0 0.9960.996 0.9520.952 0.4030.403 -0.977-0.977 0 0 0 0 -0.093-0.093 0.10.1 12.112.1 9.79.7 10.710.7 0 0 0.10.1 0.56250.5625 0.56250.5625 0.56250.5625 0.56250.5625 0.56250.5625 0.56250.5625 第八讲测量程序设计在测绘中应用 初步性 可得法方程为可得法方程为 （五）（五）.平差值计算与精度评定。由法方程可解得参数改正数为：平差值计算与精度评定。由法方程可解得参数改正数为： 2965. 0024

36、9. 00481. 00074. 0 2709. 00775. 01583. 0 9888. 00137. 0 7970. 0 , 8997. 0 3647. 3 8065. 6 5297. 3 1 2 2 1 1 对称 bb N y x y x 0 2148. 1 9415. 8 1620. 6 7469. 2 4194. 32866. 01437. 00229. 0 2866. 03174. 43363. 08608. 0 1437. 03363. 00459. 10861. 0 0229. 08608. 00861. 04270. 1 2 2 1 1 y x y x 第八讲测量程序设计在

37、测绘中应用 初步性 x lxBV 21,P P 由上式可解出由上式可解出，再根据，再根据 得各个改正数为得各个改正数为 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. 30 tn PVV T 2965. 0,2709. 0,9888. 0,7970. 0 2 2 2 2 1 1 1 1 YYXXYYXX QQQQ 2.4

38、833(cm) 4.1622(cm) 2 2 2 2 0 2 1 1 1 1 0 1 YYXX P YYXX P QQ QQ 第八讲测量程序设计在测绘中应用 初步性 利用利用Matlab编写上题的代码为：编写上题的代码为： clear B=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

39、 0.0716 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.

40、1. ;9.7;10.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/100 invN=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*

41、sqrt(invN(3,3)+invN(4,4) 第八讲测量程序设计在测绘中应用 初步性 v二、二、MATLABMATLAB在数字测量技术中简单应用在数字测量技术中简单应用 v1.1.坐标正、反算坐标正、反算 v1 1）坐标正算）坐标正算 v 坐标正算即已知坐标正算即已知A A点坐标点坐标 ,A ,A点至点至B B点水平点水平 距离距离 以及坐标方位角以及坐标方位角 ，求，求B B点坐标点坐标 的过程。计算公式如下：的过程。计算公式如下： v AB S AB ABAABABAB ABAABABAB yySyy xxSxx sin cos 第八讲测量程序设计在测绘中应用 初步性 vax=inpu

42、t(请输入请输入A的纵坐标的纵坐标:) vay=input(请输入请输入A点的横坐标点的横坐标:) vs=input(请输入请输入AB长度长度:) vad=input(请输入请输入AB的坐标方位角度的坐标方位角度:) 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两

43、已知点的坐标分别为两已知点的坐标分别为 和和 v，则直线，则直线ABAB的坐标方位角的坐标方位角 和水平距离和水平距离 为：为： v v ),( AA yx ),( BB yx 22 22 )()( arctanarctan ABABABABAB AB AB AB AB AB yxyyxxS x y xx yy 第八讲测量程序设计在测绘中应用 初步性 vclear vclc vx=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 on vfo

44、r i=1:5 v plot(x(i),y(i),o) vend vfor 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); end 第八讲测量程序设计在测绘中应用 初步性 vs %给出各边长给出各边长 vfor i=1:4 v if a(i)=0 m(i)=atand(abs(b(i)/a(i); v %结果是度，得到是与坐标纵轴结果是度，得到是与坐标纵轴x轴构成的锐角轴构成的锐角 k(i)=m(i); if a(i)0 elseif a(i)0 c(i)=180-k(i); elseif a(i)0 v el

45、se v c(i)=270; v end v end v d=fix(c) v e=fix(c-d)*60) v f=(c-d)*60-e)*60 v g=d e f v v 第八讲测量程序设计在测绘中应用 初步性 2.可以进行地形图绘制（若有可以进行地形图绘制（若有11个碎部点坐标）个碎部点坐标） v 115.14 v129.56 v138.66 v2112.55 v2210.51 v232.70 v3122.03 v3222.97 v3332.94 v4255.58 v4324.55 第八讲测量程序设计在测绘中应用 初步性 vfid=fopen(E:dengzhixian.txt,r);

46、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=linspace(yb,ye,200); vX,Y=meshgrid(Xl,Yl); %坐标矩阵坐标矩阵 第八讲测量程序设计在测绘中应用

47、初步性 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 off 第八讲测量程序设计在测绘中应用 初步性 v3.导线网近似平差计算导线网近似平差计算 v 规范规定，对于二级及以下等级的图根导线规范规定，对于二级

48、及以下等级的图根导线 允许对单一导线、单结点导线网采用近似平差方允许对单一导线、单结点导线网采用近似平差方 法进行计算。导线近似平差的基本思路是将角度法进行计算。导线近似平差的基本思路是将角度 观测值和边长观测值分别进行平差计算。观测值和边长观测值分别进行平差计算。 v附合导线坐标近似平差计算步骤如下：附合导线坐标近似平差计算步骤如下： v（1 1）方位角计算；）方位角计算； v（2 2）由转角公式计算其角度闭合差；）由转角公式计算其角度闭合差； v（3 3）将角度闭合差除以测边数分配到各观测角中；将角度闭合差除以测边数分配到各观测角中； v（4 4）计算调整后的导线边的方位角；）计算调整后的导线边的方位角； 第八讲测量程序设计在测绘中应用 初步性 v（5）计算坐标增量；）计算坐标增量； v（6）计算坐标增量闭合差并分配；）计算坐标增量闭合差并分配； v（7）用调整后的方位角计算各点坐标。）用调整后的方位角计算各点坐标。 v对于闭合导线和支导线的情况可类似计算对于闭合导线和支导线的情况可类似计算 v注意：在进行控制网精平时，一般采用间接注意：在进行控制网精平时，一般采用间接 平差，误差方程形式统一、规律性强、便于平差，误差方程形式统一、规律性强、便于 程序开发；而使用条件平差等，其方程形式程序开发；而使用条件平差等，其方程形式 多样、规律性差，不利于程序