空間后方交會編程實習匯報一實習目的用程序設計語言（VisualC++或者C語言）編寫一種完整的單片空間后方交會程序，通過對提供的試驗數據進行計算，輸出像片的外方位元素并評估精度。本試驗的目的在于讓學生深入理解單片空間后方交會的原理，體會在有多出觀測狀況下，用最小二乘平差措施編程實現解求影像外方位元素的過程。通過上機調試程序加強動手能力的培養，通過對試驗成果的分析，增強學生綜合運用所學知識處理實際問題的能力。二實習內容運用一定數量的地面控制點，根據共線條件方程求解像片外方位元素。三實習數據已知航攝儀的內方位元素：fk=153.24mm，x0=y0=0.0mm，攝影比例尺為1:50000；4個地面控制點的地面坐標及其對應像點的像片坐標：點號像片坐標(mm)地面點坐標(m)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四實習原理假如我們懂得每幅影像的6個外方位元素，就能確定被攝物體與航攝影像的關系。因此，怎樣獲取影像的外方位元素，一直是攝影測量工作者所探討的問題。可采用的措施有：運用雷達、全球定位系統（GPS）、慣性導航系統（INS）以及星相攝影機來獲取影像的外方位元素；也可以運用影像覆蓋范圍內一定數量的控制點的空間坐標與攝影坐標，根據共線條件方程，反求該影像的外方位元素，這種措施稱為單幅影像的空間后方交會。單像空間后方交會的基本思想是：以單幅影像為基礎，從該影像所覆蓋地面范圍內若干控制點的已知地面坐標和對應點的像坐標量測值出發，根據共線條件方程，解求該影像在航空攝影時刻的外方位元素Xs，Ys，Zs，t，w，k。五實習流程（1）獲取已知數據。從攝影資料中查取影像比例尺1/m，平均攝影距離（航空攝影的航高、內方位元素x0，y0，f；獲取控制點的空間坐標Xt，Yt，Zt。（2）量測控制點的像點坐標并進行必要的影像坐標系統誤差改正，得到像點坐標。（3）確定未知數的初始值。單像空間后方交會必須給出待定參數的初始值，在豎直航空攝影且地面控制點大體對稱分布的狀況下，可按如下措施確定初始值：Z0s=H=m*f+ΣZi/4；X0s=Σxi/n;Y0s=ΣYi/n;t=ω=κ=0;式中：m為攝影比例尺分母；(4)計算旋轉矩陣R。運用角元素的近似值按下式計算方向余弦值a1,a2,a3,b1,b2,b3,c1,c2,c3，構成R陣。（5）逐點計算像點坐標的近似值。運用未知數的近似值按共線條件方程計算控制點像點坐標的近似值（x）、（y）。（6）按下式逐點計算誤差方程式的系數和常數項，構成誤差方程。(7)計算法方程的系數矩陣ATA與常數項ATL，構成法方程；（8）解求外方位元素。根據法方程，解求外方位元素的改正數，并與對應的近似值求和，得到外方位元素新的近似值。（9）檢查計算與否收斂。將所求得的外方位元素的改正數與規定的限差比較，一般對t、ω、κ、Xs、Ys、Zs的改正數Δt，Δω，Δκ，ΔXs，ΔYs，ΔZs予以限差，當改正數不不小于限差時，迭代結束。否則用新的近似值反復（4）——（8）環節計算，直到滿足規定為止。（10）空間后方交會的精度估計：按上述措施所求得的影像外方位元素的精度可以通過法方程式中未知數的系數矩陣的逆陣（ATA）-1來解求，此時視像點坐標為等精度不有關觀測值。由于ATA）-1中第i個主對角線上的元素Qii就是法方程式中第i個未知數的權倒數，若單位權中誤差為m0，則第i個未知數的中誤差為：mi=0當參與空間后方交會的控制點有n個時，則單位權中誤差可按下式計算：六重要代碼與詳解voidR(doublet,doublew,doublek,double*a,double*b,double*c){ a[0]=cos(t)*cos(k)-sin(t)*sin(w)*sin(k);a[1]=-cos(t)*sin(k)-sin(t)*sin(w)*cos(k);a[2]=-sin(t)*cos(w);b[0]=cos(w)*sin(k);b[1]=cos(w)*cos(k);b[2]=-sin(w);c[0]=sin(t)*cos(k)+cos(t)*sin(w)*sin(k);c[1]=-sin(t)*sin(k)+cos(t)*sin(w)*cos(k);c[2]=cos(t)*cos(w);}//子函數計算旋轉矩陣R。運用角元素的近似值按下式計算方向余弦值a1,a2,a3,b1,b2,b3,c1,c2,c3，構成R陣。voidmain(){inti,m,num;doublet,w,k,Xs,Ys,Zs,f;//六個外方位元素與焦距fdoublex[N]={-86.15,-53.40,-14.78,10.46},y[N]={-68.99,82.21,-76.63,64.43};//4個像點坐標doubleX[N]={36589.41,37631.08,39100.97,40426.54},Y[N]={25273.32,31324.51,24934.98,30319.81},Z[N]={2195.17,728.69,2386.50,757.31};//四個控制點的空間坐標m=50000;f=153.24;//影像比例尺1/m與焦距f//以上重要是已知數據的定義doubleV[6]={0},a[3],b[3],c[3],Xo[N],Yo[N],Zo[N],A[6*M],B[6*M],l[M],C[36],D[6],E[8];for(i=0;i<N;i++){x[i]=x[i]/1000.0;y[i]=y[i]/1000.0;//像點坐標單位換算}//如下為確定未知數的初始值。Xs=(X[0]+X[1]+X[2]+X[3])/N;Ys=(Y[0]+Y[1]+Y[2]+Y[3])/N;Zs=m*f+(Z[0]+Z[1]+Z[2]+Z[3])/N;//Xs，Ys，Zs為攝站點的空間坐標初始值f=f/1000.0;t=w=k=0.0;//角元素的近似值for(num=1;num<10;num++)//num控制循環迭代次數{ R(t,w,k,a,b,c);//計算旋轉矩陣R。運用角元素的近似值按下式計算方向余弦值a1,a2,a3,b1,b2,b3,c1,c2,c3，構成R陣。for(i=0;i<N;i++)//下面的循環用來逐點計算誤差方程式的系數和常數項，構成誤差方程。{//如下是用共線條件方程計算控制點像點坐標的近似值（x）、（y）。Xo[i]=-f*(a[0]*(X[i]-Xs)+b[0]*(Y[i]-Ys)+c[0]*(Z[i]-Zs))/(a[2]*(X[i]-Xs)+b[2]*(Y[i]-Ys)+c[2]*(Z[i]-Zs));Yo[i]=-f*(a[1]*(X[i]-Xs)+b[1]*(Y[i]-Ys)+c[1]*(Z[i]-Zs))/(a[2]*(X[i]-Xs)+b[2]*(Y[i]-Ys)+c[2]*(Z[i]-Zs));Zo[i]=a[2]*(X[i]-Xs)+b[2]*(Y[i]-Ys)+c[2]*(Z[i]-Zs);//Zo便于背面計算A[12*i+0]=(a[0]*f+a[2]*x[i])/Zo[i];A[12*i+1]=(b[0]*f+b[2]*x[i])/Zo[i];A[12*i+2]=(c[0]*f+c[2]*x[i])/Zo[i];A[12*i+3]=y[i]*sin(w)-(x[i]*(x[i]*cos(k)-y[i]*sin(k))/f+f*cos(k))*cos(w);A[12*i+4]=-f*sin(k)-x[i]*(x[i]*sin(k)+y[i]*cos(k))/f;A[12*i+5]=y[i];A[12*i+6]=(a[1]*f+a[2]*y[i])/Zo[i];A[12*i+7]=(b[1]*f+b[2]*y[i])/Zo[i];A[12*i+8]=(c[1]*f+c[2]*y[i])/Zo[i];A[12*i+9]=-x[i]*sin(w)-(y[i]*(x[i]*cos(k)-y[i]*sin(k))/f-f*sin(k))*cos(w);A[12*i+10]=-f*cos(k)-y[i]*(x[i]*sin(k)+y[i]*cos(k))/f;A[12*i+11]=-x[i];l[2*i]=x[i]-Xo[i];l[2*i+1]=y[i]-Yo[i];//計算l}zhuanzhi(A,B,8,6);//轉置Axiangchen(B,A,C,6,8,6);//求ATAxiangchen(B,l,D,6,8,1);//求常數項ATLqiuni(C,6);//求ATA的逆xiangchen(C,D,V,6,6,1);//求改正數VXs+=V[0];Ys+=V[1];Zs+=V[2];t+=V[3];w+=V[4];k+=V[5];//成果改正 if((fabs(V[0])<=0.00001)&&(fabs(V[1])<=0.00001)&&(fabs(V[2])<=0.00001)&&(fabs(V[3])<=0.00001)&&(fabs(V[4])<=0.00001)&&(fabs(V[5])<=0.00001))break;//限差控制。檢查計算與否收斂。將所求得的外方位元素的改正數與規定的限差比較。當改正數不不小于限差時，迭代結束。否則用新的近似值反復計算，直到滿足規定為止。}//如下是計算單位權中誤差doubles=0,m0,v[8];xiangchen(A,V,E,8,6,1);//計算AXfor(i=0;i<8;i++){v[i]=E[i]-l[i];//計算AX-L s+=v[i]*v[i];//計算[vv]}m0=sqrt(s/2);//單位權中誤差按下式計算：，這里n=4.R(t,w,k,a,b,c);//輸出成果printf("像主點的空間坐標為：\n");printf("Xs=%.2f\n",Xs);printf("Ys=%.2f\n",Ys);printf("Zs=%.2f\n",Zs);printf("旋轉矩陣R為:\n");for(i=0;i<3;i++) printf("%.5f",a[i]);printf("\n");for(i=0;i<3;i++) printf("%.5f",b[i]);printf("\n");for(i=0;i<3;i++) printf("%.5f",c[i]);printf("\n");printf("單位權中誤差為:",m0);cout<<m0<<endl;}七實習成果八實習總結通過本次實習我深刻地理解單片空間后方交會的原理，尤其是共線方