攝影測量學(xué)單像空間后方交會編程實(shí)習(xí)報告

1、攝影測量學(xué)單像空間后方交會編程實(shí)習(xí)報告班級： 130x姓名： XX學(xué)號： 2013302590xxx 指導(dǎo)老師：荻一、實(shí)習(xí)目的通過對提供的數(shù)據(jù)進(jìn)行計算，輸出像片的外方位元素并評定精度。深入理解 單像空間后方交會的思想，體會在有多余觀測情況下，用最小二乘平差方法編程 實(shí)現(xiàn)解求影像外方位元素的過程。通過嘗試編程實(shí)現(xiàn)加強(qiáng)編程處理問題的能力和 對實(shí)習(xí)內(nèi)容的理解，通過對實(shí)驗(yàn)結(jié)果的分析，增強(qiáng)綜合運(yùn)用所學(xué)知識解決實(shí)際問 題的能力。了解攝影測最平差的基本過程，學(xué)握空間后方交會的定義和實(shí)現(xiàn)算法。實(shí)習(xí)內(nèi)容根據(jù)學(xué)習(xí)的單像空間后方交會的知識，用程序設(shè)計語言（C 卄或 C 語言）編寫一個完整的單像空間后方交會程序，通過

2、對提供的數(shù)據(jù)進(jìn)行計算，輸出像片的外方位元素并評定精度。三、實(shí)習(xí)數(shù)據(jù)已知航攝儀的內(nèi)方位元素: 2153. 24mm, xo=yo=0,攝影比例尺為1:15000 ：4 個地而控制點(diǎn)的地面坐標(biāo)及其對應(yīng)像點(diǎn)的像片坐標(biāo)：!1!實(shí)習(xí)原理像片坐標(biāo)地面點(diǎn)坐標(biāo)x(mm)y(mm)X(m)Y(m)Z(m)1?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如果我們知道每幅影像的6 個外方位

3、元素，就能確定被攝物體與航攝影像的關(guān)系。因此，如何獲取影像的外方位元素，一直是攝影測量匸作者所探討的問題?？刹扇〉姆椒ㄓ校豪美走_(dá)、全球定位系統(tǒng)（GPS ） 、慣性導(dǎo)航系統(tǒng) （ INS）以及星相攝影機(jī)來獲取影像的外方位元素；也可以利用影像覆蓋范偉1 內(nèi)一定數(shù)量的控制點(diǎn)的空間坐標(biāo)與攝影坐標(biāo)，根據(jù)共線條件方程，反求該影像的外方位元素，這種方法稱為單幅影像的空間后方交會。單像空間后方交會的基本思想是：以單幅影像為基礎(chǔ)，從該影像所覆蓋地面范用內(nèi)若干控制點(diǎn)的己知地而坐標(biāo)和相應(yīng)點(diǎn)的像坐標(biāo)量測值出發(fā)，根據(jù)共線條件方程，解求該影像在航空攝影時刻的外方位元素Xs, Ys, Zs,OUK 。五、實(shí)習(xí)流程1.獲取已

4、知數(shù)據(jù)。從攝影資料中查取影像比例尺1/m, 平均攝影距離 ( 航空攝影的航高、內(nèi)方位元素xo, y 0, f ;獲取控制點(diǎn)的空間坐標(biāo)人，X， Zto2.量測控制點(diǎn)的像點(diǎn)坐標(biāo)并進(jìn)行必要的影像坐標(biāo)系統(tǒng)誤差改正，得到像點(diǎn)坐標(biāo)。3.確定未知數(shù)的初始值。單像空間后方交會必須給出待定參數(shù)的初始值，在豎直航空攝影且地面控制點(diǎn)大體對稱分布的情況下，可按如下方法確定初始值：Z o = H = m ? fni-1總i-1滬=3° = 0式中： m 為攝影比例尺分母, n 為控制點(diǎn)個數(shù)；4.計算旋轉(zhuǎn)矩陣Ro 利用角元素的近似值按公式計算方向余弦值ah ac, as, bi, b：,b3 , c b C2,

5、 cs, 組成R 陣。5.逐點(diǎn)計算像點(diǎn)坐標(biāo)的近似值。利用耒知數(shù)的近似值按共線條件方程計算控制點(diǎn)像點(diǎn)坐標(biāo)的近似值(x), (y)o6. 逐點(diǎn)計算誤差方程式的系數(shù)和常數(shù)項(xiàng)，組成誤差方程。7. 計算法方程的系數(shù)矩陣 A% 與常數(shù)項(xiàng) AI, 組成法方程。8.解求外方位元素。根據(jù)法方程，解求外方位元素的改正數(shù)，并與相應(yīng)的近似值求和，得到外方位元素新的近似值。9.檢查計算是否收斂。將所求得的外方位元素的改正數(shù)與規(guī)定的限差比較,通常對 ?U),K的改正數(shù) q)4 34K 給予限差，當(dāng)改正數(shù)小于限差時，迭代結(jié)束。否則用新的近似值重復(fù) (4) 一一(8) 步驟計算，直到滿足要求為止。10. 空間后方交會的精度估

