后方交会程序实现(c语言版)

合集下载

空间后方交会程序设计

单片空间后方交会程序设计专业：测绘工程班级：测绘101姓名：张雯学号：2010220200512 用C 或VC++语言实现单片后方交汇的计算。

以单幅影像为基础，从该影象所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，求解该影象在航空摄影时刻的像片外方位元素Xs ，Ys ，Zs,ф,ω,κ.共线条件方程如下：x-x0=-f*[a1(X-Xs)+b1(Y-Ys)+c1(Z-Zs)]/[a3(X-Xs)+b3(Y-Ys)+c3(Z-Zs)] y-y0=-f*[a2(X-Xs)+b2(Y-Ys)+c2(Z-Zs)]/[a3(X-Xs)+b3(Y-Ys)+c3(Z-Zs)] x,y 为像点的像平面坐标；x0,y0,f 为影像的外方位元素；Xs,Ys,Zs 为摄站点的物方空间坐标；X,Y,Z 为物方点的物方空间坐标；输入原始数据归算像点坐标x,,y 计算和确定初值X s 0, Y s 0, Z s 0, φ0,ω0,κ0 组成旋转矩阵R 计算（x ），（y ）和l x ，l y 逐点组成误差方程并法化所有点完否？解法方程，求未知数改正数计算改正后外方位元素未知数改正数<限差否？整理并输出计算结果正常输出迭代次数小于限差否？输出中间结果和出错信息4个地面控制点的地面坐标及其对应像点的像片坐标：结果输出：已知条件:像点坐标x,y:-53.4 82.21-14.78 -76.6310.46 64.43153.24 -86.15地点坐标Xa,Ya,Za:37631.1 31324.5 728.6939101 24935 2386.540426.5 30319.8 757.31-68.99 36589.4 25273.3摄影比例尺：1:50000内方位元素:x0=y0=0 f=153.24mm计算结果:旋转矩阵：0.997709 0.0675334 0.00398399-0.0675254 0.997715 -0.00211178-0.0041175 0.00183792 0.99999像点坐标位:(单位：mm)-86.15 -68.99-53.41 82.21-14.78 -76.6310.47 64.43外方位元素：Xs=39795.435 精度为:1.1254813Ys=27476.479 精度为:1.2437625Zs=7572.6929 精度为:0.48380521q=-0.0039840098 精度为:0.00018182003 w=0.0021117837 精度为:0.00015959235 k=-0.067576934 精度为:7.2440432e-005 迭代次数：4Press any key to continue。

摄影测量学单像空间后方交会编程实习报告(精品资料).doc

【最新整理，下载后即可编辑】摄影测量学单像空间后方交会编程实习报告班级：130x姓名：xx学号：2013302590xxx指导老师：李欣一、实习目的通过对提供的数据进行计算，输出像片的外方位元素并评定精度。

深入理解单像空间后方交会的思想，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。

通过尝试编程实现加强编程处理问题的能力和对实习内容的理解，通过对实验结果的分析，增强综合运用所学知识解决实际问题的能力。

了解摄影测量平差的基本过程，掌握空间后方交会的定义和实现算法。

二、实习内容根据学习的单像空间后方交会的知识，用程序设计语言（C++或C语言）编写一个完整的单像空间后方交会程序，通过对提供的数据进行计算，输出像片的外方位元素并评定精度。

三、实习数据已知航摄仪的内方位元素：fk =153.24mm，x=y=0，摄影比例尺为1:15000；4个地面控制点的地面坐标及其对应像点的像片坐标：四、实习原理如果我们知道每幅影像的6个外方位元素，就能确定被摄物体与航摄影像的关系。

因此，如何获取影像的外方位元素，一直是摄影测量工作者所探讨的问题。

可采取的方法有：利用雷达、全球定位系统（GPS）、惯性导航系统（INS）以及星相摄影机来获取影像的外方位元素；也可以利用影像覆盖范围内一定数量的控制点的空间坐标与摄影坐标，根据共线条件方程，反求该影像的外方位元素，这种方法称为单幅影像的空间后方交会。

单像空间后方交会的基本思想是：以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标量测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素Xs,Ys,Zs,ϕ,ω,κ。

五、实习流程1. 获取已知数据。

从摄影资料中查取影像比例尺1/m ，平均摄影距离（航空摄影的航高、内方位元素x 0，y 0，f ；获取控制点的空间坐标X t ，Y t ，Z t 。

2. 量测控制点的像点坐标并进行必要的影像坐标系统误差改正，得到像点坐标。

单片空间后方交会程序设计

单片空间后方交会程序设计
1 目的
用程序设计语言（VC或者VB）编写一个完整的单片空间后方交会程序，通过对提供的试验数据进行计算，输出像片的外方位元素并评定精度。

本实验的目的在于让学生深入理解单片空间后方交会的原理，体会在有多余观测情况下，用最小二乘平差方法编程实现解求影像外方位元素的过程。

通过上机调试程序加强动手能力的培养，通过对实验结果的分析，增强学生综合运用所学知识解决实际问题的能力。

2. 内容
利用一定数量的地面控制点，根据共线条件方程求解像片外方位元素。

3. 数据准备
已知航摄仪的内方位元素：f
k =153.24mm，x
=y
=0.0mm，摄影比例尺为1:50000；
4个地面控制点的地面坐标及其对应像点的像片坐标：
4. 操作步骤
上机调试程序并打印结果。

摄影测量后方交会VC实现代码及实习报告

摄影测量课间实习单片空间后方交会班级：09031姓名：吴煜晖学号：20093025901232011-10-8一、实习原理单像空间后方交会的基本思想是：以单幅影像为基础，从该影像所覆盖地面范围内若干控制点的已知地面坐标和相应点的像坐标测值出发，根据共线条件方程，解求该影像在航空摄影时刻的外方位元素XS ，YS，ZS，ψ，ω，κ。

由于空间后方交会所采用的数学模型共线方程是非线性函数，为了便于外方位元素的解求，须首先对共线方程进行线性化，然后进行后方交会，最后在精度评定。

二、实习过程1、实习所用数据本次实习数据采用课本P.44 27题所给数据。

如下图：2、程序流程图及界面设计本程序程序框图如下：本人采用Visual C++6.0编此程序，利用MFC来设计程序。

主程序页面设计如下：子窗口（即进行计算后所得结果页面）设计如下：3、程序代码本程序代码较多，在此讲部分重要代码列出，其余代码参见程序源代码。

对话框类头文件内声明的类成员及函数（后来增加的）：public:void houfangjiaohui();double fi,w,k; //影像外方位角double Xs0,Ys0,Zs0; //后方交会所求解double a[3],b[3],c[3];double mx[6],m0;void inv(double *a,int n); /*正定矩阵求逆*/void transpose(double *m1, double *m2, int m, int n); //矩阵转置void mult(double *m1, double *m2, double *result, int i_1, int j_12, int j_2); //矩阵相乘部分核心代码，在此讲后方交会计算代码给出：void CFinalworkDlg::houfangjiaohui(){double t_fk;fi=w=k=0.0;int i; //中间变量double t_x[4],t_y[4],t_X[4],t_Y[4],t_Z[4];//将单位换成米t_x[0]=m_x1/1000.0; t_x[1]=m_x2/1000.0; t_x[2]=m_x3/1000.0;t_x[3]=m_x4/1000.0;t_y[0]=m_y1/1000.0; t_y[1]=m_y2/1000.0; t_y[2]=m_y3/1000.0;t_y[3]=m_y4/1000.0;t_X[0]=m_X1; t_X[1]=m_X2; t_X[2]=m_X3; t_X[3]=m_X4;t_Y[0]=m_Y1; t_Y[1]=m_Y2; t_Y[2]=m_Y3; t_Y[3]=m_Y4;t_Z[0]=m_Z1; t_Z[1]=m_Z2; t_Z[2]=m_Z3; t_Z[3]=m_Z4;Xs0=Ys0=Zs0=0.0;for(i=0;i<4;i++){Xs0+=t_X[i];Ys0+=t_Y[i];}//确定未知数初始值Xs0/=4;Ys0/=4;t_fk=m_fk/1000.0;Zs0=t_fk*m_mk;double x00=m_x00;double y00=m_y00;double A[8*6];double AT[6*8];double ATA[6*6];double L[8];double ATL[6*1];double Xo[4],Yo[4],Zo[4],Xom,Yom,Zom;double Delta[6];while(1){a[0]=cos(fi)*cos(k)-sin(fi)*sin(w)*sin(k);a[1]=-cos(fi)*sin(k)-sin(fi)*sin(w)*cos(k);a[2]=-sin(fi)*cos(w);b[0]=cos(w)*sin(k);b[1]=cos(w)*cos(k);b[2]=-sin(w);c[0]=sin(fi)*cos(k)+cos(fi)*sin(w)*sin(k);c[1]=-sin(fi)*sin(k)+cos(fi)*sin(w)*cos(k);c[2]=cos(fi)*cos(w);for(i=0;i<4;i++){Xom=a[0]*(t_X[i]-Xs0)+b[0]*(t_Y[i]-Ys0)+c[0]*(t_Z[i]-Zs0);Yom=a[1]*(t_X[i]-Xs0)+b[1]*(t_Y[i]-Ys0)+c[1]*(t_Z[i]-Zs0);Zom=a[2]*(t_X[i]-Xs0)+b[2]*(t_Y[i]-Ys0)+c[2]*(t_Z[i]-Zs0);Xo[i]=-t_fk*Xom/Zom;Yo[i]=-t_fk*Yom/Zom;Zo[i]=Zom;A[i*12+0]=(a[0]*t_fk+a[2]*(t_x[i]-x00))/Zo[i];A[i*12+1]=(b[0]*t_fk+b[2]*(t_x[i]-x00))/Zo[i];A[i*12+2]=(c[0]*t_fk+c[2]*(t_x[i]-x00))/Zo[i];A[i*12+3]=(t_y[i]-y00)*sin(w)-cos(w)*(t_fk*cos(k)+((t_x[i]-x00)*cos(k)-(t _y[i]-y00)*sin(k))*(t_x[i]-x00)/t_fk);A[i*12+4]=-t_fk*sin(k)-(t_x[i]-x00)*((t_x[i]-x00)*sin(k)+(t_y[i]-y00)*cos(k ))/t_fk;A[i*12+5]=t_y[i]-y00;A[i*12+6]=(a[1]*t_fk+a[2]*(t_y[i]-y00))/Zo[i];A[i*12+7]=(b[1]*t_fk+b[2]*(t_y[i]-y00))/Zo[i];A[i*12+8]=(c[1]*t_fk+c[2]*(t_y[i]-y00))/Zo[i];A[i*12+9]=-(t_x[i]-x00)*sin(w)-((t_y[i]-y00)*((t_x[i]-x00)*cos(k)-(t_y[i]-y0 0)*sin(k))/t_fk-t_fk*sin(k))*cos(w);A[i*12+10]=-t_fk*cos(k)-(t_y[i]-y00)*((t_x[i]-x00)*sin(k)+(t_y[i]-y00)*cos(k))/t_fk;A[i*12+11]=-(t_x[i]-x00);L[i*2+0]=(t_x[i]-x00-Xo[i]);L[i*2+1]=(t_y[i]-y00-Yo[i]);}transpose(A,AT,8,6);mult(AT,A,ATA,6,8,6);mult(AT,L,ATL,6,8,1);inv(ATA,6);mult(ATA,ATL,Delta,6,6,1);Xs0+=Delta[0];Ys0+=Delta[1];Zs0+=Delta[2];fi+=Delta[3];w+=Delta[4];k+=Delta[5];if((fabs(Delta[0])>0.00001||fabs(Delta[1])>0.00001||fabs(Delta[2])>0.00001|| fabs(Delta[3])>0.00001||fabs(Delta[4])>0.00001||fabs(Delta[5])>0.00001)==0) {double V[8];double s=0.0;mult(A,Delta,V,8,6,1);for(i=0;i<8;i++){V[i]-=L[i];s+=V[i]*V[i];}m0=sqrt(s/2.0); //单位权中误差for(i=0;i<6;i++) //精度评定{mx[i]=sqrt(ATA[7*i])*m0;}break;}}}三、实习结果程序主界面如下：在该界面中，我们可以手动输入各项数据，也可以清空所有数据以及将所有空恢复默认数据（即题目所给数据）。

