第一篇：攝影測量空間后前交會VC++

#include #include using namespace std;const int N=4;const int n=6;/*---------------矩陣相乘---------------------*/

void mult(double *m1,double *m2,double *result,int i_1,int j_12,int j_2){

int i,j,k;for(i=0;i

for(j=0;j

result[i*j_2+j]=0.0;

for(k=0;k

result[i*j_2+j]+=m1[i*j_12+k]*m2[j+k*j_2];

} return;} /*---------------矩陣求逆---------------------*/ void inverse(double c[n][n]){

int i,j,h,k;double p;double q[n][12];for(i=0;i

for(j=0;j

q[i][j]=c[i][j];

for(i=0;i

for(j=n;j<12;j++){

if(i+6==j)

q[i][j]=1;

else q[i][j]=0;

}

for(h=k=0;k

for(i=k+1;i

if(q[i][h]==0)

continue;

p=q[k][h]/q[i][h];

for(j=0;j<12;j++){

q[i][j]*=p;

q[i][j]-=q[k][j];

} } for(h=k=n-1;k>0;k--,h--)// 消去對角線以上的數(shù)據(jù)

for(i=k-1;i>=0;i--){

if(q[i][h]==0)

continue;

p=q[k][h]/q[i][h];

for(j=0;j<12;j++){

q[i][j]*=p;

q[i][j]-=q[k][j];

} } for(i=0;i

p=1.0/q[i][i];

for(j=0;j<12;j++)

q[i][j]*=p;} for(i=0;i

for(j=0;j

c[i][j]=q[i][j+6];} /*---------------矩陣轉(zhuǎn)置---------------------*/

void transpose(double *m1,double *m2,int m,int n){ //矩陣轉(zhuǎn)置

int i,j;

for(i=0;i

for(j=0;j

m2[j*m+i]=m1[i*n+j];

return;

} void main(){ double Xs,Ys,Zs,q,w,k;double Xsr,Ysr,Zsr,qr,wr,kr;// double a[3],b[3],c[3];double ar[3],br[3],cr[3];// double x0,y0,f;double x[N],y[N];double X[N],Y[N],Z[N];double x1[N],y1[N];double xr[N],yr[N];// double x1r[N],y1r[N];// double m;double L[2*N];double Lr[2*N];// double XX[6];double XXr[6];// double A[2*N][6];double B[2*N][6];// double X0[N],Y0[N],Z0[N],At[6][2*N],result1[6][6],result2[6][1];double X0r[N],Y0r[N],Z0r[N],Bt[6][2*N],result1r[6][6],result2r[6][1];// double R1[3][3],R2[3][3],XI1[3],XII1[3],XI2[3],XII2[3],XI3[3],XII3[3],XI4[3],XII4[3],XI5[3],XII5[3];double N1,N2,N3,N4,N5,BX,BY,BZ,G[5][3];double RS1[3][1],RSR1[3][1],RS2[3][1],RSR2[3][1],RS3[3][1],RSR3[3][1],RS4[3][1],RSR4[3][1],RS5[3][1],RSR5[3][1];int i,j,n=0,nr=0;double sum=0;double m1,n1,m2,n2,m3,n3,m4,n4,m5,n5,p1,q1,p2,q2,p3,q3,p4,q4,p5,q5;

/*---------------輸入點(diǎn)地面坐標(biāo)---------------------*/ X[0]=5083.205;X[1]=5780.02;X[2]=5210.879;X[3]=5909.264;Y[0]=5852.099;Y[1]=5906.365;Y[2]=4258.446;Y[3]=4314.283;Z[0]=527.925;Z[1]=571.549;Z[2]=461.81;Z[3]=455.48;m1=0.051758;n1=0.081555;p1=-0.039953;q1=0.078463;m2=0.014618;n2=-0.000231;p2=-0.076016;q2=0.000036;m3=0.04988;n3=-0.000792;p3=-0.042201;q3=-0.001022;m4=0.086243;n4=-0.001346;p4=-0.007706;q4=-0.002112;m5=0.048135;n5=-0.079962;p5=-0.044438;q5=-0.079736;/*---------------輸入點(diǎn)像片坐標(biāo)---------------------*/ x[0]=16.012;x[1]=88.56;x[2]=13.362;x[3]=82.24;y[0]=79.963;y[1]=81.134;y[2]=-79.37;y[3]=-80.027;

xr[0]=-73.93;xr[1]=-5.252;xr[2]=-79.122;xr[3]=-9.887;yr[0]=78.706;yr[1]=78.184;yr[2]=-78.879;yr[3]=-80.089;

/*-----------------設(shè)定外方位元素初始值--------------*/ x0=0;y0=0;f=152.00;m=10000;Xs=0;Ys=0;Zs=f*m/1000;Xsr=0;Ysr=0;Zsr=f*m/1000;// q=0;w=0;k=0;qr=0;wr=0;kr=0;// XX[3]=1;XXr[3]=1;// /*------------------迭代計(jì)算左片Xs,Ys,Zs,q,w,k--------------------------*/ while((XX[3]>6/206265 || XX[4]>6/206265 || XX[5]>6/206265)&&n<100){// /*----------------旋轉(zhuǎn)矩陣R-----------------------*/ a[0]=cos(q)*cos(k)-sin(q)*sin(w)*sin(k);a[1]=-cos(q)*sin(k)-sin(q)*sin(w)*cos(k);a[2]=-sin(q)*cos(w);

b[0]=cos(w)*sin(k);b[1]=cos(w)*cos(k);b[2]=-sin(w);c[0]=sin(q)*cos(k)+cos(q)*sin(w)*sin(k);c[1]=-sin(q)*sin(k)+cos(q)*sin(w)*cos(k);c[2]=cos(q)*cos(w);/*-----------------像點(diǎn)坐標(biāo)計(jì)算值------------------*/ for(i=0;i

A[2*i][0]=((a[0]*f+a[2]*(x[i]-x0)))/Z0[i];

A[2*i][1]=((b[0]*f+b[2]*(x[i]-x0)))/Z0[i];

A[2*i][2]=((c[0]*f+c[2]*(x[i]-x0)))/Z0[i];

A[2*i][3]=(y[i]-y0)*sin(w)-((x[i]-x0)*((x[i]-x0)*cos(k)-y[i]*sin(k))/f+f*cos(k))*cos(w);

A[2*i][4]=-f*sin(k)-(x[i]-x0)*((x[i]-x0)*sin(k)+(y[i]-y0)*cos(k))/f;

A[2*i][5]=y[i]-y0;

L[2*i]=x[i]-x1[i];

A[1+2*i][0]=((a[1]*f+a[2]*(y[i]-y0)))/Z0[i];

A[1+2*i][1]=((b[1]*f+b[2]*(y[i]-y0)))/Z0[i];

A[1+2*i][2]=((c[1]*f+c[2]*(y[i]-y0)))/Z0[i];

A[1+2*i][3]=-(x[i]-x0)*sin(w)-((y[i]-y0)*((x[i]-x0)*cos(k)-(y[i]-y0)*sin(k))/f-f*sin(k))*cos(w);

A[1+2*i][4]=-f*cos(k)-(y[i]-y0)*((x[i]-x0)*sin(k)+(y[i]-y0)*cos(k))/f;

A[1+2*i][5]=-x[i]+x0;

L[1+2*i]=y[i]-y1[i];}

/*-------------------解法方程--------------------*/ transpose(&A[0][0],&At[0][0],2*N,6);//求A轉(zhuǎn)置矩陣

mult(&At[0][0],&A[0][0],&result1[0][0],6,2*N,6);//A矩陣 X A轉(zhuǎn)置矩陣-> result1 inverse(result1);// result1 求逆

mult(&At[0][0],L,&result2[0][0],6,2*N,1);// A轉(zhuǎn)置矩陣 X L-> result2 mult(&result1[0][0],&result2[0][0],&XX[0],6,6,1);// result1 X result2-> XX Xs+=XX[0];

} Ys+=XX[1];Zs+=XX[2];q+=XX[3];w+=XX[4];k+=XX[5];n++;/*------------------迭代計(jì)算有片Xsr,Ysr,Zsr,qr,wr,kr--------------------------*/

while((XXr[3]>6/206265 || XXr[4]>6/206265 || XXr[5]>6/206265)&&nr<100){//右點(diǎn)算Xs,Ys,Zs,q,w,k /*----------------旋轉(zhuǎn)矩陣R-----------------------*/

ar[0]=cos(qr)*cos(kr)-sin(qr)*sin(wr)*sin(kr);

ar[1]=-cos(qr)*sin(kr)-sin(qr)*sin(wr)*cos(kr);

ar[2]=-sin(qr)*cos(wr);

br[0]=cos(wr)*sin(kr);

br[1]=cos(wr)*cos(kr);

br[2]=-sin(wr);

cr[0]=sin(qr)*cos(kr)+cos(qr)*sin(wr)*sin(kr);

cr[1]=-sin(qr)*sin(kr)+cos(qr)*sin(wr)*cos(kr);

cr[2]=cos(qr)*cos(wr);

/*-----------------像點(diǎn)坐標(biāo)計(jì)算值------------------*/

for(i=0;i

X0r[i]=ar[0]*(X[i]-Xsr)+br[0]*(Y[i]-Ysr)+cr[0]*(Z[i]-Zsr);

Y0r[i]=ar[1]*(X[i]-Xsr)+br[1]*(Y[i]-Ysr)+cr[1]*(Z[i]-Zsr);

Z0r[i]=ar[2]*(X[i]-Xsr)+br[2]*(Y[i]-Ysr)+cr[2]*(Z[i]-Zsr);

x1r[i]=x0-f*X0r[i]/Z0r[i];

y1r[i]=y0-f*Y0r[i]/Z0r[i];

}

/*-------------誤差方程中各偏導(dǎo)數(shù)的值--------------*/

for(i=0;i

B[2*i][0]=((ar[0]*f+ar[2]*(xr[i]-x0)))/Z0r[i];

B[2*i][1]=((br[0]*f+br[2]*(xr[i]-x0)))/Z0r[i];

B[2*i][2]=((cr[0]*f+cr[2]*(xr[i]-x0)))/Z0r[i];

B[2*i][3]=(yr[i]-y0)*sin(wr)-((xr[i]-x0)*((xr[i]-x0)*cos(kr)-yr[i]*sin(kr))/f+f*cos(kr))*cos(wr);

B[2*i][4]=-f*sin(kr)-(xr[i]-x0)*((xr[i]-x0)*sin(kr)+(yr[i]-y0)*cos(kr))/f;

B[2*i][5]=yr[i]-y0;

Lr[2*i]=xr[i]-x1r[i];

B[1+2*i][0]=((ar[1]*f+ar[2]*(yr[i]-y0)))/Z0r[i];

B[1+2*i][1]=((br[1]*f+br[2]*(yr[i]-y0)))/Z0r[i];

B[1+2*i][2]=((cr[1]*f+cr[2]*(yr[i]-y0)))/Z0r[i];

B[1+2*i][3]=-(xr[i]-x0)*sin(wr)-((yr[i]-y0)*((xr[i]-x0)*cos(kr)-(yr[i]-y0)*sin(kr))/f-f*sin(kr))*cos(wr);

}

} B[1+2*i][4]=-f*cos(kr)-(yr[i]-y0)*((xr[i]-x0)*sin(kr)+(yr[i]-y0)*cos(kr))/f;B[1+2*i][5]=-xr[i]+x0;Lr[1+2*i]=yr[i]-y1r[i];/*-------------------解法方程--------------------*/ transpose(&B[0][0],&Bt[0][0],2*N,6);//求A轉(zhuǎn)置矩陣

mult(&Bt[0][0],&B[0][0],&result1r[0][0],6,2*N,6);//A矩陣 X A轉(zhuǎn)置矩陣-> result1 inverse(result1r);// result1 求逆

mult(&Bt[0][0],Lr,&result2r[0][0],6,2*N,1);// A轉(zhuǎn)置矩陣 X L-> result2 mult(&result1r[0][0],&result2r[0][0],&XXr[0],6,6,1);// result1 X result2-> XX Xsr+=XXr[0];Ysr+=XXr[1];Zsr+=XXr[2];qr+=XXr[3];wr+=XXr[4];kr+=XXr[5];nr++;BX=Xsr-Xs;//BX,BY,BZ BY=Ysr-Ys;BZ=Zsr-Zs;for(i=0;i<3;i++){//R1,R2賦值

R1[0][i]=a[i];

R1[1][i]=b[i];

R1[2][i]=c[i];

R2[0][i]=ar[i];

R2[1][i]=br[i];

R2[2][i]=cr[i];} RS1[0][0]=m1;RS1[1][0]=n1;RS1[2][0]=-f/1000;// 給 將要與R1，R2相乘的 RS矩陣和RSR矩陣 賦值

RS2[0][0]=m2;RS2[1][0]=n2;RS2[2][0]=-f/1000;RS3[0][0]=m3;RS3[1][0]=n3;RS3[2][0]=-f/1000;RS4[0][0]=m4;RS4[1][0]=n4;RS4[2][0]=-f/1000;RS5[0][0]=m5;RS5[1][0]=n5;RS5[2][0]=-f/1000;

RSR1[0][0]=p1;RSR1[1][0]=q1;RSR1[2][0]=-f/1000;RSR2[0][0]=p2;RSR2[1][0]=q2;RSR2[2][0]=-f/1000;RSR3[0][0]=p3;RSR3[1][0]=q3;RSR3[2][0]=-f/1000;RSR4[0][0]=p4;RSR4[1][0]=q4;RSR4[2][0]=-f/1000;RSR5[0][0]=p5;RSR5[1][0]=q5;RSR5[2][0]=-f/1000;

mult(&R1[0][0],&RS1[0][0],&XI1[0],3,3,1);//R1 X RS mult(&R1[0][0],&RS2[0][0],&XI2[0],3,3,1);mult(&R1[0][0],&RS3[0][0],&XI3[0],3,3,1);mult(&R1[0][0],&RS4[0][0],&XI4[0],3,3,1);mult(&R1[0][0],&RS5[0][0],&XI5[0],3,3,1);mult(&R2[0][0],&RSR1[0][0],&XII1[0],3,3,1);//R2 X RSR mult(&R2[0][0],&RSR2[0][0],&XII2[0],3,3,1);mult(&R2[0][0],&RSR3[0][0],&XII3[0],3,3,1);mult(&R2[0][0],&RSR4[0][0],&XII4[0],3,3,1);mult(&R2[0][0],&RSR5[0][0],&XII5[0],3,3,1);N1=(BX*XII1[2]-BZ*XII1[0])/(XI1[0]*XII1[2]-XI1[2]*XII1[0]);//給N1 賦值 有五個(gè)點(diǎn) N2=(BX*XII2[2]-BZ*XII2[0])/(XI2[0]*XII2[2]-XI2[2]*XII2[0]);N3=(BX*XII3[2]-BZ*XII3[0])/(XI3[0]*XII3[2]-XI3[2]*XII3[0]);N4=(BX*XII4[2]-BZ*XII4[0])/(XI4[0]*XII4[2]-XI4[2]*XII4[0]);N5=(BX*XII5[2]-BZ*XII5[0])/(XI5[0]*XII5[2]-XI5[2]*XII5[0]);G[0][0]=Xs+N1*XI1[0];//五個(gè)點(diǎn)XA,YA,ZA計(jì)算G[i][j] i:第幾個(gè)點(diǎn)，j：是x，y還是z G[0][1]=Ys+N1*XI1[1];G[0][2]=Zs+N1*XI1[2];G[1][0]=Xs+N2*XI1[0];G[1][1]=Ys+N2*XI1[1];G[1][2]=Zs+N2*XI1[2];G[2][0]=Xs+N3*XI1[0];G[2][1]=Ys+N3*XI1[1];G[2][2]=Zs+N3*XI1[2];G[3][0]=Xs+N4*XI1[0];G[3][1]=Ys+N4*XI1[1];G[3][2]=Zs+N4*XI1[2];G[4][0]=Xs+N5*XI1[0];G[4][1]=Ys+N5*XI1[1];G[4][2]=Zs+N5*XI1[2];/*----------------旋轉(zhuǎn)矩陣R-----------------------*/ cout<<“--------左--------”<

”<

cout<<“Ys

”<

cout<<“Zs

”<

cout<<“迭代次數(shù)為:”<

”<

cout<<“Ysr

”<

cout<<“Zsr

”<

5”<

%-.5f

%.5f

%-.5f

%.5f

%-.5f n”,XI1[0],XI2[0],XI3[0],XI4[0],XI5[0]);printf(“ Y1

%-.5f

%.5f

%-.5f

%.5f

%-.5f n”,XI1[1],XI2[1],XI3[1],XI4[1],XI5[1]);printf(“ Z1

%-.5f

%.5f

%-.5f

%.5f

%-.5f n”,XI1[2],XI2[2],XI3[2],XI4[2],XI5[2]);printf(“ X2

%-.5f

%.5f

%-.5f

%.5f

%-.5f n”,XII1[0],XII2[0],XII3[0],XII4[0],XII5[0]);printf(“ Y2

%-.5f

%.5f

%-.5f

%.5f

%-.5f n”,XII1[1],XII2[1],XII3[1],XII4[1],XII5[1]);printf(“ Z2

%-.5f

%.5f

%-.5f

%.5f

%-.5f n”,XII1[2],XII2[2],XII3[2],XII4[2],XII5[2]);

printf(“ N1

%-.5f

%.5f

%-.5f

%.5f

%-.5f nn”,N1,N2,N3,N4,N5);printf(“ XA

%-.5f

%.5f

%-.5f

%.5f

%-.5f n”,G[0][0],G[1][0],G[2][0],G[3][0],G[4][0]);printf(“ YA

%-.5f

%.5f

%-.5f

%.5f

%-.5f n”,G[0][1],G[1][1],G[2][1],G[3][1],G[4][1]);printf(“ ZA

%-.5f

%.5f

%-.5f

%.5f

%-.5f n”,G[0][2],G[1][2],G[2][2],G[3][2],G[4][2]);

}

第二篇：攝影測量實(shí)習(xí)報(bào)告-單片空間后方交會

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

實(shí)習(xí)內(nèi)容:單片空間后方交會編程

實(shí)習(xí)者:李友兵

學(xué)號:0810050121

指導(dǎo)老師:張金平老師

實(shí)習(xí)時(shí)間:2011.05.30——2011.06.03

一、實(shí)習(xí)目的與任務(wù)

此次攝影測量實(shí)習(xí)主要是要自主編程實(shí)現(xiàn)單像空間后方交會，通過已知的內(nèi)方位元素和控制點(diǎn)像點(diǎn)坐標(biāo)和地面坐標(biāo)求解六個(gè)外方位元素，在此過程中深入理解單像空間后方交會的原理和對編程的熟悉和理解（我用的是C語言編程），和對時(shí)間的合理運(yùn)用，對知識的綜合運(yùn)用，培養(yǎng)理論的實(shí)際運(yùn)用能力，任務(wù)是在一個(gè)星期內(nèi)自主完成。

二、單片空間后方交會理論基礎(chǔ)

單像空間后方交會：是通過以像點(diǎn)平面坐標(biāo)為觀測值，以控制點(diǎn)坐標(biāo)為已知值，利用共線條件方程和最小二乘原理，運(yùn)用間接平差方法，通過迭代求解6個(gè)外方位元素。程序設(shè)計(jì)的思路與流程

三、程序設(shè)計(jì)的思路與流程（編程框架和步驟）

1、根據(jù)內(nèi)方位元素和已知的數(shù)據(jù)將控制點(diǎn)的框標(biāo)坐標(biāo)轉(zhuǎn)換為像平面坐標(biāo)系坐標(biāo)：x=x′-x0，y=y′-y0

2、確定未知數(shù)的初始值：角元素初始值κ0=ω0=Φ0=0；線元素初始值Zs0=H=mf，Xs0=(X1+X2+X3+X4)/4,Ys0=（Y1+Y2+Y3+Y4）/4

3、利用角元素初始值計(jì)算方向余弦值組成旋轉(zhuǎn)矩陣R a1=cosΦ*cosκ-sinΦ*sinω*sinκ,a2=-cosΦ*sinκ-sinΦ*sinω*cosκ，a3=-sinΦ*cosω

b1=cosω*sinκ,b2=cosω*cosκ,b3=-sinω

c1=sinΦ*cosκ+cosΦ*sinω*sinκ,c2=-sinΦ*sinκ+cosΦ*sinω*cosκ,c3=cosΦ*cosω

4、逐點(diǎn)計(jì)算控制點(diǎn)的像點(diǎn)坐標(biāo)的近似值，共線方程為: x=(-f*(a1*(X-Xs)+b1*(Y-Ys)+c1*(Z-Zs)))/(a3*(X-Xs)+b3*(Y-Ys)+c3*(Z-Zs))

y=(-f*(a2*(X-Xs)+b2*(Y-Ys)+c2*(Z-Zs)))/(a3*(X-Xs)+b3*(Y-Ys)+c3*(Z-Zs))

(x)=(-f*(a1*(X-Xs)+b1*(Y-Ys)+c1*(Z-Zs)))/(a3*(X-Xs)+b3*(Y-Ys)+c3*(Z-Zs))(y)=(-f*(a2*(X-Xs)+b2*(Y-Ys)+c2*(Z-Zs)))/(a3*(X-Xs)+b3*(Y-Ys)+c3*(Z-Zs))

5、組成誤差方程：Vx=a11*dXs+a12*dYs+a13*dZs+a14*dΦ+a15*dω+a16*dκ+(x)-x Vy=a21*dXs+a22*dYs+a23*dZs+a24*dΦ+a25*dω+a26*dκ+(y)-y 系數(shù)求解(是共線方程分別外方位元素求導(dǎo)，是共線方程線性化的系數(shù))：變量代換

A= a1*(X-Xs)+b1*(Y-Ys)+c1*(Z-Zs)B= a2*(X-Xs)+b2*(Y-Ys)+c2*(Z-Zs)C= a3*(X-Xs)+b3*(Y-Ys)+c3*(Z-Zs)a11=(a1*f+a3*x)/C,a12=(b1*f+b3*x)/C,a13=(c1*f+c3*f)/C,a14=y*sinω-((x*(x*cosκ-y*sinκ))/f+f*cosκ)cosω,a15=-f*sinκ-(x*(x*sinκ+y*cosκ)/f),a16=y a21=(a2*f+a3*x)/C,a22=(b2*f+b3*x)/C,a23=（c2*f+c3*x）/C,a24=-xsinω-((x*(x*cosκ-y*sinκ))/f-f*sinκ)cosω,a25=-f*cosκ-(y*(x*sinκ+y*cosκ)/f),a26 =-x 誤差方程的常系數(shù)（是像點(diǎn)坐標(biāo)觀測值與計(jì)算的近似值的差值）： lx=x-(x),ly=y-(y)

6、組成法方程，解求外方位元素改正數(shù)X=(ATA)-1ATL(A為誤差方程的系數(shù)矩陣，L為誤差方程的常系數(shù)矩陣，通過步驟5求得，此處先求ATA再求矩陣的逆矩陣，解得的改正數(shù)加上相應(yīng)的近似值得到外方位元素新的近似值)

7、檢查計(jì)算是否收斂：將求得的外方位元素改正數(shù)與規(guī)定的限差比較，大于限差繼續(xù)迭代，小于限差則終止。

四、各子函數(shù)詳細(xì)設(shè)計(jì)的關(guān)鍵技術(shù)參數(shù)

子函數(shù)（輸入函數(shù)、Input，矩陣求積Matrixmultiply，計(jì)算函數(shù)Resection，矩陣轉(zhuǎn)置Matrixtranspose，矩陣求逆Matrixinverse，輸出函數(shù)Output）主要用了Matrixtranspose，矩陣的行變列，列變行,參數(shù)為要轉(zhuǎn)置的矩陣，轉(zhuǎn)置后的矩陣，要轉(zhuǎn)置矩陣的行列數(shù)，Matrixinverse,求矩陣的代數(shù)余子式，參數(shù)有要求逆的矩陣和，逆矩陣的行數(shù)，Matrixmultiply，一矩陣的行乘以二矩陣的列，參數(shù)為一矩陣，二矩陣，所求矩陣，一的行，一的列，二的列

五、像片外方位元素解算結(jié)

六、實(shí)習(xí)體會

從周一開始進(jìn)行整個(gè)實(shí)習(xí)框架進(jìn)行構(gòu)建，編寫了實(shí)習(xí)的思路既步驟；然后是梳理框架，要用到的子函數(shù)，子函數(shù)參數(shù)，怎么編寫，有不懂的東西網(wǎng)上參考資料；再是自己動(dòng)手編程，發(fā)現(xiàn)編程要求很細(xì)膩，出不得一點(diǎn)差錯(cuò)，程序編寫出來后，有很多的錯(cuò)誤，要求逐一進(jìn)行改正，修改；最后才得以運(yùn)行。通過本次實(shí)習(xí)，我深刻理解了單片空間后方交會原理，進(jìn)一步熟悉理解了C語言編程，認(rèn)識到了時(shí)間的搭配的重要性和參考資料的必要性，當(dāng)遇到難題是要迎難而上，達(dá)到突破，當(dāng)完成是能感到一絲絲欣慰和成就感。

附錄:源代碼

#include “stdio.h” #include “math.h” #include “Matrixmultiply.c” #include “Matrixtranspose.c” #include “Matrixinverse.c” void main(){ int i,j,k,f=0;double x0=0.00018, y0=0.00026,fk=0.15324;

//內(nèi)方位元素

double m=40000;//估算比例尺

double R[3][3],XG[6][1],AT[6][8],ATA[6][6],ATL[6][1];double Xs=0.0, Ys=0.0, Zs=0.0,Q=0.0,W=0.0,K=0.0;

double X,Y,Z,L[8][1],A[8][6];

double B[4][5]={-0.08615,-0.06899,36589.41,25273.32,2195.17,-0.05340,0.08221,37631.08,31324.51,728.69,-0.01478,-0.07663,39100.97,24934.98,2386.80,0.01046,0.06443,40426.54,30319.81,757.31};

for(i=0;i<4;i++){

Xs=Xs+B[i][2];

Ys=Ys+B[i][3];

Zs=Zs+B[i][4];} Xs=Xs/4;Ys=Ys/4;Zs=m*fk;//求得外方位線元素的初始值

do//迭代計(jì)算

{

f++;//迭代次數(shù)

//組成旋轉(zhuǎn)矩陣

R[0][0]=cos(Q)*cos(K)-sin(Q)*sin(W)*sin(K);

R[0][1]=-cos(Q)*sin(K)-sin(Q)*sin(W)*cos(K);

R[0][2]=-sin(Q)*cos(W);

R[1][0]=cos(W)*sin(K);

R[1][1]=cos(W)*cos(K);

R[1][2]=-sin(W);

R[2][0]=sin(Q)*cos(K)+cos(Q)*sin(W)*sin(K);

R[2][1]=-sin(Q)*sin(K)+cos(Q)*sin(W)*cos(K);

R[2][2]=cos(Q)*cos(W);

//計(jì)算系數(shù)陣和常數(shù)項(xiàng)

for(i=0, k=0,j=0;i<=3;i++,k++,j++)

{

X=R[0][0]*(B[i][2]-Xs)+R[1][0]*(B[i][3]-Ys)+R[2][0]*(B[i][4]-Zs);

Y=R[0][1]*(B[i][2]-Xs)+R[1][1]*(B[i][3]-Ys)+R[2][1]*(B[i][4]-Zs);

Z=R[0][2]*(B[i][2]-Xs)+R[1][2]*(B[i][3]-Ys)+R[2][2]*(B[i][4]-Zs);//為了計(jì)算簡單而進(jìn)行變量代換

L[j][0]=B[i][0]-(x0-fk*X/Z);

L[j+1][0]=B[i][1]-(y0-fk*Y/Z);

j++;

A[k][0]=(R[0][0]*fk+R[0][2]*(B[i][0]-x0))/Z;

A[k][1]=(R[1][0]*fk+R[1][2]*(B[i][0]-x0))/Z;

A[k][2]=(R[2][0]*fk+R[2][2]*(B[i][0]-x0))/Z;

A[k][3]=(B[i][1]-y0)*sin(W)-((B[i][0]-x0)*((B[i][0]-x0)*cos(K)-(B[i][1]-y0)*sin(K))/fk+fk*cos(K))*cos(W);

A[k][4]=-fk*sin(K)-(B[i][0]-x0)*((B[i][0]-x0)*sin(K)+(B[i][1]-y0)*cos(K))/fk;

A[k][5]=B[i][1]-y0;

A[k+1][0]=(R[0][1]*fk+R[0][2]*(B[i][1]-y0))/Z;

A[k+1][1]=(R[1][1]*fk+R[1][2]*(B[i][1]-y0))/Z;

A[k+1][2]=(R[2][1]*fk+R[2][2]*(B[i][1]-y0))/Z;

A[k+1][3]=-(B[i][0]-x0)*sin(W)-((B[i][1]-y0)*((B[i][0]-x0)*cos(K)-(B[i][1]-y0)*sin(K))/fk-fk*sin(K))*cos(W);

A[k+1][4]=-fk*cos(K)-(B[i][1]-y0)*((B[i][0]-x0)*sin(K)+(B[i][1]-y0)*cos(K))/fk;

A[k+1][5]=-(B[i][0]-x0);

k++;

}

Matrixtranspose(A,AT,8,6);//此為轉(zhuǎn)置函數(shù)的調(diào)用,求AT

Matrixmultiply(AT,A,ATA,6,8,6);//此為矩陣相乘函數(shù)的調(diào)用，求ATA

Matrixinverse(ATA,6);//此為求矩陣逆函數(shù)的調(diào)用，求ATA的逆

Matrixmultiply(AT,L,ATL,6,8,1);//此為矩陣相乘函數(shù)的調(diào)用，求ATL

Matrixmultiply(ATA,ATL,XG,6,6,1);//此為矩陣相乘函數(shù)的調(diào)用，求外方位元素改正數(shù)

Xs=Xs+XG[0][0];Ys=Ys+XG[1][0];Zs=Zs+XG[2][0];

Q=Q+XG[3][0];W=W+XG[4][0];K=K+XG[5][0];//初始值加外方位元素改正數(shù)進(jìn)行迭代

}while(XG[3][0]>=0.00000001||XG[4][0]>=0.00000001||XG[5][0]>=0.00000001);//當(dāng)限差滿足要求時(shí)要再一次進(jìn)行旋轉(zhuǎn)矩陣的求解

R[0][0]=cos(Q)*cos(K)-sin(Q)*sin(W)*sin(K);

R[0][1]=-cos(Q)*sin(K)-sin(Q)*sin(W)*cos(K);

R[0][2]=-sin(Q)*cos(W);

R[1][0]=cos(W)*sin(K);

R[1][1]=cos(W)*cos(K);

R[1][2]=-sin(W);

R[2][0]=sin(Q)*cos(K)+cos(Q)*sin(W)*sin(K);

R[2][1]=-sin(Q)*sin(K)+cos(Q)*sin(W)*cos(K);

R[2][2]=cos(Q)*cos(W);printf(“迭代次數(shù):%d”,f);

//屏幕輸出誤差方程系數(shù)陣、常數(shù)項(xiàng)、改正數(shù)

printf(“nn誤差方程系數(shù)矩陣A為:nn”);for(i=0;i<6;i++){

for(j=0;j<6;j++)

printf(“%13.5e ”,A[i][j]);

printf(“n”);}

printf(“n常數(shù)項(xiàng)L為:nn”);

for(i=0;i<8;i++){

for(j=0;j<1;j++)

printf(“%13.5e ”,L[i][j]);

printf(“n”);} printf(“n改正數(shù)XG為:nn”);

for(i=0;i<6;i++){

for(j=0;j<1;j++)

printf(“%13.5e ”,XG[i][j]);

printf(“n”);}

printf(“n相片的外方位元素為:nn”);

printf(“ Xs=%13.7e, Ys=%13.7e, Zs=%13.7e nn”,Xs,Ys,Zs);

printf(“ Q=%13.7e, W=%13.7e, K=%13.7e n”,Q,W,K);printf(“n旋轉(zhuǎn)矩陣R為:nn”);

for(i=0;i<3;i++){

for(j=0;j<3;j++)

printf(“%13.5e ”,R[i][j]);

printf(“n”);}

} //子函數(shù)

#include

#include

#include

int Matrixinverse(a,n)

int n;

double a[];

{ int *is,*js,i,j,k,l,u,v;

double d,p;

is=malloc(n*sizeof(int));

js=malloc(n*sizeof(int));

for(k=0;k<=n-1;k++)

{ d=0.0;

for(i=k;i<=n-1;i++)

for(j=k;j<=n-1;j++)

{ l=i*n+j;p=fabs(a[l]);

if(p>d){ d=p;is[k]=i;js[k]=j;}

}

if(d+1.0==1.0)

{ free(is);free(js);printf(“err**not invn”);

return(0);

}

if(is[k]!=k)

for(j=0;j<=n-1;j++)

{ u=k*n+j;v=is[k]*n+j;

p=a[u];a[u]=a[v];a[v]=p;

}

if(js[k]!=k)

for(i=0;i<=n-1;i++)

{ u=i*n+k;v=i*n+js[k];

p=a[u];a[u]=a[v];a[v]=p;

}

l=k*n+k;

a[l]=1.0/a[l];

for(j=0;j<=n-1;j++)

if(j!=k)

{ u=k*n+j;a[u]=a[u]*a[l];}

for(i=0;i<=n-1;i++)

if(i!=k)

for(j=0;j<=n-1;j++)

if(j!=k)

{ u=i*n+j;

a[u]=a[u]-a[i*n+k]*a[k*n+j];

}

for(i=0;i<=n-1;i++)

if(i!=k)

{ u=i*n+k;a[u]=-a[u]*a[l];}

}

for(k=n-1;k>=0;k--)

{ if(js[k]!=k)

for(j=0;j<=n-1;j++)

{ u=k*n+j;v=js[k]*n+j;

p=a[u];a[u]=a[v];a[v]=p;

}

if(is[k]!=k)

for(i=0;i<=n-1;i++)

{ u=i*n+k;v=i*n+is[k];

p=a[u];a[u]=a[v];a[v]=p;

}

}

free(is);free(js);

return(1);

} //子函數(shù)

void Matrixmultiply(a,b,c,m,n,k)int m,n,k;double a[],b[],c[];{ int i,j,l,u;for(i=0;i

for(j=0;j

{

u=i*k+j;c[u]=0.0;

for(l=0;l

c[u]+=a[i*n+l]*b[l*k+j];

}

return;} //子函數(shù)

void Matrixtranspose(a,b,m,n)int m,n;double a[],b[];{ int i,j,u;for(i=0;i

for(j=0;j

{

u=j*m+i;b[u]=0.0;b[j*m+i]=a[i*n+j];

} } return;}

第三篇：攝影測量

4D指的是DEM、DOM、DLG、DRG。意義如下：

數(shù)字高程模型（Digital Elevation Model 簡稱DEM）是在高斯投影平面上規(guī)則格網(wǎng)點(diǎn)平面坐標(biāo)（x,y）及其高程（z）的數(shù)據(jù)集。Dem的水平間隔可隨地貌類型不同而改變。根據(jù)不同的高程精度，可分為不同等級產(chǎn)品。

數(shù)字正射影像圖（Digital Orthophoto Map簡稱DOM）是利用數(shù)字高程模型對掃描處理的數(shù)字化的航空相片 / 遙感相片（單色 / 彩色），經(jīng)逐象元進(jìn)行糾正，再按影像鑲嵌，根據(jù)圖幅范圍剪裁生成的影像數(shù)據(jù)。一般帶有公里格網(wǎng)、圖廓內(nèi) / 外整飾和注記的平面圖。

數(shù)字線劃地圖（Digital Line Graphic簡稱DLG）是現(xiàn)有地形圖上基礎(chǔ)地理要素的矢量數(shù)據(jù)集，且保存要素間空間關(guān)系和相關(guān)的屬性信息。

數(shù)字柵格地圖（Digital Raster Graphic簡稱DRG）是紙質(zhì)地形圖的數(shù)字化產(chǎn)品。每幅圖經(jīng)掃描、糾正、圖幅處理及數(shù)據(jù)壓縮處理后，形成在內(nèi)容、幾何精度和色彩上與地形圖保持一致的柵格文件。

第四篇：攝影測量

像點(diǎn)位移：當(dāng)像片傾斜、地面起伏時(shí)，地面點(diǎn)在航攝像片上構(gòu)像相對于理想情況下的構(gòu)像所產(chǎn)生的位置差異稱像點(diǎn)位移

引起原因：1.像片傾斜引起的像點(diǎn)位移

2.地形起伏引起的像點(diǎn)位移

像片的內(nèi)方位元素：攝影物鏡后節(jié)點(diǎn)與像片之間相互位置的參數(shù)

像片外方位元素：已建立的攝影光束，確定像片攝影瞬間在地面直角坐標(biāo)系中空間位置和姿態(tài)的參數(shù)

第五篇：攝影測量

像平面坐標(biāo)系:用以表示像點(diǎn)在像平面上的位置，通常采用右手坐標(biāo)系x,y軸的選擇按需要而定，在解析和數(shù)字?jǐn)z影測量中，常根據(jù)框標(biāo)來確定像平面坐標(biāo)系，稱為像框標(biāo)坐標(biāo)系。像空間坐標(biāo)系:以攝影中心S為坐標(biāo)原點(diǎn)，x，y軸與像平面坐標(biāo)系的x，y軸平行，z軸與主光軸重合，形成像空間右手直角坐標(biāo)系S-xyz。

像空間輔助坐標(biāo)系:此坐標(biāo)系的原點(diǎn)仍選在攝影中心S，坐標(biāo)軸系的選擇視需要而定，通常有三種選取方法。其一是取鉛垂方向?yàn)閦軸，航向?yàn)閄軸，構(gòu)成右手直角坐標(biāo)系，見圖(a)。其二是以每條航線內(nèi)第一張像片的像空間坐標(biāo)系作為像空間輔助坐標(biāo)系，見圖(b)。其三是以每個(gè)像片對的左片攝影中心為坐標(biāo)原點(diǎn)，攝影基線方向?yàn)閄軸，以攝影基線及左片主光軸構(gòu)成的面作為XZ平面，構(gòu)成右手直角坐標(biāo)系，如圖(c)。用S-XYZ表示。

地面測量坐標(biāo)系：空間大地坐標(biāo)基準(zhǔn)下的高斯-克呂克6°帶或3°帶投影的平面直角坐標(biāo)與定義的從某一基準(zhǔn)面量起的高程，兩者結(jié)合而成的空間左手直角坐標(biāo)系。

攝影測量坐標(biāo)系：是指描述攝影測量模型的空間直角坐標(biāo)系。其原點(diǎn)選在某攝站或某一已知點(diǎn)，X軸大體與航線方向一致，Z軸與鉛垂線方向一致且向上為正的右旋空間直角坐標(biāo)系。像點(diǎn)位移：航空像片是地面的中心投影，根據(jù)中心投影的原理，無論是帶有起伏狀態(tài)的地形，還是高出地面的任何物體，反映到航空像片上的像點(diǎn)與其平面位置相比，一般都會產(chǎn)生位置的移動(dòng)，這種像點(diǎn)位置的移動(dòng)，叫做像點(diǎn)位移。主要是由像片傾斜、地面點(diǎn)相對于基準(zhǔn)面的高差和物理因素（如攝影材料變形、壓平誤差、攝影物鏡畸變、大氣折光和地球曲率等）產(chǎn)生。

為什么研究像點(diǎn)位移的規(guī)律可以清楚知道透鏡成像的大小，虛實(shí)6-=等基本情況，對于應(yīng)用透鏡解決生活和實(shí)際問題是必不可少的像點(diǎn)位移的規(guī)律：

單向空間后方交會：已知至少3個(gè)地面控制點(diǎn)的坐標(biāo)A，B，C，與其影像上對應(yīng)的三個(gè)像點(diǎn)的影像坐標(biāo)a，b，c，根據(jù)共線方程，反求該像片的6個(gè)外方位元素。

立體像對雙像前方交會：現(xiàn)已知這兩張像片的內(nèi)外方位元素，設(shè)想將該像片的內(nèi)外方位元素值置于攝影時(shí)的位置，顯然同名射線S1a1和S2a2必然交于地面點(diǎn)A。這種由立體像對中兩張像片的內(nèi)外方位元素和像點(diǎn)坐標(biāo)來確定相應(yīng)地面點(diǎn)的地面坐標(biāo)的方法。

立體像對雙像前方交會的目的：利用立體像對上的同名像點(diǎn)，才能得到兩條同名射線在空間相交的點(diǎn)。

解析法相對定向：通過計(jì)算相對定向元素建立地面立體模型。

解析法相對定向的共面條件：B·（S1a1 X S2a2）=0

模型連接：將單個(gè)模型連接成為航帶模型,要將各模型的不同的比例尺歸化為統(tǒng)一的比例尺,通常以相鄰像對重疊范圍內(nèi)三個(gè)連接點(diǎn)的高程應(yīng)相等為條件,建立統(tǒng)一的以第一個(gè)模型的比例尺為基準(zhǔn)的航帶模型.模型連接的實(shí)質(zhì)是求出相鄰模型間的比例歸化系數(shù).航帶網(wǎng)整體平差的實(shí)質(zhì)是以一條航帶模型為平差單元,解求航帶的非線性改正系數(shù),即多項(xiàng)式系數(shù).