6、計：按上述方法所求得的影像外方位元素的精度可以通過法方程式中未知數(shù)的系數(shù)矩陣的逆陣(A%) 來解求，此時視像點(diǎn)坐標(biāo)為等精度不相關(guān)觀測值。因 (AT A)- 1 中第 i 個主對角線上的元素就是法方程式中第 i 個未知數(shù)的權(quán)倒數(shù)， 若單位權(quán)中誤差為 mo, 則第 i 個未知數(shù)的中誤差為：=何 mo當(dāng)參加空間后方交會的控制點(diǎn)有n 個時，則單位權(quán)中誤差可按下式計算:m0 =±w2n-6流程圖如下：整理并輸岀計算結(jié)果非正常結(jié)束正常結(jié)束# include <stdio.h> include -Matrix-h 11矩陣運(yùn)算頭文件，來自網(wǎng)絡(luò)，己上傳，外方位元素double Xs,Ys

7、 Zsqwkdouble xOOj;內(nèi)方位元素double x4,y4;影像坐標(biāo)double X04“ 04,Z04;地而坐標(biāo)double m;比例尺double a2,a2,a3,bl,b2,b3，crc2? c3; 旋轉(zhuǎn)矩陣int n=0;迭代次數(shù)地面坐標(biāo)X00=36589.41;X0l=37631.08;X02=39100.97;X03=40426.54;Y00=25273.32;Y0 =31324.51;Y02=24934.98;Y03=30319.81;Z00=2195.17;Z0l=728.69;Z02=2386.50;Z03=75731;影像坐標(biāo)刈 0)=0.08615;x =0