利用编程计算器进行距离后方交会的严密平差计算

本论述采用测量业中广泛使用的ＣＡＳＩＯ编程计算器，按条件平差的方法进行严密平差，计算出交会点的坐标，并进行精度检核与评定。在待测点Ｐ处安置仪
器，观测Ｐ点到左边Ｌ点的距离ＪＵＬＩ（Ｌ）、中间Ｍ点的距离ＪＵＬＩ（Ｍ）和右边Ｒ点的距离ＪＵＬＩ（Ｒ）。然后根据
Ｌ点的坐标ｘ（Ｌ）、Ｙ（Ｌ），Ｍ点的坐标ｘ（Ｍ）、Ｙ（Ｍ）
和Ｒ点的坐标ｘ（Ｒ）、Ｙ（Ｒ），就可以平差计算出Ｐ点
的坐标ｘｐ、Ｙｐ，并进行精度评定。
信息技术
２０１３年（第４２卷）第３期
利用编程计算器进行距离后方交会的严密平差计算
孙芳琳
（兰州城市建设学校，甘肃兰州７３００４６）
摘
要：测量工作中可采用后方交会法进行观测。早期经常采用角度后方交会法进行观测，但精度低，测量位中误差ＤＩＡＮＷＥＩ
ＷＵＣＨＡ及交会点Ｐ的坐标Ｘｐ、Ｙｐ。
５算例
已知三点Ｌ（４２１２．７６３，６６１５．３９０）、Ｍ（４６２３．８７５，５４４３．２８５）、Ｒ（５７４２．２８５，５８３５．４１０）。测得Ｐ点到左
少使用。随着测距仪器的普及，距离后方交会法在测量工作中已经被广泛使用。当采用已知两个点的普通距离后方交会进行观测，其计算比较简单，常用的全站仪上也都有计算程序，但两点的距离后方交会没有必要的检核，容易出现错误，点位误差大小也无法得知。如果采用已知三个点的距离后方交会，利用编程计算器进行严密平差计算，不但可以检查出错误，而且可以计算出点位误差的大小，能够使观测成果精确可靠。关键词：控制点；后方交会；编程计算器；严密平差；精度评定

武汉大学《摄影测量学》复习题库

2. 内容
熟悉 1818 立体坐标量测仪的基本结构，立体观察，坐标量测。左右视差（p）读数鼓
上下视差（q）读数鼓
x 读数鼓
x 手轮 y 手轮
3. 资料准备
一个 18cm×18cm 的立体影像对
y 读数鼓
左右视差手轮上下视差环
左像片
右像片
4. 操作步骤
仪器归零：各个手轮应放在零读数位置上，左、右测标分别对准左、右像片盘的中心即仪器坐标系与像片坐标系重合。
Z 2195.17
728.69 2386.50
757.31
4. 操作步骤
上机调试程序并打印结果。
“POS 辅助光束法平差系统 WuCAPS”使用
1. 目的
通过参观 POS 辅助光束法区域网平差程序系统 WuCAPS,使学生初步了解摄影测量区域网平差的基本功能和一般作业流程。
2. 内容
指导教师讲解摄影测量区域网平差的基本概念、主要功能及一般作业流程。学生按照要求，完成一些简单的操作，例如，内定向、相对定向、绝对定向、航带法区域网平差、光束法区域网平差、 GPS 辅助光束法区域网平差、POS 辅助光束法区域网平差等。
像片定向：移动 x 手轮，单眼观察测标的移动看是否沿像片上的 x 轴向运动，若测标不在 x 轴向上，则需要用κ 螺旋旋转像片，使测标保持在 x 轴上移动。
坐标量测：移动 x，y，p，q 手轮，使测标立体切准量测像点，并记下相应读数鼓上的读数。坐标计算：x1=x-x0，y1=y-y0，x2=x1-(p-p0)，y2=y1-(q-q0)，其中，x0，y0，p0，q0为仪器零位置。
的方法？
5. 什么是共线条件方程式？试推导其数学表达式，并说明它在摄影测量中的应
用。

后方交会程序实现c语言版.doc

#include<stdio.h> #include<stdlib.h>#include<malloc.h>#include<math.h>#include<conio.h>#define N 4#define T 1.41421void turn(double *A,double A2[],int m,int n) // 计算矩阵的转置{ int i,j; for(i=0;i<m;i++)for(j=0;j<n;j++)A2[j*m+i]=A[i*n+j];}void mulAB(double *A,double *B,double *C,int am,int an,int bm,int bn) //计算两矩阵相乘{ int i,j,l,u;if(an!=bm){ printf("error!cannot do the multiplication.\n");return; }for(i=0;i<am;i++)for(j=0;j<bn;j++) {u=i*bn+j;C[u]=0.0;for(l=0;l<an;l++)C[u]+=A[i*an+l]*B[l*bn+j];} return;}double *inv(double *a,int n) // 计算矩阵的逆，本程序的难,采用高斯约旦全选主元法点{int *is,*js,i,j,k,l,u,v;double d,p;is=(int*)malloc(n*sizeof(int));js=(int*)malloc(n*sizeof(int));for (k=0; k<=n-1; k++){d=0.0;for (i=k;i<n;i++) for (j=k;j<n;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("error not inv\n");return NULL;}if (is[k]!=k)for (j=0;j<n;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;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;j++) if (j!=k){u=k*n+j;a[u]=a[u]*a[l];}for (i=0;i<n;i++)if (i!=k)for (j=0;j<n;j++)if (j!=k)a[u]=a[u]-a[i*n+k]*a[k*n+j];}for (i=0;i<n;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;a[u]=a[v];a[v]=p;}if (is[k]!=k)for (i=0;i<n;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 a;}main()// 主函数，空间后方交会的计算{FILE *fp; // 定义一个文件指针long m;int i,j=0,it;double G[1000];doublef,t,w,k,limit=0,S1=0.0,S2=0.0,S3=0.0,x[N]={0},y[N]={0},x0[N]={0},y0[N]={0},X[N]={0},Y[N] ={0},Z[N]={0},Xs0,Ys0,Zs0;doublea[3],b[3],c[3],A[2*N*6],AT[6*2*N],ATA[6*6],*ATA_=NULL,l[2*N],ATl[6],V[6]={0};double F[6],Qx[6][6],Mi[6][6];if((fp=fopen("e:\\shuju.txt","r"))==NULL) // 使 fp 指向被打开的文件{ // 并判断文件是否打开成功printf("\nerror on open shuju.txt\n");getch();exit(1);} for(i=0;i<N;i++)fscanf(fp,"%lf%lf%lf%lf%lf",&x[i],&y[i],&X[i],&Y[i],&Z[i]); 件中的数据赋给主函数定义的变量printf(" 原始数据 :\n");printf("x\t\ty\t\tX\t\tY\t\tZ\t\t\n");// 输出文件中的原始数据for(i=0;i<N;i++)printf("%.3lf\t\t%.3lf\t\t%.3lf\t%.3lf\t%.3lf\n",x[i],y[i],X[i],Y[i],Z[i]);printf("\n 请输入摄影机主距和摄影比例尺分母； f,m:");scanf("%lf%ld",&f,&m); //输入 f,mf=f/1000.0;for(i=0;i<N;i++){ x[i]/=1000.0;y[i]/=1000.0;S1+=X[i];S2+=Y[i];S3+=Z[i];}Xs0=S1/N;Ys0=S2/N;Zs0=m*f+S3/N;t=0.0;w=0.0;k=0.0; // 计算外方位元素的初始值while(1)shuju.txt// 将文printf("\n ------------------------------- 第%d 次计算----------------------------- \n",j+1);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); // 计算旋转矩阵for(i=0;i<N;i++){x0[i]=-f*(a[0]*(X[i]-Xs0)+b[0]*(Y[i]-Ys0)+c[0]*(Z[i]-Zs0))/(a[2]*(X[i]-Xs0)+b[2]*(Y[i]-Ys0)+c [2]*(Z[i]-Zs0));y0[i]=-f*(a[1]*(X[i]-Xs0)+b[1]*(Y[i]-Ys0)+c[1]*(Z[i]-Zs0))/(a[2]*(X[i]-Xs0)+b[2]*(Y[i]-Ys0)+c[2]*(Z[i]-Zs0));//计算像点坐标近似值G[i]=a[2]*(X[i]-Xs0)+b[2]*(Y[i]-Ys0)+c[2]*(Z[i]-Zs0);}for(i=0;i<N;i++){A[i*12+0]=(a[0]*f+a[2]*x[i])/G[i];A[i*12+1]=(b[0]*f+b[2]*x[i])/G[i];A[i*12+2]=(c[0]*f+c[2]*x[i])/G[i];A[i*12+3]=y[i]*sin(w)-(x[i]*(x[i]*cos(k)-y[i]*sin(k))/f+f*cos(k))*cos(w);A[i*12+4]=-f*sin(k)-x[i]*(x[i]*sin(k)+y[i]*cos(k))/f;A[i*12+5]=y[i];A[i*12+6]=(a[1]*f+a[2]*y[i])/G[i];A[i*12+7]=(b[1]*f+b[2]*y[i])/G[i];A[i*12+8]=(c[1]*f+c[2]*y[i])/G[i];A[i*12+9]=-x[i]*sin(w)-(y[i]*(x[i]*cos(k)-y[i]*sin(k))/f-f*sin(k))*cos(w);A[i*12+10]=-f*cos(k)-y[i]*(x[i]*sin(k)+y[i]*cos(k))/f;A[i*12+11]=-x[i]; l[i*2+0]=x0[i]-x[i];l[i*2+1]=y0[i]-y[i]; //计算误差方程的系数阵以及lx,ly}// printf("output matrix: A\n");// printmatrix(A,2*N,6);// printf("output matrix: l\n");// printmatrix(l,2*N,1);turn(A,A T,2*N,6); // 计算AT// printf("output matrix: AT\n");// printmatrix(AT,6,2*N);mulAB(AT,A,ATA,6,2*N,2*N,6); // 计算ATA ，组法方程ATA_=inv(ATA,6); // 计算ATA 的逆，中间量int p;int cnt=-1; for(it=0;it<36;it++){p=it%6;if(it%6==0){cnt++;}Qx[cnt][p++]=A TA_[it];}for(int it=0;it<6;it++){for(int jt=0;jt<6;jt++){ if(it!=jt){Qx[it][jt]=0;// 提取Qx 的主对角线元素=Qii}// printf("%-10.3lf ",Qx[it][jt]);}// printf("\n");}mulAB(A T,l,ATl,6,2*N,2*N,1); // 计算常数项ATL// printf("outpit matrinx: ATl\n");// printmatrix(ATl,6,1);mulAB(A TA_,ATl,V,6,6,6,1); // 解法方程，求改正数，//printf("output matrix: V\n");// printmatrix(V ,6,1);Xs0+=V[0];Ys0+=V[1];Zs0+=V[2];t+=V[3];w+=V[4];k+=V[5];for(i=0;i<6;i++){F[i]=V[i]/T;//m0}printf(" 第%d 次计算的外方位元素为：\n",++j); printf("Xs=%.5lf\tYs=%.5lf\tZs=%.5lf\nt=%.5lf\tw=%.5lf\tk=%.5lf\n",Xs0,Ys0,Zs0,t,w,k);if(Xs0-limit<=0.0001&&Xs0-limit>=-0.0001){ // 控制迭代次数break;}limit=Xs0;}printf("\n 外方位元素为：\n");printf("Xs=%.5lf\tYs=%.5lf\tZs=%.5lf\nt=%.5lf\tw=%.5lf\tk=%.5lf\n",Xs0,Ys0,Zs0,t,w,k);printf("\n 中误差为:\n");for(i=0;i<6;i++){for(j=0;j<6;j++){开根号Mi[i][j]=F[i]*(sqrt(Qx[i][j]));//mi=m0*Qiiprintf("%-13.10lf",Mi[i][j]);}printf("\n");}fclose(fp);return0; }。

