攝影測(cè)量實(shí)習(xí)報(bào)告

1、實(shí)驗(yàn)報(bào)告一：?jiǎn)蜗窨臻g后方交會(huì)程序設(shè)計(jì)報(bào)告內(nèi)容：實(shí)驗(yàn)報(bào)告主要是針對(duì)實(shí)驗(yàn)?zāi)康暮蛯?shí)驗(yàn)過程進(jìn)行描述，并對(duì)所得實(shí)驗(yàn) 結(jié)果進(jìn)行分析和評(píng)價(jià)，指出實(shí)驗(yàn)注意事項(xiàng)和對(duì)本實(shí)驗(yàn)提出建議。(一)已知數(shù)據(jù)與作業(yè)要求已知數(shù)據(jù)：XyXYZ1-86.15-68.9936589.4125273.322195.172-53.4082.2137631.0831324.51728.693-14.78-76.6339100.9724934.982386.50410.4664.4340426.5430319.81757.31內(nèi)方位元素？?= 153.24，?= ?= 0 ;攝影比例尺分母 m=4000Q要求：已知4對(duì)點(diǎn)的影像坐標(biāo)和地面坐標(biāo)

2、，計(jì)算近似垂直下空間后方交會(huì)的 解，得到6個(gè)外方位元素的值。(二) 單像空間后方交會(huì)計(jì)算原理空間后方交會(huì)的數(shù)學(xué)模型是共線條件方程，其具體形式可寫為：?(? ? + ?(? ? + ?(? ?(?- ?+ ?(? ? + ?(? ?)?(?- ? + ?(? ? + ?込？2 ?= -? ?(?- ? + ?(?- ? + ?- ?式中(?為像點(diǎn)的像平面坐標(biāo)，(?？?為像點(diǎn)對(duì)應(yīng)地面點(diǎn)的物方空間坐標(biāo)，？為相 片主距，？?為外方位線元素，？?= 1,2,3)為相片外方位角元素所組成的 9個(gè)方向余弦。對(duì)任一控制點(diǎn)，其地面坐標(biāo)(?？?穢和對(duì)應(yīng)像點(diǎn)坐標(biāo)(?都是已知的，代入 共線條件方程可以列出兩個(gè)方程式，

3、因此從理論上講只要有三個(gè)控制點(diǎn)就可以列出 6個(gè)方程，從而求解出6個(gè)外方位元素。但為了避免粗差，提高測(cè)量精度，應(yīng)有多 余觀測(cè)，因此，在實(shí)際應(yīng)用中一般需要4個(gè)或更多的控制點(diǎn)。(三) 單像空間后方交會(huì)的計(jì)算過程1. 獲取已知數(shù)據(jù)：從航攝資料中獲取像片比例尺 1/m、內(nèi)方位元素？?,？?，獲 取控制點(diǎn)的空間坐標(biāo)X，丫，Zo2. 量測(cè)控制點(diǎn)的像點(diǎn)坐標(biāo)：量測(cè)控制點(diǎn)框標(biāo)坐標(biāo)，并經(jīng)像主點(diǎn)改正得到像點(diǎn)坐 標(biāo)x，y。3. 確定未知數(shù)的初始值：在豎直攝影情況下，角元素的初始值給0,線元素中?0= ? ?Q ?可取四個(gè)角上控制點(diǎn)坐標(biāo)的平均值?=韋？ ?= ?= ? = 7?1 = ? = ORO式中：仃為攝影比例尺分

4、母；n為控制點(diǎn)個(gè)數(shù)4. 用三個(gè)角元素的初始值計(jì)算方向余弦，并組成旋轉(zhuǎn)矩陣?=:cos?cos? Sin?sin?Sin?=-COS ?si n ? sin? ?Si n ? ?cos ?=-Sin ?cos? ? ?=cos?Sin ?= ? ? ? ?=cos?cos ? ?2? ?=-Sin ?1?=Sin ?cos? cos?sin ?sin ?=-Sin ?Sin ?+ cos?sin ?cos?=cos?cos?5.用所取未知數(shù)的初始值和控制點(diǎn)的地面坐標(biāo)，代入共線方程，逐點(diǎn)計(jì)算像點(diǎn)(x)、 (y)。坐標(biāo)的近似值,?=?= ?(?-?)+?(?-?+?iX?-?_?(?-?)+?(?-

5、?+?3(?-?)“？?(？?-?)+?(?-?+?2(?-?)(?)=(?)=(?2 =(?2 =(?3 =(?3 =(?4 =(?4 =“？?(？- ? + ?(?- ? + ?(?- ? ?(?- ? + ?(?- ? + ?(?- ? “？(？- ? + ?(?-鋼 + ?(?- ? ?(?- ? + ?(?-漏 + ?(?- ? ?(?- ? + ?(?- ? + ?(?- ? -?(? - ? + ?(?-鬲 + ?孔？- E ?(?- ? + ?(?-鋼 + ?(?- ? _?(? - ? + ?(?-鋼 + ?亂？- ? _?(?- ? + ?(?-錮 + ?(?- ? _?(

6、? - ? + ?(?-鋼 + ?亂？- ?) “？(？- ? + ?(?-鋼 + ?(?- ? ?(? - ? + ?(?-鋼 + ?(?- ? ?(?- ? + ?(?-窮 + ?(?- ? -?(? - ? + ?(?- ?) + ?- E _?(?- ? + ?(?-笏 + ?(?- ? _?(? - ? + ?(?-鋼 + ?亂？- ?6. 用像點(diǎn)的觀測(cè)值和由5中計(jì)算的近似值，計(jì)算每個(gè)點(diǎn)的常數(shù)項(xiàng) ？ ?+?+ ?(?初？刃+?(?衿？?)+?寅？務(wù)？?)? ?(?-?)+?(?-?+?3<?7?(?2123 4)?(?7?)+?(?*?+?(?7?)(?，3,?(?S-?)+?

7、(?-?+?(?)7. 計(jì)算法方程的系數(shù)矩陣？?與常數(shù)項(xiàng)？?(?-?)+?(?-?)+?3(?-?)? Z= -?2?= 01?=-?4 =:-?(1+祠?5 =?1?=?1 =0?=- ?=-?4 =Z -?5 =-?(1 +Tr)?6 =-?I6X 08 ?2?3?4?5?6?2?3?4?5?6?2?3?4?5?6?2?3?4?5?46?2?3?5?6?2?3?4?5?6?2?3?4?5?6?2?3?4?5?6 ?1?1?1? ? ?礒 000?注意：A矩陣中的（?不是已知的像點(diǎn)坐標(biāo)，而是將每次迭代后所得到的外方 位元素代入共線條件方程計(jì)算出來的，A矩陣在迭代過程中也是變化的8. 求解法方

8、程，得到未知數(shù)的改正數(shù)。?= （?）1 ?9. 檢查計(jì)算結(jié)果是否收斂：將改正數(shù)與限差比較，小于限差則計(jì)算終止，否則 用新的近似值重復(fù)4-9的計(jì)算，直到滿足要求為止。?= ?M+?需?=T?"1+ ?=?-1+?勇=?畀-1+?液? =?參-1+ ?層? =?-1+?式中：？?弋表迭代次數(shù)。10. 空間后方教會(huì)的精度估算。 ?= ?1?+?2 ?3?+?4 ?+?5 ?5 K ? - ? = ?1?22 ?+?3?+?4 ?+?25 ?26 ? - ?由單位權(quán)中誤差?= ?和6個(gè)參數(shù)的協(xié)因數(shù)陣？?= (?)1得出參數(shù)？?的中誤差為？?= ?丿?爲(wèi)？（四）算法流程圖（五）計(jì)算結(jié)果MATL

9、A計(jì)算結(jié)果輸出截圖:JCs =3.979a44340D520e÷C 4fai =0. W39S5559117g04¥3 =N 7476464S4D27101e4C4Onega OJ 002Ll3670028136ZS =7. 57368833LD23273t÷C3kaa -OJ 065r64B10370,5IJlII -R =1.125401304950231L 2436735942644O- 6377 &56790616O. 9770900f98137C. 0675j406U65390. 003985629934395-0. C6526C4663190.

10、 957715272B03382-0. 0021137265110190. (MQlaISO7341 L-QI 004119272433176Ol 001839750147646OJ 999989823405350Ol 0001596811695 詬OL DOOo"243527C#計(jì)算結(jié)果輸出截圖（六）實(shí)驗(yàn)心得與體會(huì)（七）程序?qū)崿F(xiàn)代碼見附錄附錄：?jiǎn)蜗窨臻g后方交會(huì)代碼 matlab 程序 : clc;clear;format long%單位統(tǒng)一為米xy=1/1000*-86.15 -68.99;-53.40 82.21;-14.78 -76.63;10.46 64.43;% 導(dǎo)入像點(diǎn)坐

11、標(biāo)XYZ=36589.41 25273.32 2195.17;37631.08 31324.51 728.69;.39100.97 24934.98 2386.50;40426.54 30319.81 757.31; % 導(dǎo)入物點(diǎn)坐標(biāo) f=0.15324;x0=0;y0=0;m=40000; %攝影比例尺分母fai=0;omega=0;kappa=0; % 初始外方位角元素Xs=0;Ys=0;for i0=1:4;Xs=Xs+XYZ(i0,1);Ys=Ys+XYZ(i0,2);endZs=m*f;Xs=Xs/4;Ys=Ys/4; %初始外方位線元素R=zeros(3);L=zeros(8,1)

12、;delta=zeros(6,1);A=zeros(8,6);x=zeros(4,1);y=zeros(4,1);t=0;while 1%旋轉(zhuǎn)矩陣 R(1,1)=cos(fai)*cos(kappa)-sin(fai)*sin(omega)*sin(kappa);R(1,2)=-cos(fai)*sin(kappa)-sin(fai)*sin(omega)*cos(kappa); R(1,3)=-sin(fai)*cos(omega);R(2,1)=cos(omega)*sin(kappa);R(2,2)=cos(omega)*cos(kappa);R(2,3)=-sin(omega); R(

13、3,1)=sin(fai)*cos(kappa)+cos(fai)*sin(omega)*sin(kappa); R(3,2)=-sin(fai)*sin(kappa)+cos(fai)*sin(omega)*cos(kappa); R(3,3)=cos(fai)*cos(omega);for i1=1:4 %共線條件方程 x(i1,1)=-f*(R(1,1)*(XYZ(i1,1)-Xs)+R(2,1)*(XYZ(i1,2)-Ys)+R(3,1)*(XYZ(i1,3)-Zs)/. (R(1,3)*(XYZ(i1,1)-Xs)+R(2,3)*(XYZ(i1,2)-Ys)+R(3,3)*(XYZ(

14、i1,3)-Zs);y(i1,1)=-f*(R(1,2)*(XYZ(i1,1)-Xs)+R(2,2)*(XYZ(i1,2)-Ys)+R(3,2)*(XYZ(i1,3)- Zs)/.(R(1,3)*(XYZ(i1,1)-Xs)+R(2,3)*(XYZ(i1,2)-Ys)+R(3,3)*(XYZ(i1,3)-Zs); %誤差方程常數(shù)項(xiàng)L(2*i1-1,1)=xy(i1,1)-x(i1,1); L(2*i1,1)=xy(i1,2)-y(i1,1);%系數(shù)矩陣 A(2*i1-1,1)=-f/(Zs-XYZ(i1,3);A(2*i1-1,2)=0; A(2*i1-1,3)=-x(i1,1)/(Zs-XY

15、Z(i1,3);A(2*i1-1,4)=-f*(1+x(i1,1)*x(i1,1)f2);A(2*i1-1,5)=-x(i1,1)*y(i1,1)/f;A(2*i1-1,6)=y(i1,1);A(2*i1,1)=0; A(2*i1,2)=-f(Zs-XYZ(i1,3);A(2*i1,3)=-y(i1,1)(Zs-XYZ(i1,3);A(2*i1,4)=-x(i1,1)*y(i1,1)f; A(2*i1,5)=-f*(1+y(i1,1)*y(i1,1)f2);A(2*i1,6)=-x(i1,1);enddelta=(A'*A)A'*L; % 改正數(shù)Xs=Xs+delta(1,1)

16、;Ys=Ys+delta(2,1);Zs=Zs+delta(3,1);fai=fai+delta(4,1);omega=omega+delta(5,1); kappa=kappa+delta(6,1);t=t+1;if (abs(delta(4,1)<10(-6)&&(abs(delta(5,1)<10(-6)&&( abs(delta(6,1)<10(- 6)breakendendV=A*delta-L;VV=0;for i2=1:8VV=VV+V(i2)*V(i2);endm0=sqrt(VV2); %單位權(quán)中誤差mi=zeros(6,1)

17、; %未知數(shù)中誤差 ffcinv=inv(A'*A); %法方程的逆for i3=1:6 mi(i3)=m0*sqrt(ffcinv(i3,i3);end %輸出外方位元素，旋轉(zhuǎn)矩陣和未知數(shù)中誤差XsYsZs fai omega kappa R miC#程序：using System;using System.Collections.Generic;using System.Linq;using System.Text;using System.Threading.Tasks;namespace resectionclass Program/轉(zhuǎn)置static double, Trans

18、pose(double, a)double, b = new doublea.GetLength(1), a.GetLength(0);for (int i = 0; i < a.GetLength(0); i+)for (int j = 0; j < a.GetLength(1); j+)bj, i = ai, j;return b;/相加static double, addition(double, a, double, b)double, c = new doublea.GetLength(0), a.GetLength(1);if (a.GetLength(0) = b.G

19、etLength(0) && (a.GetLength(1) = b.GetLength(1)for (int i = 0; i < a.GetLength(0); i+)for (int j = 0; j < a.GetLength(1); j+)ci, j = ai, j + bi, j;return c;elseSystem.Exception e = new Exception(" 兩矩陣維數(shù)不等，檢查數(shù)據(jù)");throw e;/相減static double, subtraction(double, a, double, b)doubl

20、e, c = new doublea.GetLength(0), a.GetLength(1);if (a.GetLength(0) = b.GetLength(0) && (a.GetLength(1) = b.GetLength(1)for (int i = 0; i < a.GetLength(0); i+)for (int j = 0; j < a.GetLength(1); j+)ci, j = ai, j - bi, j;return c;elseSystem.Exception e = new Exception(" 兩矩陣維數(shù)不等，檢查數(shù)

21、據(jù)");throw e;/相乘static double, multiply(double, a, double, b)double, c = new doublea.GetLength(0), b.GetLength(1);if (a.GetLength(1) != b.GetLength(0)System.Exception e = new Exception(" 左矩陣列數(shù)與右矩陣行數(shù) 不等，檢查數(shù)據(jù) ");throw e;elsedouble temp;for (int i = 0; i < a.GetLength(0); i+)for (int j

22、 = 0; j < b.GetLength(1); j+)temp = 0;for (int k = 0; k < a.GetLength(1); k+)temp += ai, k * bk, j;ci, j = temp;return c;/行交換static double, swap(int temp, int temp1, double, a)double b=new doublea.GetLength(1);double c=new doublea.GetLength(1);for (int i = 0; i < a.GetLength(1); i+)bi = ate

23、mp, i;ci = atemp1, i;for (int i = 0; i < a.GetLength(1); i+)atemp1, i=bi; atemp, i= ci;return a;/行列式static double Det(double, a)double c = 1;if (a.GetLength(0) != a.GetLength(1)據(jù)");System.Exception e = new Exception(" 矩陣行列數(shù)不等，檢查數(shù) throw e;for (int k = 0; k < a.GetLength(0); k+)for (in

24、t i = k + 1; i < a.GetLength(0); i+) double temp2; if (ai, k != 0) temp2 = ai, k / ak, k; for (int j = 0; j < a.GetLength(1); j+) ai, j -= ak, j * temp2;/上三角矩陣for (int i = 0; i < a.GetLength(0); i+)c *= ai, i;return c;/ 求逆static double, inverse(double, a)double, b = new doublea.GetLength(0)

25、, 2 * a.GetLength(0); for (int i = 0; i < a.GetLength(0); i+)for (int j = 0; j < a.GetLength(0); j+)bi, j = ai, j;for (int j = a.GetLength(0); j < 2 * a.GetLength(0); j+)if (j - a.GetLength(0) = i)bi, j = 1;elsebi, j = 0;/以上步驟為構(gòu)造 (A,E) 矩陣double temp = Math.Abs(b0, 0);int temp1 = 0;for (int

26、i = 0; i < a.GetLength(0); i+)if (temp <= Math.Abs(bi, 0)temp = bi, 0;temp1 = i;/以上選出第一列中最大值b = swap(0, temp1, b);for (int k = 0; k < b.GetLength(0); k+)for (int i = k + 1; i < b.GetLength(0); i+)double temp2;if (bi, k != 0)temp2 = bi, k / bk, k;for (int j = 0; j < b.GetLength(1); j+)

27、 bi, j -= bk, j * temp2; /上三角矩陣 double c = Det(a); if (c = 0) System.Exception e = new Exception(" 矩陣行列式為零，檢查數(shù)據(jù)");throw e;/檢查矩陣行列式是否為 0for (int k = b.GetLength(0) - 1; k >= 0; k-)for (int i = k - 1; i >= 0; i-)double temp3;if (bi, k != 0)temp3 = bi, k / bk, k;for (int j = 0; j < b

28、.GetLength(1); j+) bi, j -= bk, j * temp3;/已處理成對(duì)角線矩陣for (int i = 0; i < b.GetLength(0); i+)double temp4 = bi, i;for (int j = 0; j < b.GetLength(1); j+)bi, j /= temp4;/已處理成 (E,A 逆 )的形式for (int i = 0; i < b.GetLength(0); i+)for (int j = b.GetLength(0); j < b.GetLength(1); j+)ai, j - b.GetL

29、ength(0) = bi, j;/此時(shí) a 為原始輸入矩陣的逆 return a;static void Main(string args)double, xy = new double, -86.15, -68.99 , -53.40, 82.21 , - 14.78, -76.63 , 10.46, 64.43 ;double, XYZ = new double, 36589.41, 25273.32, 2195.17 , 37631.08, 31324.51, 728.69 , 39100.97, 24934.98, 2386.50 , 40426.54, 30319.81, 757.

30、31 ;double f = 0.15324, m = 40000;double fai = 0, omega = 0, kappa = 0;/ 角元素 double Xs = 0, Ys = 0, Zs = 0;/ 線元素 double, R = new double3, 3;/ 旋轉(zhuǎn)矩陣 double, L = new double2*xy.GetLength(0), 1;/ 常數(shù)項(xiàng) double, A = new double2 * xy.GetLength(0), 6;/ 系數(shù)矩陣 double, x = new doublexy.GetLength(0), 1; double, y

31、= new doublexy.GetLength(0), 1; double, Diff = new double6,1;/ 改正數(shù)矩陣 double, V = new double2 * xy.GetLength(0),1;/ 誤差矩陣 double, FFCinv = new double6, 6;/ 法方程的逆矩陣 double, mi = new double6, 1;/ 未知數(shù)中誤差 double vv = 0,m0=0;/vv 為殘差平方和 m0 為單位權(quán)中誤差 int t = 0;/ 迭代次數(shù) for (int i = 0; i < 4; i+) for (int j =

32、0; j < 2; j+)xyi, j /= 1000;for (int i = 0; i < 4; i+)Xs = Xs + XYZi, 0;Ys = Ys + XYZi, 1;Zs = m * f;Xs = Xs / 4;Ys = Ys / 4; while (true) /旋轉(zhuǎn)矩陣R0, 0 = Math.Cos(fai) * Math.Cos(kappa) - Math.Sin(fai) * Math.Sin(omega) * Math.Sin(kappa);R0, 1 = -Math.Cos(fai) * Math.Sin(kappa) - Math.Sin(fai) M

33、ath.Sin(omega) * Math.Cos(kappa);R0, 2 = -Math.Sin(fai) * Math.Cos(omega);R1, 0 = Math.Cos(omega) * Math.Sin(kappa);R1, 1 = Math.Cos(omega) * Math.Cos(kappa);R1, 2 = -Math.Sin(omega);R2, 0 = Math.Sin(fai) * Math.Cos(kappa) + Math.Cos(fai) Math.Sin(omega) * Math.Sin(kappa);R2, 1 = -Math.Sin(fai) * Ma

34、th.Sin(kappa) + Math.Cos(fai) Math.Sin(omega) * Math.Cos(kappa);R2, 2 = Math.Cos(fai) * Math.Cos(omega);for (int i = 0; i < 4; i+) /共線條件方程 xi, 0 = -f * (R0, 0 * (XYZi, 0 - Xs) + R1, 0 * (XYZi, 1 - Ys) + R2, 0 * (XYZi, 2 - Zs) /(R0, 2 * (XYZi, 0 - Xs) + R1, 2 * (XYZi, 1 - Ys) + R2, 2 * (XYZi, 2 -

35、Zs);yi, 0 = -f * (R0, 1 * (XYZi, 0 - Xs) + R1, 1 * (XYZi, 1 - Ys) + R2, 1 * (XYZi, 2 - Zs) /(R0, 2 * (XYZi, 0 - Xs) + R1, 2 * (XYZi, 1 - Ys) + R2, 2 * (XYZi, 2 - Zs);/誤差方程常數(shù)項(xiàng) L2 * i, 0 = xyi, 0 - xi, 0;L2 * i+1, 0 = xyi, 1 - yi, 0;/系數(shù)矩陣A2 * i, 0 = -f / (Zs - XYZi, 2);A2 * i, 1 = 0;A2 * i, 2 = -xi, 0 / (Zs - XYZi, 2);A2 * i, 3 = -f * (1 + xi, 0 * xi, 0 / (f * f);A2 * i, 4 = -xi, 0 * yi, 0 / f;A2 * i, 5 = yi, 0;A2 * i+1, 0 = 0;A2 * i+1, 1 = -f / (Zs - XYZi