8、05340;x2=-0.01478;x3=0.01046;y0=-0.06899;yl=0.08221;y2=-0.07663;丫 =0.06443;確定內(nèi)外方位元素初始(fix0=00=0;f=153.24/1000;m=15000;Xs=O;Ys=O;q=03A/=0；k=0;for (int i=0;i<4;i+)Xs += XOi;Ys += YOi; Xs /= 4;Ys /= 4;Zs=f*m;Matrix A(& 6);Matrix L(& 1);Matrix AT(6, 8);Matrix ATA(6# 6);Matrix ATAJ6, 6);Matrix

9、 Xx(6,1);迭代計算do旋轉(zhuǎn)矩陣Ral = cos(q)*cos(k) - sin(q)*sin(w)*sin(k);a2 = -cos(q)*sin(k) - sin(q)*sin(w)*cos(k);a3 = -sin(q)*cos(w);bl = cos(w)*sin(k);b2 = cos(w)*cos 伙 ) ；b3 = -sin(w);cl = sin(q)*cos(k) + cos(q)wsin(w)*sin(k);c2 = -sin(q)*sin(k) + cos(q)*sin(w)*cos(k);計算像點(diǎn)坐標(biāo)for (int i = 0; i < 4; i+)do

10、uble X = al * ( X0(i - Xs) + bl * (Y0i - Ys) + cl * (Z0i? Zs); double Y = a2* ( X0i - Xs) + b2 * (Y0i - Ys) + c2 * (Z0l - Zs);double Z = a3 ? ( X0i - Xs) + b3 * (Y0i - Ys) + c3 * (Z0i? Zs); double xx =xi-x0; double yy = yi-yO;A2 ? i0 = (al*f + a3*xx) / Z;A2*il = (bl*f + b3*xx)/Z;A2#i2 = (cl*f + c3*x

11、x)/Z;A2 ? i3 = yy*sin(w) - (xx*(xx*cos(k) - w*sin(k) / f + f*cos(k)*cos(w); A2 * i4 = -f*sin(k) - xx*(xx*sin(k) + yy*cos(k) / f;A2*i5=yy;A2* i +10 = (a2*f + a3*yy) / Z;A2 拿 i + 11 = (b2*f + b3*yy) / Z;A2# i + 12 = (c2 卄 +<3 拿 yy) / Z;A2 * i + 13 = - xxsin(w) ? (yy*(xx*cos(k) -yysin(k) /f -f*sin(k

12、)*cos(w);A2 ? i + 14 = -f*cos(k) -yy*(xx*sin(k) + yy*cos(k) /f;A2 # i + 1二? xx;L2*i0)=xi-(-f*X/Z);L2*i + l0 =yi-(-f*Y/Z);法方程的解X=(ATA)_ATL AT = A.Transpose();ATA=AT*A;轉(zhuǎn)豐ATA_ = ATA.Inverse();Xx = (ATA_ * AT) * L;求逆法方程解Xs +=XxOO; Ys += XxlO;Zs +=Xx20; q+=Xx0; w += Xx40; k +=Xx50;/printffdXnn);H%.5lf%.5

13、lfn,l/z/printf( %lf %lf %lf %? 5lfXs, Ys, Zs, qw, k);while <(abs(Xx30) > 0.000001 | abs(Xx40) > 0.000001 | abs(Xx50) > 0.000001) && n<100);求中誤差double m0,mi6;單位權(quán)中謀差的絕對值mO以及第i 個未知數(shù)的中i 靈差mimO = sqrt(A*Xx - L).Transpose()*(A*Xx ? L).ToDouble() / (2 * 4 - 6); Matrix Q(6,6);Q = (A.T

14、ranspose() A).Inverse();for (int i = 0; i < 6; i+)mii = (sqrt(Qii)*mO;/printf("96lfn"/mO);輸出printf ( “迭代次數(shù) dn Hzn);printf ( ”旋轉(zhuǎn)矩陣n“)；printf(%5lf %.5lf%? 5lfn ：a：La 乙 a3);H,b2 b3);printf( %.5lf %.5lf%.5lfn bl/printfC'%.5lf %.5lf%.5lfn" / cl,c2,c3);printf("n 外方位元索解 :nWs=%.2l

15、f Ys=%.2lf Zs=%.2lfnq=%.5lf k=%.5lfnn". Xs, Ys, Zs,w=%.5lfq, w, k);printf ( ”單位權(quán)中誤差的絕對值:%n” ,mO);printfC'Xs 的精度 :%printff'Ys 的精度 56lfmn",mil);printffZs 的精度 56lfmn ,/mi2);printff'q 的精度 :%3);printf( M w 的精度 lfn ,mi4);printf("k 的精度) ；保存到文件，結(jié)果文件默認(rèn)保存在程序冃錄FILE*fp;fp=fopen("

16、resulttxt","w");/ fprintfffp ；'迭代次數(shù) :%n“ ,n);fprintf(fp, M 旋轉(zhuǎn)矩陣R:rT);fprintf(fp,"%.5lf %.5lf%.5lfn",al,a2,a3);fprintf(fp, .5lf”%.5lfn” ,bl,b2,b3);fprintf(fp z"%.5lf%.5lf%.5lfn ,/d/c2/ c3);w=%.5lffprintfffp/'Xn 外方位元素解 : nWs=%.2lfYs=%.2lf Zs=%.2lfnq=%.5lfk=%.5lfnn

17、". Xs, Ys, Zs, q, w, k);fprintfffp ；* 單位權(quán)中誤差的絕對值S6lfmn",mO);fprintfffp,"Xs 的精度 : lfmn",miO);fprintf(fp/'Ys 的精度 56lfmn" zmil);fprintf(fp,"Zs 的精度僥fprintf(fp ；'q 的精度 :%3);fprintf(fp/'w的精度 56lfn",mi4);fprintf(fp, M k 的精度 56lfn H,mi5);七、運(yùn)算結(jié)果1.運(yùn)行結(jié)果2. 文件保存結(jié)果 resuk.txt ?記事本文件 (F) 瞬( E)