C#空间后方交会

double f = 153.24;
int N;
public string filpath = null;
public string filename = null;
public Form1()
{
InitializeComponent();
}
private void button1_Click(object sender, EventArgs e)
MessageBox.Show("迭代超过20次，迭代失败"); break; }
Math.Sin(K); Math.Cos(K);
Math.Sin(K);
else { t0 = jisuan.X1[0];
t1 =jisuan.X1[1]; t2 = jisuan.X1[2]; t3=jisuan.X1 [3]; t4=jisuan.X1 [4]; t5=jisuan.X1 [5]; Xs =Xs+t0; Ys =Ys+t1; Zs =Zs+t2; F = F + t3; W = W + t4; K = K + t5; FWK[0, 0] = Math.Cos(F) * Math.Cos(K) - Math.Sin(F) * Math.Sin(W) *
public partial class Form1 : Form {
double []X=new double [10]; double []Y=new double [10]; double []Z=new double [10];
double []x=new double [10]; double []y=new double [10]; double Xs=0, Ys=0, Zs=0,m; double[,] FWK = new double [,] {{1,0,0},{0,1,0},{0,0,1}}; double F=0, W=0, K=0;

(空间后方交会的计算过程)空间后方交会 PPT

3）确定未知数的初始值：在竖直摄影的情况下，角元
素的初始值为0,即 00k00 ；线元素
中用，四个ZS角0 上H的m 控f制，X点S0,坐YS0标取的值平可均值，即：XS0
1 4
4 i1
Xtpi
大家好
YS0
1 4
4
Y15 tpi
i1
4) 计算旋转矩阵R：利用角元素的近似值计算
方向元素，组成旋转矩阵R。
或写成： V x a 1d 1S X a 1d 2 S Y a 1d 3S Z a 1d 4 a 1d 5 a 1d 6 lk x V y a 2 d 1S X a 2d 2 S Y a 2d 3S Z a 2d 4 a 2d 5 a 2d 6 l k y
其中：
llx y x y ((x y)) xy ffa a a a 3 12 3 ((((X X X X A A A A X X X X S SS S)))) b b b b 1 33 2 ((Y ((Y Y Y A A A A Y Y Y Y S SS S )))) c cc c 1 3(3 2 (Z ((Z Z Z A A A A Z Z Z Z S SS S ))))
大家好
14
空间后方交会的计算过程
1）获取已知数据：从摄影资料中查取像片比例尺，平
均航高，内方位元素；从外业测量成果中，获取控制点的地面测量坐标并转换为地面摄影测量坐标。
2）量测控制点的坐标：将控制点标刺在像片上，利用
立体坐标量测仪量测控制点的像框标坐标系坐标，
并经像主点坐标改正，得到像点坐标x,y；
空间后方交会的数学模型是共线方程，即中心投影的构像方程：
xf yf
a1(XAXS)b1(YAYS)c1(ZAZS) a3(XAXS)b3(YAYS)c3(ZAZS) a2(XAXS)b2(YAYS)c2(ZAZS)

单片空间后方交会C#源代码

单片空间后方交会C#源代码主方法：private void Cal_Click(object sender, EventArgs e) {string[] lines = RichText.Text.Split('\n');long m = lines.Length;m = m - 1;//真实数据行数double[] Coor_x = new double[m];//已知点x坐标double[] Coor_y = new double[m];//已知点x坐标double[] Coor_X = new double[m];//已知点X坐标double[] Coor_Y = new double[m];//已知点Y坐标double[] Coor_Z = new double[m];//已知点Z坐标///赋值for (int i = 0; i < m; i++){string[] FJstring = Regex.Split(lines[i+1], ",");Coor_x[i] = 0.001*(Convert.T oDouble(FJstring[0])); Coor_y[i] = 0.001 *( Convert.T oDouble(FJstring[1])); Coor_X[i] = Convert.ToDouble(FJstring[2]);Coor_Y[i] = Convert.ToDouble(FJstring[3]);Coor_Z[i] = Convert.ToDouble(FJstring[4]);}if (textBox_m.Text == ""){MessageBox.Show("请输入参数！");}if (textBox_m.Text != ""){double M = double.Parse(textBox_m.Text);//比例尺double f = 0.001 * (double.Parse(textBox_f.Text));//焦距double x0 = 0.001 * double.Parse(textBox_x0.Text);//内方位元素x0double y0 = 0.001 * double.Parse(textBox_y0.Text);//内方位元素y0double X0 = 0, Y0 = 0, Z0 = 0;//外方位坐标元素初始值double min = (double.Parse(textBox_k.Text));//焦距double angle1 = 0, angle2 = 0, angle3 = 0;//外方位角元素初始值for (int i = 0; i < m; i++){//累加X0 = Coor_X[i] + X0;Y0 = Coor_Y[i] + Y0;Z0 = Coor_Z[i] + Z0;}X0 = X0 / m;Y0 = Y0 / m;Z0 = Z0 / m + M * f;double[,] A = new double[2 * m, 6];//A，二维：（2 * m）*6 double[,] l = new double[2 * m, 1];//l，二维：（2 * m）*1double[,] V = new double[2 * m, 6];//V，二维：（2 * m）*6 double[,] XXX = new double[6, 1];//XXX，6*1 所求参数值，即外方位元素改正数///循环解求误差方程，以限差为标准结束循环Matrix ad = new Matrix();Cls_AL ae = new Cls_AL();int num = 0;//记录循环次数for (int i = 0; ; i++){//求A，lae.Cal_AL(X0, Y0, Z0, angle1, angle2, angle3, m, f, x0, y0, Coor_x, Coor_y, Coor_X, Coor_Y, Coor_Z, A, l);//求XXXXXX = ad.Cal_Chengfa(ad.Cal_Chengfa(ad.Cal_Qiuni(ad.Cal_Chengfa(ad .Cal_Zhuanzhi(A), A)), ad.Cal_Zhuanzhi(A)), l);///改正外方位元素X0 = X0 + XXX[0, 0];Y0 = Y0 + XXX[1, 0];Z0 = Z0 + XXX[2, 0];angle1 = angle1 + XXX[3, 0];angle2 = angle2 + XXX[4, 0];angle3 = angle3 + XXX[5, 0];num++;if ((Math.Abs(XXX[3, 0]) < min) && (Math.Abs(XXX[4, 0]) < min) && (Math.Abs(XXX[5, 0]) < min)){ break; }}///精度评定;V = ad.Cal_Jianfa(ad.Cal_Chengfa(A, XXX), l);double[,] c = new double[1, 1];//单位权中误差c = ad.Cal_Chengfa(ad.Cal_Zhuanzhi(V), V);double cc = Math.Sqrt(c[0, 0] / (2 * m - 6));double[,] AA = new double[6, 6];//外方位元素改正数的协因数阵QQ AA = ad.Cal_Qiuni(ad.Cal_Chengfa(ad.Cal_Zhuanzhi(A), A));RichText.Text = ("外方位元素为：" + '\n' + "Xs = " + Convert.T oString(X0) + '\n' + "Ys = " + Convert.ToString(Y0) + '\n' + "Zs = " + Convert.T oString(Z0) + '\n' + "angel_1 = " + Convert.T oString(angle1) + '\n' + "angel_2 = " + Convert.T oString(angle2) + '\n' + "angel_3 = " +Convert.T oString(angle3) + '\n' + "单位权中误差为：" + Convert.T oString(cc) + '\n' + "Xs精度为： " + Convert.T oString(cc * Math.Sqrt(AA[0, 0])) + '\n' + "Ys精度为：" + Convert.ToString(cc * Math.Sqrt(AA[1, 1])) + '\n' + "Zs精度为：" + Convert.ToString(cc * Math.Sqrt(AA[2, 2])) + '\n' + "迭代次数为：" + Convert.T oString(num));}}动态类库：using System;using System.Collections.Generic;using System.Linq;using System.Text;namespace dll{////// Matrix类包含三个方法：矩阵乘法，矩阵转置，矩阵求逆///。

C语言空间后方交会源代码

#include<stdio.h>#include<math.h>#define n 4 //控制点个数#define PI 3.14159265struct coordinate{double x; //像点坐标double y;double Xt; //控制点坐标double Yt;double Zt;};// void inverse(double c[6][6]) //矩阵求逆// {// int i,j,h,k;// double p;// double q[6][12];// for(i=0;i<6;i++)//构造高斯矩阵// for(j=0;j<6;j++)// q[i][j]=c[i][j];// for(i=0;i<6;i++)// for(j=6;j<12;j++)// {// if(i+6==j)// q[i][j]=1;// else// q[i][j]=0;// }// for(h=k=0;k<n-1;k++,h++)//消去对角线以下的数据// for(i=k+1;i<6;i++)// {// if(q[i][h]==0)// continue;// p=q[k][h]/q[i][h];// // p=q[i][h]/q[k][h];// for(j=0;j<12;j++)// {// q[i][j]*=p;// q[i][j]-=q[k][j];// }// }// for(h=k=5;k>0;k--,h--) // 消去对角线以上的数据// for(i=k-1;i>=0;i--)// {// if(q[i][h]==0)// continue;// p=q[k][h]/q[i][h];// // p=q[i][h]/q[k][h];// for(j=11;j>0;j--)// {// q[i][j]*=p;// q[i][j]-=q[k][j];// }// }// for(i=0;i<6;i++)//将对角线上数据化为1// {// p=1.0/q[i][i];// for(j=0;j<12;j++)// q[i][j]*=p;// }// for(i=0;i<6;i++) //提取逆矩阵// for(j=0;j<n;j++)// {// c[i][j]=q[i][j+6];// }// }void ContraryMatrix(double *const pMatrix, double *const _pMatrix, const int &dim) {int ii,jj,kk;int flag=0;double *tMatrix = new double[2*dim*dim];for (int i=0; i<dim; i++){for (int j=0; j<dim; j++)tMatrix[i*dim*2+j] = pMatrix[i*dim+j];}for (i=0; i<dim; i++){for (int j=dim; j<dim*2; j++)tMatrix[i*dim*2+j] = 0.0;tMatrix[i*dim*2+dim+i] = 1.0;}//Initialization over!for (i=0; i<dim; i++)//Process Cols{double base = tMatrix[i*dim*2+i];if (fabs(base) < 1E-300){for (ii=i;ii<dim;ii++){if(tMatrix[ii*dim*2+i]!=0){flag=1;for (jj=0;jj<2*dim;jj++){tMatrix[i*dim*2+jj]+=tMatrix[ii*dim*2+jj];}}}base = tMatrix[i*dim*2+i];if (flag==0){printf( "求逆矩阵过程中被零除，无法求解!") ;}// exit(0);}for (int j=0; j<dim; j++)//row{if (j == i) continue;double times = tMatrix[j*dim*2+i]/base;for (int k=0; k<dim*2; k++)//col{tMatrix[j*dim*2+k] = tMatrix[j*dim*2+k] - times*tMatrix[i*dim*2+k];}}for (int k=0; k<dim*2; k++){tMatrix[i*dim*2+k] /= base;}}for (i=0; i<dim; i++){for (int j=0; j<dim; j++)_pMatrix[i*dim+j] = tMatrix[i*dim*2+j+dim];}delete[] tMatrix;}void main(){double f,xo,yo; //内方位元素int m; //摄影比例尺分母f=0.15324;xo=0;yo=0;m=50000;struct coordinate coor[n]={{-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.50},{0.01046,0.06443,40426.54,30319.81,757.31}};int i,j; //循环变量double Xs,Ys,Zs; //外方位线元素double phi,omega,kappa; //外方位角元素double R[3][3]; //旋转矩阵Rdouble (x)[n],(y)[n]; //控制点像点坐标的近似值double A[8][6]; //误差方程中的矩阵Adouble ATA[6][6]; //矩阵A的转置矩阵与A的乘积double L[8][1];double ATL[6][1]; //矩阵A的转置矩阵与L的乘积double _ATA[6][6];double X[6][1]; //未知数矩阵double V[8][1]; //改正数矩阵double DXs,DYs,DZs,Dphi,Domega,Dkappa; //6个外方位元素的改正值//与精度计算有关的量double m0; //单位权中误差double Q[6][6]; //协方差矩阵double m1,m2,m3,m4,m5,m6; //未知数Xs,Ys,Zs,phi,omega,kappa的中误差矩阵Xs=0;Ys=0;for(i=0;i<n;i++){Xs+=coor[i].Xt;Ys+=coor[i].Yt;}//外方位元素的初始值Xs/=n;Ys/=n;Zs=f*m;phi=0;omega=0;kappa=0; //在竖直摄影情况下do{//计算旋转矩阵R的各个元素R[0][0]=cos(phi)*cos(kappa)-sin(phi)*sin(omega)*sin(kappa);R[0][1]=-cos(phi)*sin(kappa)-sin(phi)*sin(omega)*cos(kappa);R[0][2]=-sin(phi)*cos(omega);R[1][0]=cos(omega)*sin(kappa);R[1][1]=cos(omega)*cos(kappa);R[1][2]=-sin(omega);R[2][0]=sin(phi)*cos(kappa)+cos(phi)*sin(omega)*sin(kappa);R[2][1]=-sin(phi)*sin(kappa) + cos(phi)*sin(omega)*cos(kappa);R[2][2]=cos(phi)*cos(omega);//将未知数的近似值代人共线方程式，计算(x),(y)for(i=0;i<n;i++){(x)[i]=-f*(R[0][0]*(coor[i].Xt-Xs)+R[1][0]*(coor[i].Yt-Ys)+R[2][0]*(coor[i].Zt-Zs))/(R[0][2]*(coor[ i].Xt-Xs)+R[1][2]*(coor[i].Yt-Ys)+R[2][2]*(coor[i].Zt-Zs));(y)[i]=-f*(R[0][1]*(coor[i].Xt-Xs)+R[1][1]*(coor[i].Yt-Ys)+R[2][1]*(coor[i].Zt-Zs))/(R[0][2]*(coor[i].Xt -Xs)+R[1][2]*(coor[i].Yt-Ys)+R[2][2]*(coor[i].Zt-Zs));}//计算矩阵A的各个元素for(i=0;i<n;i++){A[2*i][0]=(R[0][0]*f + R[0][2] * (coor[i].x - xo)) / (R[0][2]*(coor[i].Xt - Xs) + R[1][2]*(coor[i].Yt - Ys) + R[2][2]*(coor[i].Zt - Zs));A[2*i][1]=(R[1][0]*f + R[1][2] * (coor[i].x - xo)) / (R[0][2]*(coor[i].Xt - Xs) + R[1][2]*(coor[i].Yt - Ys) + R[2][2]*(coor[i].Zt - Zs));A[2*i][2]=(R[2][0]*f + R[2][2] * (coor[i].x - xo)) / (R[0][2]*(coor[i].Xt - Xs) + R[1][2]*(coor[i].Yt - Ys) + R[2][2]*(coor[i].Zt - Zs));A[2*i][3]=(coor[i].y-yo)*sin(omega)-((coor[i].x-xo)/f*((coor[i].x-xo)*cos(kappa)-(coor[i].y-yo)* sin(kappa))+f*cos(kappa))*cos(omega);A[2*i][4]=-f*sin(kappa)-(coor[i].x-xo)/f*((coor[i].x-xo)*sin(kappa)+(coor[i].y-yo)*cos(kappa));A[2*i][5]=coor[i].y-yo;A[2*i+1][0]=(R[0][1]*f + R[0][2] * (coor[i].y - yo)) / (R[0][2]*(coor[i].Xt - Xs) + R[1][2]*(coor[i].Yt - Ys) + R[2][2]*(coor[i].Zt - Zs));A[2*i+1][1]=(R[1][1]*f + R[1][2] * (coor[i].y - yo)) / (R[0][2]*(coor[i].Xt - Xs) + R[1][2]*(coor[i].Yt - Ys) + R[2][2]*(coor[i].Zt - Zs));A[2*i+1][2]=(R[2][1]*f + R[2][2] * (coor[i].y - yo)) / (R[0][2]*(coor[i].Xt - Xs) + R[1][2]*(coor[i].Yt - Ys) + R[2][2]*(coor[i].Zt - Zs));A[2*i+1][3]=-(coor[i].x-xo)*sin(omega)-((coor[i].y-yo)/f*((coor[i].x-xo)*cos(kappa)-(coor[i].y-y o)*sin(kappa))-f*sin(kappa))*cos(omega);A[2*i+1][4]=-f*cos(kappa)-(coor[i].y-yo)/f*((coor[i].x-xo)*sin(kappa)+(coor[i].y-yo)*cos(kappa ));A[2*i+1][5]=-(coor[i].x-xo);}for (i=0;i<6;i++){for (j=0;j<6;j++){printf("%f ",A[i][j]);if (j%5==0&&j!=0){printf("\n");}}}printf("'\n");//计算ATA的各个元素for(i=0;i<6;i++)for(j=0;j<6;j++)ATA[i][j]=A[0][i]*A[0][j]+A[1][i]*A[1][j]+A[2][i]*A[2][j]+A[3][i]*A[3][j]+A[4][i]*A[4][j]+A[5][i]*A[5][j]+A[6][i]*A[6][j]+A[7][i]*A[7][j];for (i=0;i<6;i++){for (j=0;j<6;j++){_ATA[i][j]=ATA[i][j];printf("%f ",ATA[i][j]);if (j%5==0&&j!=0){printf("\n");}}}printf("\n");//计算矩阵L的各个元素for(i=0;i<n;i++){L[2*i][0]=coor[i].x+f*(R[0][0]*(coor[i].Xt-Xs)+R[1][0]*(coor[i].Yt-Ys)+R[2][0]*(coor[i].Zt-Zs))/(R [0][2]*(coor[i].Xt - Xs) + R[1][2]*(coor[i].Yt - Ys) + R[2][2]*(coor[i].Zt - Zs));L[2*i+1][0]=coor[i].y+f*(R[0][1]*(coor[i].Xt-Xs)+R[1][1]*(coor[i].Yt-Ys)+R[2][1]*(coor[i].Zt-Zs))/(R[0 ][2]*(coor[i].Xt - Xs) + R[1][2]*(coor[i].Yt - Ys) + R[2][2]*(coor[i].Zt - Zs));}//计算矩阵ATL的各个元素值for(i=0;i<6;i++)ATL[i][0]=A[0][i]*L[0][0]+A[1][i]*L[1][0]+A[2][i]*L[2][0]+A[3][i]*L[3][0]+A[4][i]*L[4][0]+A[5][i]*L[5][0]+A[6][i]*L[6][0]+A[7][i]*L[7][0];/* for (i=0;i<6;i++){for (j=0;j<1;j++){ATL[i][j]=0;for (int k=0;k<8;k++){ATL[i][j]+=A[k][i]*L[k][j];}}}*///调用函数计算矩阵ATA的逆矩阵// inverse(ATA);ContraryMatrix(*_ATA, *ATA, 6);for (i=0;i<6;i++){for (j=0;j<6;j++){printf("%f ",ATA[i][j]);if (j%5==0&&j!=0){printf("\n");}}}//计算矩阵X的各个元素值for(i=0;i<6;i++){X[i][0]=ATA[i][0]*ATL[0][0]+ATA[i][1]*ATL[1][0]+ATA[i][2]*ATL[2][0]+ATA[i][3]*ATL[3][0] +ATA[i][4]*ATL[4][0]+ATA[i][5]*ATL[5][0];printf("%f ",X[i][0]);}// for (i=0;i<6;i++)// {// for (j=0;j<1;j++)// {// X[i][j]=0;// for (int k=0;k<6;k++)// {// X[i][j]+=ATA[i][k]*ATL[k][j];//// }// printf("%f",X[i][j]);// }// }DXs=X[0][0];DYs=X[1][0];DZs=X[2][0];Dphi=X[3][0];Domega=X[4][0];Dkappa=X[5][0];Xs+=DXs;Ys+=DYs;Zs+=DZs;phi+=Dphi;omega+=Domega;kappa+=Dkappa;//计算矩阵V的各个元素for(i=0;i<n;i++){V[2*i][0]=A[2*i][0]*DXs+A[2*i][1]*DYs+A[2*i][2]*DZs+A[2*i][3]*Dphi+A[2*i][4]*Domega+A[ 2*i][5]*Dkappa-L[2*i][0];V[2*i+1][0]=A[2*i+1][0]*DXs+A[2*i+1][1]*DYs+A[2*i+1][2]*DZs+A[2*i+1][3]*Dphi+A[2*i+1][ 4]*Domega+A[2*i+1][5]*Dkappa-L[2*i+1][0];}/* //像点坐标改正后的值for(i=0;i<n;i++){coor[i].x+=V[2*i][0];coor[i].y+=V[2*i+1][0];}*/}//判断限差的条件while(!(fabs(Dphi) < 0.1/60/180*PI && fabs(Domega) < 0.1/60/180*PI && fabs(Dkappa) < 0.1/60/180*PI));//外方位元素计算完毕//有关精度的计算for(i=0;i<6;i++)Q[i][i]=ATA[i][i];m0=sqrt((V[0][0]*V[0][0]+V[1][0]*V[1][0]+V[2][0]*V[2][0]+V[3][0]*V[3][0]+V[4][0]*V[4][0]+V [5][0]*V[5][0]+V[6][0]*V[6][0]+V[7][0]*V[7][0])/(2*n-6));//计算各未知数的中误差m1=sqrt(Q[0][0])*m0;m2=sqrt(Q[1][1])*m0;m3=sqrt(Q[2][2])*m0;m4=sqrt(Q[3][3])*m0;m5=sqrt(Q[4][4])*m0;m6=sqrt(Q[5][5])*m0;printf("改正后的相对坐标及其对应的地面点坐标（x,y,Xt,Yt,Zt）:\n");for(i=0;i<n;i++){printf("%lf\t%lf\t%lf\t%lf\t%lf",(x)[i],(y)[i],coor[i].Xt,coor[i].Yt,coor[i].Zt);printf("\n");}printf("旋转矩阵R：\n");for(i=0;i<3;i++){for(j=0;j<3;j++)printf("%lf\t",R[i][j]);printf("\n");}printf("外方位元素(Xs,Ys,Zs,phi,omega,kappa)及改正值：\n");printf("%lf+-%lf\n%lf+-%lf\n%lf+-%lf\n%lf+-%lf\n%lf+-%lf\n%lf+-%lf\n",Xs,DXs,Ys,DYs,Zs,DZs, phi,Dphi,omega,Domega,kappa,Dkappa);printf("单位权中误差:\n");printf("%0.9f\n",m0);printf("外方位元素(Xs,Ys,Zs,phi,omega,kappa)的中误差:\n");printf("%lf\n%lf\n%lf\n%lf\n%lf\n%lf\n",m1,m2,m3,m4,m5,m6);//计算完毕，输出结果,以文件方式保存printf("计算结果存储在文件中\n");FILE *fp;fp=fopen("计算结果.txt","w");fprintf(fp,"改正后的相对坐标及其对应的地面点坐标（x,y,Xt,Yt,Zt）:\n");for(i=0;i<n;i++){fprintf(fp,"%lf\t%lf\t%lf\t%lf\t%lf",coor[i].x,coor[i].y,coor[i].Xt,coor[i].Yt,coor[i].Zt);fprintf(fp,"\n");}fprintf(fp,"旋转矩阵R：\n");for(i=0;i<3;i++){for(j=0;j<3;j++)fprintf(fp,"%lf\t",R[i][j]);fprintf(fp,"\n");}fprintf(fp,"外方位元素(Xs,Ys,Zs,phi,omega,kappa)及改正值：\n");fprintf(fp,"%lf+-%lf\n%lf+-%lf\n%lf+-%lf\n%lf+-%lf\n%lf+-%lf\n%lf+-%lf\n",Xs,DXs,Ys,DYs,Zs,D Zs,phi,Dphi,omega,Domega,kappa,Dkappa);fprintf(fp,"单位权中误差:\n");fprintf(fp,"%0.9f\n",m0);fprintf(fp,"外方位元素(Xs,Ys,Zs,phi,omega,kappa)的中误差:\n");fprintf(fp,"%lf\n%lf\n%lf\n%lf\n%lf\n%lf\n",m1,m2,m3,m4,m5,m6);fclose(fp);}。

后方交会代码

后⽅交会代码///<summary>///测点信息///</summary>///<param name="list"></param>///<returns></returns> public class PointInfo{///<summary>///是否基准点 0基准点 1 测量点 2全站仪///</summary>public int Type { get; set; }///<summary>/// X坐标///</summary>public double Y { get; set; }///<summary>/// Y坐标///</summary>public double X { get; set; }///<summary>///⽔平位置///</summary>public double Hz { get; set; }}///<summary>///后⽅交会算法返回坐标///</summary>///<param name="list"></param>///<returns></returns>private static double[] Resection(List<PointInfo> list){var dataList = list.Where(x => x.Type == 0).ToList();if (dataList.Count < 3) return null;//⾄少三个基准点var pointA = dataList[0];var pointB = dataList[1];var pointC = dataList[2];//三个基准点坐标double xa = pointA.X;double xb = pointB.X;double xc = pointC.X;double ya = pointA.Y;double yb = pointB.Y;double yc = pointC.Y;double alpha = pointC.Hz - pointB.Hz;double beta = pointA.Hz - pointC.Hz;double gamma = pointB.Hz - pointA.Hz;double cotA = ((xb - xa) * (xc - xa) + (yb - ya) * (yc - ya)) / ((xb - xa) * (yc - ya) - (yb - ya) * (xc - xa));double cotB = ((xc - xb) * (xa - xb) + (yc - yb) * (ya - yb)) / ((xc - xb) * (ya - yb) - (yc - yb) * (xa - xb));double cotC = ((xa - xc) * (xb - xc) + (ya - yc) * (yb - yc)) / ((xa - xc) * (yb - yc) - (ya - yc) * (xb - xc));double pA = 1.0 / (cotA - Math.Pow(Math.Tan(alpha), -1));double pB = 1.0 / (cotB - Math.Pow(Math.Tan(beta), -1));double pC = 1.0 / (cotC - Math.Pow(Math.Tan(gamma), -1));double xp = (pA * xa + pB * xb + pC * xc) / (pA + pB + pC);double yp = (pA * ya + pB * yb + pC * yc) / (pA + pB + pC);return new double[2] { xp, yp };}。

基于C#混合编程的角锥体法空间后方交会实现

空间后方交会是摄影测量学中一个很重要的概念,它是利用控制点坐标来求解摄影测量学中外方位元素的过程。

当解算出外方位元素后,可根据相片的外方位元素值,确定摄影瞬间相片在物方坐标系下的位置和姿态[1]。

本文采用计算机语言,利用角锥体法空间后方交会法实现求解外方位元素,并通过算例证明了程序的正解性。

1角锥体法原理基于角锥体法的空间后方交会是一种解算外方位元素的有效方法,它是根据摄影像方与物方空间所对应光束线之间夹角相等的原理,先后分别解算出外方位线元素和外方位角元素,由于是分别解算线元素与角元素,从而有效克服了共线方程法解算线元素与角元素存在的相关性。

图1为角锥体模型,A 、B 、C 为3个已知地面控制点,a ,b ,c 为相应的像点,S 为摄影中心。

从模型可以看出,根据3个像点坐标,相机焦距f ,在三角形aSb 与三角形ASB 中,∠aSb =∠A SB 。

同理,在三角形aSc 与三角形ASC ,三角形bSc 与三角形BSC 中,∠aSc =∠A SC ,∠bSc =∠BSC 。

在三角形ASB 中,A ,B 坐标已知,∠ASB 已求得,根据余弦公式可列出余弦方程。

同理在三角形ASC 和BSC 中也可列出2个余弦方程,3个方程组成一个含3个方程3个距离未知数的非线性方程组,由此方程组可解出投影中心到3个地面控制点的距离,由3个距离,3个控制点坐标,进而交会解出投影中心的坐标Xs ,Ys ,Zs 。

由投影中心坐标,3个地面控制点坐标,投影中心与像点和相应地面点的距离之比,根据共线方程,求解旋转矩阵R 。

根据R 中的元素进而求得外方位角元素。

2混合编程具体实现2.1matlab M 文件编写及相应组件DLL 生成打开运行matlab,建立一个M 文件yuxuan.m 。

文件用fanction 函数格式语法写入,fanction 函数格式为fanction +输出函数名=输入函数名(输入参数1,输入参数2,输入参数3,…,输入参数n )[2]。

后方交会程序实现(c语言版)

#include<stdio.h>#include<stdlib.h>#include<malloc.h>#include<math.h>#include<conio.h>#define N 4#define T 1.41421void turn(double *A,double A2[],int m,int n) //计算矩阵的转置{ int i,j;for(i=0;i<m;i++)for(j=0;j<n;j++)A2[j*m+i]=A[i*n+j];}void mulAB(double *A,double *B,double *C,int am,int an,int bm,int bn) //计算两矩阵相乘{int i,j,l,u;if(an!=bm){printf("error!cannot do the multiplication.\n");return;}for(i=0;i<am;i++)for(j=0;j<bn;j++){u=i*bn+j;C[u]=0.0;for(l=0;l<an;l++)C[u]+=A[i*an+l]*B[l*bn+j];}return;}double *inv(double *a,int n) //计算矩阵的逆,本程序的难点,采用高斯约旦全选主元法{int *is,*js,i,j,k,l,u,v;double d,p;is=(int*)malloc(n*sizeof(int)); js=(int*)malloc(n*sizeof(int)); for (k=0; k<=n-1; k++){d=0.0;for (i=k;i<n;i++)for (j=k;j<n;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("error not inv\n"); return NULL;}if (is[k]!=k)for (j=0;j<n;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;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;j++) if (j!=k){u=k*n+j;a[u]=a[u]*a[l]; }for (i=0;i<n;i++) if (i!=k)for (j=0;j<n;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;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;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 a;}main()//主函数,空间后方交会的计算{ FILE *fp; //定义一个文件指针long m;int i,j=0,it;double G[1000];doublef,t,w,k,limit=0,S1=0.0,S2=0.0,S3=0.0,x[N]={0},y[N]={0},x0[N]={0},y0 [N]={0},X[N]={0},Y[N] ={0},Z[N]={0},Xs0,Ys0,Zs0;doublea[3],b[3],c[3],A[2*N*6],AT[6*2*N],ATA[6*6],*ATA_=NULL,l[2*N],ATl[6] ,V[6]={0};double F[6],Qx[6][6],Mi[6][6];if((fp=fopen("e:\\shuju.txt","r"))==NULL) //使fp指向被打开的shuju.txt 文件{ //并判断文件是否打开成功printf("\nerror on open shuju.txt\n");getch();exit(1);}for(i=0;i<N;i++)fscanf(fp,"%lf%lf%lf%lf%lf",&x[i],&y[i],&X[i],&Y[i],&Z[i]); //将文件中的数据赋给主函数定义的变量printf("原始数据:\n");printf("x\t\ty\t\tX\t\tY\t\tZ\t\t\n");//输出文件中的原始数据for(i=0;i<N;i++)printf("%.3lf\t\t%.3lf\t\t%.3lf\t%.3lf\t%.3lf\n",x[i],y[i],X[i],Y[i ],Z[i]);printf("\n请输入摄影机主距和摄影比例尺分母;f,m:");scanf("%lf%ld",&f,&m); //输入f,mf=f/1000.0;for(i=0;i<N;i++){ x[i]/=1000.0;y[i]/=1000.0;S1+=X[i];S2+=Y[i];S3+=Z[i];}Xs0=S1/N;Ys0=S2/N;Zs0=m*f+S3/N;t=0.0;w=0.0;k=0.0; //计算外方位元素的初始值while(1){printf("\n---------------------------------第%d次计算------------------------------\n",j+1);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); //计算旋转矩阵for(i=0;i<N;i++){x0[i]=-f*(a[0]*(X[i]-Xs0)+b[0]*(Y[i]-Ys0)+c[0]*(Z[i]-Zs0))/(a[2]*(X[i]-Xs0)+b[2]*(Y[i]-Ys0)+c [2]*(Z[i]-Zs0));y0[i]=-f*(a[1]*(X[i]-Xs0)+b[1]*(Y[i]-Ys0)+c[1]*(Z[i]-Zs0))/(a[2]*(X[i]-Xs0)+b[2]*(Y[i]-Ys0)+c [2]*(Z[i]-Zs0));//计算像点坐标近似值G[i]=a[2]*(X[i]-Xs0)+b[2]*(Y[i]-Ys0)+c[2]*(Z[i]-Zs0);}for(i=0;i<N;i++){A[i*12+0]=(a[0]*f+a[2]*x[i])/G[i];A[i*12+1]=(b[0]*f+b[2]*x[i])/G[i];A[i*12+2]=(c[0]*f+c[2]*x[i])/G[i];A[i*12+3]=y[i]*sin(w)-(x[i]*(x[i]*cos(k)-y[i]*sin(k))/f+f*cos(k))*cos(w);A[i*12+4]=-f*sin(k)-x[i]*(x[i]*sin(k)+y[i]*cos(k))/f;A[i*12+5]=y[i];A[i*12+6]=(a[1]*f+a[2]*y[i])/G[i];A[i*12+7]=(b[1]*f+b[2]*y[i])/G[i];A[i*12+8]=(c[1]*f+c[2]*y[i])/G[i];A[i*12+9]=-x[i]*sin(w)-(y[i]*(x[i]*cos(k)-y[i]*sin(k))/f-f*sin(k))*cos(w);A[i*12+10]=-f*cos(k)-y[i]*(x[i]*sin(k)+y[i]*cos(k))/f;A[i*12+11]=-x[i];l[i*2+0]=x0[i]-x[i];l[i*2+1]=y0[i]-y[i]; //计算误差方程的系数阵以及lx,ly}// printf("output matrix: A\n");// printmatrix(A,2*N,6);// printf("output matrix: l\n");// printmatrix(l,2*N,1);turn(A,A T,2*N,6); //计算AT// printf("output matrix: AT\n");// printmatrix(AT,6,2*N);mulAB(AT,A,ATA,6,2*N,2*N,6); //计算ATA,组法方程ATA_=inv(ATA,6); //计算ATA的逆,中间量int p;int cnt=-1;for(it=0;it<36;it++){p=it%6;if(it%6==0){cnt++;}Qx[cnt][p++]=A TA_[it];}for(int it=0;it<6;it++){for(int jt=0;jt<6;jt++){if(it!=jt){Qx[it][jt]=0;//提取Qx的主对角线元素=Qii }// printf("%-10.3lf ",Qx[it][jt]);}// printf("\n");}mulAB(A T,l,ATl,6,2*N,2*N,1); //计算常数项ATL// printf("outpit matrinx: ATl\n");// printmatrix(ATl,6,1);mulAB(A TA_,ATl,V,6,6,6,1); //解法方程,求改正数, // printf("output matrix: V\n");// printmatrix(V,6,1);Xs0+=V[0];Ys0+=V[1];Zs0+=V[2];t+=V[3];w+=V[4];k+=V[5];for(i=0;i<6;i++){F[i]=V[i]/T;//m0}printf("第%d次计算的外方位元素为:\n",++j);printf("Xs=%.5lf\tYs=%.5lf\tZs=%.5lf\nt=%.5lf\tw=%.5lf\tk=%.5lf\n", Xs0,Ys0,Zs0,t,w,k);if(Xs0-limit<=0.0001&&Xs0-limit>=-0.0001){ //控制迭代次数break;}limit=Xs0;}printf("\n外方位元素为:\n");printf("Xs=%.5lf\tYs=%.5lf\tZs=%.5lf\nt=%.5lf\tw=%.5lf\tk=%.5lf\n", Xs0,Ys0,Zs0,t,w,k);printf("\n中误差为:\n");for(i=0;i<6;i++){for(j=0;j<6;j++){Mi[i][j]=F[i]*(sqrt(Qx[i][j]));//mi=m0*Qii开根号printf("%-13.10lf",Mi[i][j]);}printf("\n");}fclose(fp); return 0; }。

空间后方交会C++程序代码

摄影测量后方交会程序（c/c++）输入数据截图：结果截图：程序源代码（其中的矩阵求逆在前面已经有了，链接）：#include <stdio.h>#include <stdlib.h>#include <math.h>const double PRECISION=1e-5;typedef double DOUBLE[5];int InputData(int &Num, DOUBLE *&Data,double &m,double &f);int Resection(const int &Num,const DOUBLE *&Data,const double &m,const double &f);int InverseMatrix(double *matrix,const int &row);int main(int argc, char* argv[]){DOUBLE *Data=NULL;int Num;double f(0),m(0);if(InputData(Num,Data,m,f)){if (Data!=NULL){delete []Data;}return 1;}if(Resection(Num,Data,m,f)){if (Data!=NULL){delete []Data;}return 1;}if (Data!=NULL){delete []Data;}printf("解算完毕...\n");do{printf("计算结果保存于\"结果.txt\"文件中\n""请选择操作（输入P打开结果数据，R打开原始数据，其它退出程序）：");fflush(stdin); //刷新输入流char order=getchar();if ('P'==order || 'p'==order){system("结果.txt");}else if ('R'==order || 'r'==order){system("data.txt");}elsebreak;system("cls");}while(1);system("PAUSE");return 0;}/***********************************************函数名：InputData*函数介绍：从文件（data.txt）中读取数据,*文件格式如下：*点数 m(未知写作0)* 内方位元素(f x0 y0)*编号 x y X Y Z*下面是一个实例：4 0153.24 0 01 -86.15 -68.99 36589.41 25273.32 2195.172 -53.40 82.21 37631.08 31324.51 728.693 -14.78 -76.63 39100.97 24934.98 2386.504 10.46 64.43 40426.54 30319.81 757.31*参数：(in/out)Num(点数)，*(in/out)Data(存放数据),m,f,x0,y0*返回值：int ，0成功，1文件打开失败，2控制点个*数不足，3文件格式错误*作者：vcrs*完成时间：09-10-4**********************************************/int InputData(int &Num, DOUBLE *&Data,double &m,double &f) {double x0,y0;FILE *fp_input;if (!(fp_input=fopen("data.txt","r"))){return 1;}fscanf(fp_input,"%d%lf",&Num,&m);if (Num<4){return 2;}fscanf(fp_input,"%lf%lf%lf",&f,&x0,&y0);f/=1000;if (m<0 || f<0){return 3;}Data=new DOUBLE[Num];double *temp= new double[Num-1];double scale=0;int i;for (i=0;i<Num;i++){//读取数据，忽略编号if(fscanf(fp_input,"%*d%lf%lf%lf%lf%lf",&Data[i][0],&Data[i][1],&Data[i][2],&Data[i][3],&Data[i][4])!=5){return 3;}//单位换算成mData[i][0]/=1000.0;Data[i][1]/=1000.0;}//如果m未知则归算其值if (0==m){for (i=0;i<Num-1;i++){temp[i]=(Data[i][2]-Data[i+1][2])/(Data[i][0]-Data[i+1][0])+ (Data[i][3]-Data[i+1][3])/(Data[i][1]-Data[i+1][1]);scale+=temp[i]/2.0;}m=scale/(Num-1);}fclose(fp_input);delete []temp;return 0;}/***********************************************函数名：MatrixMul*函数介绍：求两个矩阵的积,*参数：Jz1(第一个矩阵),row(第一个矩阵行数),*Jz2(第二个矩阵),row(第二个矩阵列数),com(第一个*矩阵列数),(out)JgJz(存放结果矩阵)*返回值：void*作者：vcrs*完成时间：09-10-4**********************************************/void MatrixMul(double *Jz1,const int &row,double *Jz2,const int &line,const int &com,double *JgJz){for (int i=0;i<row;i++){for (int j=0;j<line;j++){double temp=0;for (int k=0;k<com;k++){temp+=*(Jz1+i*com+k)*(*(Jz2+k*line+j));}*(JgJz+i*line+j)=temp;}}}/***********************************************函数名：OutPut*函数介绍：向结果.txt文件输出数据*参数：Q协因数阵，m精度，m0单位权中误差，6个外*方位元素，旋转矩阵*返回值：int，0成功，1失败*作者：vcrs*完成时间：09-10-4**********************************************/int OutPut(const double *&Q,const double *&m,const double &m0, const double &Xs,const double &Ys,const double &Zs,const double &Phi,const double &Omega,const double &Kappa,const double *R){FILE *fp_out;if (!(fp_out=fopen("结果.txt","w"))){return 1;}FILE *fp_input;if (!(fp_input=fopen("data.txt","r"))){return 1;}fprintf(fp_out,"**************************************""**************************************""**************************************""*********************************\n");fprintf(fp_out,"\n空间后方交会程序(C\\C++)\n遥感信息工程学院\n班级：""00000\n学号：0000000\n姓名：vcrs\n\n");fprintf(fp_out,"**************************************""**************************************""**************************************""*********************************\n");fprintf(fp_out,"已知数据：\n\n已知点数：");int num;double temp,x,y;fscanf(fp_input,"%d%lf",&num,&temp);fprintf(fp_out,"%d\n",num);fprintf(fp_out,"摄影比例尺（0表示其值位置）：");fprintf(fp_out,"%10.0lf\n",temp);fprintf(fp_out,"内方位元素(f x0 y0)：");fscanf(fp_input,"%lf%lf%lf",&temp,&x,&y);fprintf(fp_out,"%10lf\t%10lf\t%10lf\n",temp,x,y);for (int i=0;i<num;i++){double temp[5];fscanf(fp_input,"%*d%lf%lf%lf%lf%lf",&temp[0],&temp[1],&temp[2],&temp[3],&temp[4]);fprintf(fp_out,"%3d\t%10lf\t%10lf\t%10lf\t%10lf\t%10lf\n", i+1,temp[0],temp[1],temp[2],temp[3],temp[4]);}fclose(fp_input);fprintf(fp_out,"**************************************""**************************************""*********************************\n");fprintf(fp_out,"计算结果如下：\n\n外方位元素：\n"); fprintf(fp_out,"\tXs=%10lf\n",Xs);fprintf(fp_out,"\tYs=%10lf\n",Ys);fprintf(fp_out,"\tZs=%10lf\n",Zs);fprintf(fp_out,"\tPhi=%10lf\n",Phi);fprintf(fp_out,"\tOmega=%10lf\n",Omega);fprintf(fp_out,"\tKappa=%10lf\n\n",Kappa);fprintf(fp_out,"旋转矩阵：\n");for (i=0;i<3;i++){fprintf(fp_out,"\t");for (int j=0;j<3;j++){fprintf(fp_out,"%10lf\t",*(R+i*3+j));}fprintf(fp_out,"\n");}fprintf(fp_out,"\n单位权中误差：%10lf\n\n",m0);fprintf(fp_out,"协因数阵：\n");for (i=0;i<6;i++){fprintf(fp_out,"\t");for (int j=0;j<6;j++){fprintf(fp_out,"%20lf\t",*(Q+i*6+j));}fprintf(fp_out,"\n");}fprintf(fp_out,"\n外方位元素精度：");for (i=0;i<6;i++){fprintf(fp_out,"%10lf\t",m[i]);}fprintf(fp_out,"\n");fprintf(fp_out,"**************************************""**************************************""*********************************\n");fclose(fp_out);return 0;}/***********************************************函数名：Resection*函数介绍：计算*参数：Num(点数)，Data(数据)，m,,f(焦距),x0,y0*返回值：int，0成功，其它失败*作者：vcrs*完成时间：09-10-4**********************************************/int Resection(const int &Num,const DOUBLE *&Data,const double &m, const double &f){double Xs=0,Ys=0,Zs=0;int i,j;//设置初始值for (i=0;i<Num;i++){Xs+=Data[i][2];Ys+=Data[i][3];}Xs/=Num;Ys/=Num;Zs=m*f;double Phi(0),Omega(0),Kappa(0);double R[3][3]={0.0};double *L=new double[2*Num];typedef double Double6[6];Double6 *A=new Double6[2*Num];double *AT=new double[2*Num*6];double *ATA=new double[6*6];double *ATL=new double[6];double *Xg=new double[6];//迭代计算do{//旋转矩阵R[0][0]=cos(Phi)*cos(Kappa)-sin(Phi)*sin(Omega)*sin(Kappa);R[0][1]=-cos(Phi)*sin(Kappa)-sin(Phi)*sin(Omega)*cos(Kappa);R[0][2]=-sin(Phi)*cos(Omega);R[1][0]=cos(Omega)*sin(Kappa);R[1][1]=cos(Omega)*cos(Kappa);R[1][2]=-sin(Omega);R[2][0]=sin(Phi)*cos(Kappa)+cos(Phi)*sin(Omega)*sin(Kappa);R[2][1]=-sin(Phi)*sin(Kappa)+cos(Phi)*sin(Omega)*cos(Kappa);R[2][2]=cos(Phi)*cos(Omega);for (i=0;i<Num;i++){double X=R[0][0]*(Data[i][2]-Xs)+R[1][0]*(Data[i][3]-Ys)+ R[2][0]*(Data[i][4]-Zs);double Y=R[0][1]*(Data[i][2]-Xs)+R[1][1]*(Data[i][3]-Ys)+ R[2][1]*(Data[i][4]-Zs);double Z=R[0][2]*(Data[i][2]-Xs)+R[1][2]*(Data[i][3]-Ys)+ R[2][2]*(Data[i][4]-Zs);double xxx,yyy;xxx=-f*X/Z;yyy=-f*Y/Z;//常数项L[2*i]=Data[i][0]-(-f*X/Z);L[2*i+1]=Data[i][1]-(-f*Y/Z);A[2*i][0]=(R[0][0]*f+R[0][2]*(xxx))/Z;A[2*i][1]=(R[1][0]*f+R[1][2]*(xxx))/Z;A[2*i][2]=(R[2][0]*f+R[2][2]*(xxx))/Z;A[2*i][3]=(yyy)*sin(Omega)-(((xxx)/f)*((xxx)*cos(Kappa)-(yyy)*sin(Kappa))+f*cos(Kappa))*cos(Omega);A[2*i][4]=-f*sin(Kappa)-((xxx)/f)*((xxx)*sin(Kappa)+(yyy)*cos(Kappa));A[2*i][5]=(yyy);A[2*i+1][0]=(R[0][1]*f+R[0][2]*(yyy))/Z;A[2*i+1][1]=(R[1][1]*f+R[1][2]*(yyy))/Z;A[2*i+1][2]=(R[2][1]*f+R[2][2]*(yyy))/Z;A[2*i+1][3]=-(xxx)*sin(Omega)-(((yyy)/f)*((xxx)*cos(Kappa)-(yyy)*sin(Kappa))-f*sin(Kappa))*cos(Omega);A[2*i+1][4]=-f*cos(Kappa)-((yyy)/f)*((xxx)*sin(Kappa)+(yyy)*cos(Kappa));A[2*i+1][5]=-(xxx);}//求矩阵A的转置矩阵ATfor (i=0;i<2*Num;i++){for (j=0;j<6;j++){*(AT+j*2*Num+i)=A[i][j];}}//求ATAMatrixMul(AT,6,&A[0][0],6,2*Num,ATA);if(InverseMatrix(ATA,6))return 1;MatrixMul(AT,6,L,1,2*Num,ATL);MatrixMul(ATA,6,ATL,1,6,Xg);Xs+=Xg[0];Ys+=Xg[1];Zs+=Xg[2];Phi+=Xg[3];Omega+=Xg[4];Kappa+=Xg[5];} while(fabs(Xg[0])>=PRECISION ||fabs(Xg[1])>=PRECISION || fabs(Xg[2])>=PRECISION ||fabs(Xg[3])>=PRECISION ||fabs(Xg[4])>=PRECISION || (Xg[5])>=PRECISION);//注：协因数阵，旋转矩阵等计算本应该使用最后外方位元素值，//由于变换很小忽略double *Q=ATA;double *V=new double[2*Num];MatrixMul(&A[0][0],2*Num,Xg,1,6,V);double VTV=0;for(i=0;i<2*Num;i++){V[i]-=L[i];VTV+=V[i]*V[i];}double m0=sqrt(VTV/(2*Num-6));double *mm=new double[6];for (i=0;i<6;i++){mm[i]=sqrt(*(Q+i*6+i))*m0;}OutPut(Q,mm,m0,Xs,Ys,Zs,Phi,Omega,Kappa,&R[0][0]);delete []L;delete []A;delete []AT;delete []ATA;delete []ATL;delete []Xg;delete []mm;delete []V;return 0;}void swap(double &a,double &b){double temp=a;a=b;b=temp;}/***********************************************函数名：InverseMatrix*函数介绍：求矩阵的逆（高斯-约当法）*输入参数：(in/out)matrix（矩阵首地址），*（in）row（矩阵阶数）*输出参数：matrix（原矩阵的逆矩阵）*返回值：int ，0成功，1失败*调用函数：swap(double&,double&)*作者：vcrs*完成时间：09-10-4**********************************************/int InverseMatrix(double *matrix,const int &row) {double *m=new double[row*row];double *ptemp,*pt=m;int i,j;ptemp=matrix;for (i=0;i<row;i++)for (j=0;j<row;j++){*pt=*ptemp;ptemp++;pt++;}}int k;int *is=new int[row],*js=new int[row];for (k=0;k<row;k++){double max=0;//全选主元//寻找最大元素for (i=k;i<row;i++){for (j=k;j<row;j++){if (fabs(*(m+i*row+j))>max){max=*(m+i*row+j);is[k]=i;js[k]=j;}}}if (0 == max){return 1;}//行交换if (is[k]!=k){for (i=0;i<row;i++){swap(*(m+k*row+i),*(m+is[k]*row+i));}}//列交换if (js[k]!=k){for (i=0;i<row;i++){swap(*(m+i*row+k),*(m+i*row+js[k]));}}*(m+k*row+k)=1/(*(m+k*row+k));for (j=0;j<row;j++){if (j!=k){*(m+k*row+j)*=*((m+k*row+k));}}for (i=0;i<row;i++){if (i!=k){for (j=0;j<row;j++){if(j!=k){*(m+i*row+j)-=*(m+i*row+k)**(m+k*row+j);}}}}for (i=0;i<row;i++){if(i!=k){*(m+i*row+k)*=-(*(m+k*row+k));}}}int r;//恢复行列for (r=row-1;r>=0;r--){if (js[r]!=r){for (j=0;j<row;j++){swap(*(m+r*row+j),*(m+js[r]*row+j));}}if (is[r]!=r){for (i=0;i<row;i++){swap(*(m+i*row+r),*(m+i*row+is[r]));}}}ptemp=matrix;pt=m;for (i=0;i<row;i++){for (j=0;j<row;j++){*ptemp=*pt;ptemp++;pt++;}}delete []is;delete []js;delete []m;return 0;}。

后方交会测量原理及其程序实现

后方交会测量原理及其程序实现摘要本文针对后方交会测量原理的定义、意义进行了概述，并归纳总结了后方交会测量原理的3种常见情况。

然后利用Visual Basic 6.0软件作为平台，对这3种情况进行了程序实现。

关键词后方交会；余弦定理；测边交会；边角交会；Visual Basic 6.00 引言随着计算机的广泛普及和各种计算机语言的不断完善，利用计算机语言编制一些测量中常用的程序，以此来解决实地工作中的计算问题，深获测绘工作者的青睐。

本文总结了后方交会测量原理的几种情况，并尝试利用Visual Basic 6.0编写了后方交会测量原理的计算程序，在实地工作中得到了很好的应用。

1 后方交会测量原理1.1 余切公式根据A、B、C三个已知点的坐标（XA，YA）、（X B ，YB）、（XC，YC）及两个水平观测角α、β，计算出已知点到未知点的坐标方位角（这里以αBP 为例）。

然后根据这个坐标方位角及已知坐标增量和α、β角求出已知点到未知点的坐标增量，最后求得未知点的坐标XP，YP。

由图１可知，CP、BP、AP三条直线的方程式为式中：为B点至P点的坐标方位角；为A点坐标；为B点坐标；为C点坐标；为P点坐标。

类似的计算方法还有赫尔默特公式（重心公式）等，这里就不再详述。

在此特别强调后方交会危险圆的问题：如图2所示，当未知点P在3个已知点所在的圆上移动时，总有下列式子成立，即α和β固定不变，它说明此时仅有α、β这两个角和3个已知点不能唯一确定P点位置，因此称3个已知点所在的圆为后方交会危险圆。

当P点位于危险圆时，无论用何种后方交会公式，均无法求出P点坐标。

实际作业中，即使P点不正好在危险圆上，而是接近危险圆时，计算结果也会有较大的误差。

因此，一般规定：不得在之间。

1.2 测边交会测得未知点P至两已知点A、B的水平距离为S1和S2。

根据已知点坐标可反算出AB间的坐标方位角和边长S0。

在△ABP中用余弦定理可求得。

当两个已知点A、B之间互不通视的情况下，由未知点P测得与两已知点A、B之间的夹角θ和PA的距离S1，从而便可求出为未知点P的坐标。

距离后方交会计算casio4800程序

06-9-14 发布人：yangsixo 作者：yangsixo 本文已被浏览1033 次查看作者其它文章<<返回距离后方交会带高程的4850程序JLHFJHK：A“X1”：B“Y1”：H“H1”：C“X2”：D“Y2”：Lb1 1：｛MNO｝：M“JL1”：N“JL2”： O“GC”： J=0：I=0：PoL（C－A，D－B）：J＜0＝＞J=J＋360△Q=cos-1（（MM＋II－NN）÷（2×M×I））：E=J＋Q：X=A＋M×cosE：“X=”：X ▲Y=B＋M×sinE：“Y=”：Y ▲Z=H＋K－O：“Z=”：Z ▲ goto 1K为菱镜常数；X1、Y1、为仪器前方左边坐标；X2、Y2、为仪器前方右边坐标；H1为为仪器前方左边点高程；JL1为仪器前方左边距离；JL2为仪器前方右边距离；GC为仪器前方左边实测高差；X、Y计算坐标；Z计算高程。

例：左边X1=5.659 Y1=12375.053，右边X2=-6.102 Y2=12374.736，左边距离JL1=42.635右边距离JL2=42.407.得X=0.0892 Y=12332.7834.距离后方交会计算(CASIO fx-4800P计算器)程序距离后方交会计算(CASIO fx-4800P计算器)程序一、程序功能本程序适用于在一个未知点上设测站，观测两到个已知点的距离后，解算该未知坐标。

本程序也可以在CASIO fx-4500P计算器及 CASIO fx-4850P计算器上运行。

注意：这种观测两到个已知点的距离后解算该未知坐标的方法，缺少多余观测值，也就缺少检核条件。

二、源程序Lbl 1：{ABCDEFQ}：A"XA"：B” YA”：C"XB"：D"YB"：E"D1"：F"D2"：Q：J=0：G=Pol(C-A, D-B) ：H=J+QCos-1((GG+EE-FF)÷2÷G÷E)：X"XP"=A+ECosH◢Y"YP"=B+ESinH◢Goto 1←┘注：CASIO fx-4850改如下Lbl 1：{ABCDEFQ}：A"XA"：B” YA”：C"XB"：D"YB"：E"D1"：F"D2"：Q：J=0：G=Pol(C-A, D-B) ：H=J+QCos-1((GG+EE-FF)÷2÷G÷E)："XP"：X=A+ECosH◢"YP"：Y=B+ESinH◢Goto 1←┘三、使用说明1、规定(1) 未知点为P点，已知点分别为A点、B点；(2) P点至A点的距离为D1，P点至B点的距离为D2；(3) 当A、B、P三点逆时针排列时，Q=-1；当A、B、P三点顺时针排列时，Q=1。

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

#include<stdio.h>#include<stdlib.h>#include<malloc.h>#include<math.h>#include<conio.h>#define N 4#define T 1.41421void turn(double *A,double A2[],int m,int n) //计算矩阵的转置{ int i,j;for(i=0;i<m;i++)for(j=0;j<n;j++)A2[j*m+i]=A[i*n+j];}void mulAB(double *A,double *B,double *C,int am,int an,int bm,int bn) //计算两矩阵相乘{int i,j,l,u;if(an!=bm){printf("error!cannot do the multiplication.\n");return;}for(i=0;i<am;i++)for(j=0;j<bn;j++){u=i*bn+j;C[u]=0.0;for(l=0;l<an;l++)C[u]+=A[i*an+l]*B[l*bn+j];}return;}double *inv(double *a,int n) //计算矩阵的逆，本程序的难点,采用高斯约旦全选主元法{int *is,*js,i,j,k,l,u,v;double d,p;is=(int*)malloc(n*sizeof(int));js=(int*)malloc(n*sizeof(int));for (k=0; k<=n-1; k++){d=0.0;for (i=k;i<n;i++)for (j=k;j<n;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("error not inv\n");return NULL;}if (is[k]!=k)for (j=0;j<n;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;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;j++)if (j!=k){u=k*n+j;a[u]=a[u]*a[l];}for (i=0;i<n;i++)if (i!=k)for (j=0;j<n;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;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;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 a;}main()//主函数，空间后方交会的计算{FILE *fp; //定义一个文件指针long m;int i,j=0,it;double G[1000];doublef,t,w,k,limit=0,S1=0.0,S2=0.0,S3=0.0,x[N]={0},y[N]={0},x0[N]={0},y0[N]={0},X[N]={0},Y[N] ={0},Z[N]={0},Xs0,Ys0,Zs0;doublea[3],b[3],c[3],A[2*N*6],AT[6*2*N],ATA[6*6],*ATA_=NULL,l[2*N],ATl[6],V[6]={0};double F[6],Qx[6][6],Mi[6][6];if((fp=fopen("e:\\shuju.txt","r"))==NULL) //使fp指向被打开的shuju.txt 文件{ //并判断文件是否打开成功printf("\nerror on open shuju.txt\n");getch();exit(1);}for(i=0;i<N;i++)fscanf(fp,"%lf%lf%lf%lf%lf",&x[i],&y[i],&X[i],&Y[i],&Z[i]); //将文件中的数据赋给主函数定义的变量printf("原始数据:\n");printf("x\t\ty\t\tX\t\tY\t\tZ\t\t\n");//输出文件中的原始数据for(i=0;i<N;i++)printf("%.3lf\t\t%.3lf\t\t%.3lf\t%.3lf\t%.3lf\n",x[i],y[i],X[i],Y[i],Z[i]);printf("\n请输入摄影机主距和摄影比例尺分母；f,m:");scanf("%lf%ld",&f,&m); //输入f,mf=f/1000.0;for(i=0;i<N;i++){ x[i]/=1000.0;y[i]/=1000.0;S1+=X[i];S2+=Y[i];S3+=Z[i];}Xs0=S1/N;Ys0=S2/N;Zs0=m*f+S3/N;t=0.0;w=0.0;k=0.0; //计算外方位元素的初始值while(1){printf("\n---------------------------------第%d次计算------------------------------\n",j+1);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); //计算旋转矩阵for(i=0;i<N;i++){x0[i]=-f*(a[0]*(X[i]-Xs0)+b[0]*(Y[i]-Ys0)+c[0]*(Z[i]-Zs0))/(a[2]*(X[i]-Xs0)+b[2]*(Y[i]-Ys0)+c [2]*(Z[i]-Zs0));y0[i]=-f*(a[1]*(X[i]-Xs0)+b[1]*(Y[i]-Ys0)+c[1]*(Z[i]-Zs0))/(a[2]*(X[i]-Xs0)+b[2]*(Y[i]-Ys0)+c [2]*(Z[i]-Zs0));//计算像点坐标近似值G[i]=a[2]*(X[i]-Xs0)+b[2]*(Y[i]-Ys0)+c[2]*(Z[i]-Zs0);}for(i=0;i<N;i++){A[i*12+0]=(a[0]*f+a[2]*x[i])/G[i];A[i*12+1]=(b[0]*f+b[2]*x[i])/G[i];A[i*12+2]=(c[0]*f+c[2]*x[i])/G[i];A[i*12+3]=y[i]*sin(w)-(x[i]*(x[i]*cos(k)-y[i]*sin(k))/f+f*cos(k))*cos(w);A[i*12+4]=-f*sin(k)-x[i]*(x[i]*sin(k)+y[i]*cos(k))/f;A[i*12+5]=y[i];A[i*12+6]=(a[1]*f+a[2]*y[i])/G[i];A[i*12+7]=(b[1]*f+b[2]*y[i])/G[i];A[i*12+8]=(c[1]*f+c[2]*y[i])/G[i];A[i*12+9]=-x[i]*sin(w)-(y[i]*(x[i]*cos(k)-y[i]*sin(k))/f-f*sin(k))*cos(w);A[i*12+10]=-f*cos(k)-y[i]*(x[i]*sin(k)+y[i]*cos(k))/f;A[i*12+11]=-x[i];l[i*2+0]=x0[i]-x[i];l[i*2+1]=y0[i]-y[i]; //计算误差方程的系数阵以及lx,ly}// printf("output matrix: A\n");// printmatrix(A,2*N,6);// printf("output matrix: l\n");// printmatrix(l,2*N,1);turn(A,A T,2*N,6); //计算AT// printf("output matrix: AT\n");// printmatrix(AT,6,2*N);mulAB(AT,A,ATA,6,2*N,2*N,6); //计算ATA，组法方程ATA_=inv(ATA,6); //计算ATA的逆，中间量int p;int cnt=-1;for(it=0;it<36;it++){p=it%6;if(it%6==0){cnt++;}Qx[cnt][p++]=A TA_[it];}for(int it=0;it<6;it++){for(int jt=0;jt<6;jt++){if(it!=jt){Qx[it][jt]=0;//提取Qx的主对角线元素=Qii}// printf("%-10.3lf ",Qx[it][jt]);}// printf("\n");}mulAB(A T,l,ATl,6,2*N,2*N,1); //计算常数项ATL// printf("outpit matrinx: ATl\n");// printmatrix(ATl,6,1);mulAB(A TA_,ATl,V,6,6,6,1); //解法方程，求改正数，// printf("output matrix: V\n");// printmatrix(V,6,1);Xs0+=V[0];Ys0+=V[1];Zs0+=V[2];t+=V[3];w+=V[4];k+=V[5];for(i=0;i<6;i++){F[i]=V[i]/T;//m0}printf("第%d次计算的外方位元素为：\n",++j);printf("Xs=%.5lf\tYs=%.5lf\tZs=%.5lf\nt=%.5lf\tw=%.5lf\tk=%.5lf\n",Xs0,Ys0,Zs0,t,w,k);if(Xs0-limit<=0.0001&&Xs0-limit>=-0.0001){ //控制迭代次数break;}limit=Xs0;}printf("\n外方位元素为：\n");printf("Xs=%.5lf\tYs=%.5lf\tZs=%.5lf\nt=%.5lf\tw=%.5lf\tk=%.5lf\n",Xs0,Ys0,Zs0,t,w,k);printf("\n中误差为:\n");for(i=0;i<6;i++){for(j=0;j<6;j++){Mi[i][j]=F[i]*(sqrt(Qx[i][j]));//mi=m0*Qii开根号printf("%-13.10lf",Mi[i][j]);}printf("\n");}fclose(fp);return 0;}。