空間后方交會編程實習(xí)報告

1、空間后方交會編程實習(xí)報告 一 實習(xí)目的用程序設(shè)計語言（Visual C+或者C語言）編寫一個完整的單片空間后方交會程序，通過對提供的試驗數(shù)據(jù)進(jìn)行計算，輸出像片的外方位元素并評定精度。本實驗的目的在于讓學(xué)生深入理解單片空間后方交會的原理，體會在有多余觀測情況下，用最小二乘平差方法編程實現(xiàn)解求影像外方位元素的過程。通過上機(jī)調(diào)試程序加強(qiáng)動手能力的培養(yǎng)，通過對實驗結(jié)果的分析，增強(qiáng)學(xué)生綜合運用所學(xué)知識解決實際問題的能力。二 實習(xí)內(nèi)容 利用一定數(shù)量的地面控制點，根據(jù)共線條件方程求解像片外方位元素。 三 實習(xí)數(shù)據(jù) 已知航攝儀的內(nèi)方位元素：fk=153.24mm，x0=y0=0.0mm，攝影比例尺為1:500

2、00； 4個地面控制點的地面坐標(biāo)及其對應(yīng)像點的像片坐標(biāo)：點號 像片坐標(biāo)(mm) 地面點坐標(biāo)(m) x y X Y Z 1 -86.15 -68.99 36589.41 25273.32 2195.17 2 -53.40 82.21 37631.08 31324.51 728.69 3 -14.78 -76.63 39100.97 24934.98 2386.50 4 10.46 64.43 40426.54 30319.81 757.31 四 實習(xí)原理 如果我們知道每幅影像的6個外方位元素，就能確定被攝物體與航攝影像的關(guān)系。因此，如何獲取影像的外方位元素，一直是攝影測量工作者所探討的問題。可采

3、取的方法有：利用雷達(dá)、全球定位系統(tǒng)（GPS）、慣性導(dǎo)航系統(tǒng)（INS）以及星相攝影機(jī)來獲取影像的外方位元素；也可以利用影像覆蓋范圍內(nèi)一定數(shù)量的控制點的空間坐標(biāo)與攝影坐標(biāo)，根據(jù)共線條件方程，反求該影像的外方位元素，這種方法稱為單幅影像的空間后方交會。單像空間后方交會的基本思想是：以單幅影像為基礎(chǔ)，從該影像所覆蓋地面范圍內(nèi)若干控制點的已知地面坐標(biāo)和相應(yīng)點的像坐標(biāo)量測值出發(fā)，根據(jù)共線條件方程，解求該影像在航空攝影時刻的外方位元素Xs，Ys，Zs，t，w，k。五 實習(xí)流程（1）獲取已知數(shù)據(jù)。從攝影資料中查取影像比例尺1/m，平均攝影距離（航空攝影的航高、內(nèi)方位元素x0，y0，f；獲取控制點的空間坐標(biāo)Xt

4、，Yt，Zt。（2）量測控制點的像點坐標(biāo)并進(jìn)行必要的影像坐標(biāo)系統(tǒng)誤差改正，得到像點坐標(biāo)。（3）確定未知數(shù)的初始值。單像空間后方交會必須給出待定參數(shù)的初始值，在豎直航空攝影且地面控制點大體對稱分布的情況下，可按如下方法確定初始值：Z0s=H=m*f+Zi/4；X0s=xi/n;Y0s=Yi/n;t=0;式中：m為攝影比例尺分母；(4)計算旋轉(zhuǎn)矩陣R。利用角元素的近似值按下式計算方向余弦值a1,a2,a3,b1,b2,b3,c1,c2,c3，組成R陣。（5）逐點計算像點坐標(biāo)的近似值。利用未知數(shù)的近似值按共線條件方程計算控制點像點坐標(biāo)的近似值（x）、（y）。（6）按下式逐點計算誤差方程式的系數(shù)和常數(shù)

5、項，組成誤差方程。 (7)計算法方程的系數(shù)矩陣ATA與常數(shù)項ATL，組成法方程；（8）解求外方位元素。根據(jù)法方程，解求外方位元素的改正數(shù)，并與相應(yīng)的近似值求和，得到外方位元素新的近似值。（9）檢查計算是否收斂。將所求得的外方位元素的改正數(shù)與規(guī)定的限差比較，通常對t、Xs、Ys、Zs的改正數(shù)t，Xs，Ys，Zs給予限差，當(dāng)改正數(shù)小于限差時，迭代結(jié)束。否則用新的近似值重復(fù)（4）（8）步驟計算，直到滿足要求為止。（10）空間后方交會的精度估計：按上述方法所求得的影像外方位元素的精度可以通過法方程式中未知數(shù)的系數(shù)矩陣的逆陣（ATA）-1來解求，此時視像點坐標(biāo)為等精度不相關(guān)觀測值。因為ATA）-1中第i

6、個主對角線上的元素Qii就是法方程式中第i個未知數(shù)的權(quán)倒數(shù)，若單位權(quán)中誤差為m0，則第i個未知數(shù)的中誤差為：mi=0當(dāng)參加空間后方交會的控制點有n個時，則單位權(quán)中誤差可按下式計算：六 主要代碼與詳解void R(double t,double w,double k,double *a,double *b,double *c) a0=cos(t)*cos(k)-sin(t)*sin(w)*sin(k); a1=-cos(t)*sin(k)-sin(t)*sin(w)*cos(k); a2=-sin(t)*cos(w); b0=cos(w)*sin(k); b1=cos(w)*cos(k); b2

7、=-sin(w); c0=sin(t)*cos(k)+cos(t)*sin(w)*sin(k); c1=-sin(t)*sin(k)+cos(t)*sin(w)*cos(k); c2=cos(t)*cos(w);/子函數(shù)計算旋轉(zhuǎn)矩陣R。利用角元素的近似值按下式計算方向余弦值a1,a2,a3,b1,b2,b3,c1,c2,c3，組成R陣。void main() int i,m,num;double t,w,k,Xs,Ys,Zs,f; /六個外方位元素與焦距f double xN=-86.15,-53.40,-14.78,10.46,yN=-68.99,82.21,-76.63,64.43;/4個

8、像點坐標(biāo)double XN=36589.41,37631.08,39100.97,40426.54,YN=25273.32,31324.51,24934.98,30319.81,ZN=2195.17,728.69,2386.50,757.31;/四個控制點的空間坐標(biāo)m=50000;f=153.24;/ 影像比例尺1/m 與焦距f/ 以上主要是已知數(shù)據(jù)的定義double V6=0,a3,b3,c3,XoN,YoN,ZoN,A6*M,B6*M,lM,C36,D6,E8;for(i=0;iN;i+) xi=xi/1000.0; yi=yi/1000.0;/像點坐標(biāo)單位換算 / 以下為確定未知數(shù)的初始

9、值。 Xs=(X0+X1+X2+X3)/N; Ys=(Y0+Y1+Y2+Y3)/N; Zs=m*f+(Z0+Z1+Z2+Z3)/N;/ Xs，Ys，Zs 為攝站點的空間坐標(biāo)初始值 f=f/1000.0; t=w=k=0.0; /角元素的近似值 for(num=1;num10;num+) /num控制循環(huán)迭代次數(shù)R(t,w,k,a,b,c);/ 計算旋轉(zhuǎn)矩陣R。利用角元素的近似值按下式計算方向余弦值a1,a2,a3,b1,b2,b3,c1,c2,c3，組成R陣。 for(i=0;iN;i+)/下面的循環(huán)用來逐點計算誤差方程式的系數(shù)和常數(shù)項，組成誤差方程。 /以下是用共線條件方程計算控制點像點坐標(biāo)

10、的近似值（x）、（y）。 Xoi=-f*(a0*(Xi-Xs)+b0*(Yi-Ys)+c0*(Zi-Zs)/(a2*(Xi-Xs)+b2*(Yi-Ys)+c2*(Zi-Zs); Yoi=-f*(a1*(Xi-Xs)+b1*(Yi-Ys)+c1*(Zi-Zs)/(a2*(Xi-Xs)+b2*(Yi-Ys)+c2*(Zi-Zs);Zoi=a2*(Xi-Xs)+b2*(Yi-Ys)+c2*(Zi-Zs);/Zo便于后面計算 A12*i+0=(a0*f+a2*xi)/Zoi; A12*i+1=(b0*f+b2*xi)/Zoi; A12*i+2=(c0*f+c2*xi)/Zoi; A12*i+3=yi*

11、sin(w)-(xi*(xi*cos(k)-yi*sin(k)/f+f*cos(k)*cos(w); A12*i+4=-f*sin(k)-xi*(xi*sin(k)+yi*cos(k)/f; A12*i+5=yi; A12*i+6=(a1*f+a2*yi)/Zoi; A12*i+7=(b1*f+b2*yi)/Zoi; A12*i+8=(c1*f+c2*yi)/Zoi; A12*i+9=-xi*sin(w)-(yi*(xi*cos(k)-yi*sin(k)/f-f*sin(k)*cos(w); A12*i+10=-f*cos(k)-yi*(xi*sin(k)+yi*cos(k)/f; A12*i

12、+11=-xi; l2*i=xi-Xoi; l2*i+1=yi-Yoi;/計算l zhuanzhi(A,B,8,6); /轉(zhuǎn)置A xiangchen(B,A,C,6,8,6);/求ATA xiangchen(B,l,D,6,8,1); /求常數(shù)項ATL qiuni(C,6);/求ATA的逆 xiangchen(C,D,V,6,6,1);/求改正數(shù)V Xs+=V0; Ys+=V1; Zs+=V2; t+=V3; w+=V4; k+=V5; /結(jié)果改正 if(fabs(V0)=0.00001)&(fabs(V1)=0.00001)&(fabs(V2)=0.00001)&(fabs(V3)=0.00

13、001)&(fabs(V4)=0.00001)&(fabs(V5)=0.00001) break; /限差控制。檢查計算是否收斂。將所求得的外方位元素的改正數(shù)與規(guī)定的限差比較。當(dāng)改正數(shù)小于限差時，迭代結(jié)束。否則用新的近似值重復(fù)計算，直到滿足要求為止 。 /以下是計算單位權(quán)中誤差double s=0,m0,v8; xiangchen(A,V,E,8,6,1); /計算AX for(i=0;i8;i+) vi=Ei-li;/計算AX-L s+=vi*vi;/計算vv m0=sqrt(s/2);/ 單位權(quán)中誤差按下式計算：，這里n=4. R(t,w,k,a,b,c);/輸出結(jié)果 printf(像主點

14、的空間坐標(biāo)為：n); printf(Xs=%.2fn,Xs); printf(Ys=%.2fn,Ys); printf(Zs=%.2fn,Zs); printf(旋轉(zhuǎn)矩陣R為:n); for(i=0;i3;i+) printf( %.5f,ai); printf(n); for(i=0;i3;i+) printf( %.5f,bi); printf(n); for(i=0;i3;i+) printf( %.5f,ci); printf(n); printf(單位權(quán)中誤差為:,m0); coutm0endl;七 實習(xí)結(jié)果八 實習(xí)總結(jié)通過本次實習(xí)我深刻地理解單片空間后方交會的原理，尤其是共線方程。學(xué)會了在有多余觀測情況下，用最小二乘平差方法編程實現(xiàn)解求影像外方位元素