作业报告
空间方交会
专 业： 测绘工程
班 级： 2008级X班
姓 名：
指导教师：
2010 年X月X日
1 作业务 3
2 作业思想 3
3 作业条件数 3
4 作业程 3
5 源程序 4
6 计算结果 17
7心体会建议 17
1 作业务
计算似垂直摄影情况方交会解利摄影测量空间方交会方法获取相片6外方位元素限差01
2作业思想
利摄影测量空间方交会方法求解该方法基思想利少三直面控制点坐标A（XAYAZA）B（XBYBZB）C（XCYCZC）影应三点影坐标a（xaya）b（xbyb）c（xcyc）根线方程反求该相片外方位元素XSYSZSφωκ
3作业条件数
已知摄影机距f15324mm四点点坐标相应面坐标列入表：
表1
点号
点坐标
面坐标
x（mm）
y(mm)
X(m)
Y(m)
Z(m)
1
8615
6899
3658941
2527332
219517
2
5340
8221
3763108
3132451
72869
3
1478
7663
3910097
2493498
238650
4
1046
6443
4042654
3031981
75731
4作业程
4．1 获取已知数
相片例尺1m1：10000方位元素f15324mmx0y0获取控制点面测量坐标XtYtZt
4．2 量测控制点点坐标：
次作业中已知见表1
4．3 确定未知数初始值：
似垂直摄影情况胶原素初始值0φ0 ω0 κ00线元素中ZS0Hmf15324mXS0YS0取值四控制点坐标均值：
XS0 3843700
YS0 8910662
4．4 计算旋转矩阵R：
利胶原素似值计算方余弦值组成R阵
4．5 逐点计算点坐标似值：
利未知数似值线方程式计算控制点点坐标似值（x）（y）
4．6 组成误差方程：
逐点计算误差方程式系数常数项
4．7 组成法方程式：
计算法方程系数矩阵ATA常数项ATL
4．8 求解外方位元素：
根法方程式X（AtA）1 ATL解求外方位元素改正数相应似值求外方位元素新似值
4．9 求解外方位元素：
求外方位元素改正数规定限差（01）较限差计算终止否新似值重复第4448步骤计算知道满足求止
5 源程序
#include
#include
#include
const double PRECISION1e5
typedef double DOUBLE[5]
int InputData(int &Num DOUBLE *&Datadouble &mdouble &f)
int Resection(const int &Numconst DOUBLE *&Dataconst double &mconst double &f)
int InverseMatrix(double *matrixconst int &row)
int main(int argc char* argv[])
{
DOUBLE *DataNULL
int Num
double f(0)m(0)
if(InputData(NumDatamf))
{
if (DataNULL)
{
delete []Data
}
return 1
}
if(Resection(NumDatamf))
{
if (DataNULL)
{
delete []Data
}
return 1
}
if (DataNULL)
{
delete []Data
}
printf(解算完毕＼n)
do{
printf(计算结果保存＼结果txt＼文件中＼n
请选择操作（输入P开结果数R开原始数退出程序）：)
fflush(stdin) 刷新输入流
char ordergetchar()
if (＇P＇order || ＇p＇order)
{
system(结果txt)
}
else if (＇R＇order || ＇r＇order)
{
system(datatxt)
}
else
break
system(cls)
}while(1)
system(PAUSE)
return 0
}
**********************************************
*函数名：InputData
*函数介绍：文件（datatxt）中读取数
*文件格式：
*点数 m(未知写作0)
* 方位元素(f x0 y0)
*编号 x y X Y Z
*实例：
4 0
15324 0 0
1 8615 6899 3658941 2527332 219517
2 5340 8221 3763108 3132451 72869
3 1478 7663 3910097 2493498 238650
4 1046 6443 4042654 3031981 75731
*参数：(inout)Num(点数)
*(inout)Data(存放数)mfx0y0
*返回值：int 0成功1文件开失败2控制点
*数足3文件格式错误
**********************************************
int InputData(int &Num DOUBLE *&Datadouble &mdouble &f)
{
double x0y0
FILE *fp_input
if ((fp_inputfopen(datatxtr)))
{
return 1
}
fscanf(fp_inputdlf&Num&m)
if (Num<4)
{
return 2
}
fscanf(fp_inputlflflf&f&x0&y0)
f1000
if (m<0 || f<0)
{
return 3
}
Datanew DOUBLE[Num]
double *temp new double[Num1]
double scale0
int i
for (i0i{
读取数忽略编号
if(fscanf(fp_input*dlflflflflf
&Data[i][0]&Data[i][1]&Data[i][2]
&Data[i][3]&Data[i][4])5)
{
return 3
}
单位换算成m
Data[i][0]10000
Data[i][1]10000
}
果m未知算值
if (0m)
{
for (i0i {
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]20
}
mscale(Num1)
}
fclose(fp_input)
delete []temp
return 0
}
**********************************************
*函数名：MatrixMul
*函数介绍：求两矩阵积
*参数：Jz1(第矩阵)row(第矩阵行数)
*Jz2(第二矩阵)row(第二矩阵列数)com(第
*矩阵列数)(out)JgJz(存放结果矩阵)
*返回值：void
**********************************************
void MatrixMul(double *Jz1const int &rowdouble *Jz2
const int &lineconst int &comdouble *JgJz)
{
for (int i0i{
for (int j0j {
double temp0
for (int k0k {
temp+*(Jz1+i*com+k)*(*(Jz2+k*line+j))
}
*(JgJz+i*line+j)temp
}
}
}
**********************************************
*函数名：OutPut
*函数介绍：结果txt文件输出数
*参数：Q协数阵m精度m0单位权中误差6外
*方位元素旋转矩阵
*返回值：int0成功1失败
**********************************************
int OutPut(const double *&Qconst double *&mconst double &m0
const double &Xsconst double &Ysconst double &Zs
const double &Phiconst double &Omega
const double &Kappaconst double *R)
{
FILE *fp_out
if ((fp_outfopen(结果txtw)))
{
return 1
}
FILE *fp_input
if ((fp_inputfopen(datatxtr)))
{
return 1
}
fprintf(fp_out**************************************
**************************************
**************************************
*********************************＼n)
fprintf(fp_out＼n空间方交会程序(C＼＼C++)＼n测绘班＼n
学号：20080729＼n姓名：陈闻亚＼n＼n)
fprintf(fp_out**************************************
**************************************
**************************************
*********************************＼n)
fprintf(fp_out已知数：＼n＼n已知点数：)
int numdouble tempxy
fscanf(fp_inputdlf&num&temp)
fprintf(fp_outd＼nnum)
fprintf(fp_out摄影例尺（0表示值位置）：)
fprintf(fp_out100lf＼ntemp)
fprintf(fp_out方位元素(f x0 y0)：)
fscanf(fp_inputlflflf&temp&x&y)
fprintf(fp_out10lf＼t10lf＼t10lf＼ntempxy)
for (int i0i{
double temp[5]
fscanf(fp_input*dlflflflflf
&temp[0]&temp[1]&temp[2]&temp[3]&temp[4])
fprintf(fp_out3d＼t10lf＼t10lf＼t10lf＼t10lf＼t10lf＼n
i+1temp[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＼tXs10lf＼nXs)
fprintf(fp_out＼tYs10lf＼nYs)
fprintf(fp_out＼tZs10lf＼nZs)
fprintf(fp_out＼tPhi10lf＼nPhi)
fprintf(fp_out＼tOmega10lf＼nOmega)
fprintf(fp_out＼tKappa10lf＼n＼nKappa)
fprintf(fp_out旋转矩阵：＼n)
for (i0i<3i++)
{
fprintf(fp_out＼t)
for (int j0j<3j++)
{
fprintf(fp_out10lf＼t*(R+i*3+j))
}
fprintf(fp_out＼n)
}
fprintf(fp_out＼n单位权中误差：10lf＼n＼nm0)
fprintf(fp_out协数阵：＼n)
for (i0i<6i++)
{
fprintf(fp_out＼t)
for (int j0j<6j++)
{
fprintf(fp_out20lf＼t*(Q+i*6+j))
}
fprintf(fp_out＼n)
}
fprintf(fp_out＼n外方位元素精度：)
for (i0i<6i++)
{
fprintf(fp_out10lf＼tm[i])
}
fprintf(fp_out＼n)
fprintf(fp_out**************************************
**************************************
**************************************
*********************************＼n)
fclose(fp_out)
return 0
}
**********************************************
*函数名：Resection
*函数介绍：计算
*参数：Num(点数)Data(数)mf(焦距)x0y0
*返回值：int0成功失败
**********************************************
int Resection(const int &Numconst DOUBLE *&Dataconst double &m
const double &f)
{
double Xs0Ys0Zs0
int ij
设置初始值
for (i0i{
Xs+Data[i][2]
Ys+Data[i][3]
}
XsNum
YsNum
Zsm*f
double Phi(0)Omega(0)Kappa(0)
double R[3][3]{00}
double *Lnew double[2*Num]
typedef double Double6[6]
Double6 *Anew Double6[2*Num]
double *ATnew double[2*Num*6]
double *ATAnew double[6*6]
double *ATLnew double[6]
double *Xgnew 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 (i0i {
double XR[0][0]*(Data[i][2]Xs)+R[1][0]*(Data[i][3]Ys)+
R[2][0]*(Data[i][4]Zs)
double YR[0][1]*(Data[i][2]Xs)+R[1][1]*(Data[i][3]Ys)+
R[2][1]*(Data[i][4]Zs)
double ZR[0][2]*(Data[i][2]Xs)+R[1][2]*(Data[i][3]Ys)+
R[2][2]*(Data[i][4]Zs)
double xxxyyy
xxxf*XZ
yyyf*YZ
常数项
L[2*i]Data[i][0](f*XZ)
L[2*i+1]Data[i][1](f*YZ)
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转置矩阵AT
for (i0i<2*Numi++)
{
for (j0j<6j++)
{
*(AT+j*2*Num+i)A[i][j]
}
}
求ATA
MatrixMul(AT6&A[0][0]62*NumATA)
if(InverseMatrix(ATA6))
return 1
MatrixMul(AT6L12*NumATL)
MatrixMul(ATA6ATL16Xg)
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 *QATA
double *Vnew double[2*Num]
MatrixMul(&A[0][0]2*NumXg16V)
double VTV0
for(i0i<2*Numi++)
{
V[i]L[i]
VTV+V[i]*V[i]
}
double m0sqrt(VTV(2*Num6))
double *mmnew double[6]
for (i0i<6i++)
{
mm[i]sqrt(*(Q+i*6+i))*m0
}
OutPut(Qmmm0XsYsZsPhiOmegaKappa&R[0][0])
delete []L
delete []A
delete []AT
delete []ATA
delete []ATL
delete []Xg
delete []mm
delete []V
return 0
}
void swap(double &adouble &b)
{
double tempa
ab
btemp
}
**********************************************
*函数名：InverseMatrix
*函数介绍：求矩阵逆（高斯约法）
*输入参数：(inout)matrix（矩阵首址）
*（in）row（矩阵阶数）
*输出参数：matrix（原矩阵逆矩阵）
*返回值：int 0成功1失败
*调函数：swap(double&double&)
**********************************************
int InverseMatrix(double *matrixconst int &row)
{
double *mnew double[row*row]
double *ptemp*ptm
int ij
ptempmatrix
for (i0i{
for (j0j {
*pt*ptemp
ptemp++
pt++
}
}
int k
int *isnew int[row]*jsnew int[row]
for (k0k{
double max0
全选元
寻找元素
for (iki {
for (jkj {
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 (i0i {
swap(*(m+k*row+i)*(m+is[k]*row+i))
}
}
列交换
if (js[k]k)
{
for (i0i {
swap(*(m+i*row+k)*(m+i*row+js[k]))
}
}
*(m+k*row+k)1(*(m+k*row+k))
for (j0j {
if (jk)
{
*(m+k*row+j)**((m+k*row+k))
}
}
for (i0i {
if (ik)
{
for (j0j {
if(jk)
{
*(m+i*row+j)*(m+i*row+k)**(m+k*row+j)
}
}
}
}
for (i0i {
if(ik)
{
*(m+i*row+k)*(*(m+k*row+k))
}
}
}
int r
恢复行列
for (rrow1r>0r)
{
if (js[r]r)
{
for (j0j {
swap(*(m+r*row+j)*(m+js[r]*row+j))
}
}
if (is[r]r)
{
for (i0i {
swap(*(m+i*row+r)*(m+i*row+is[r]))
}
}
}
ptempmatrix
ptm
for (i0i{
for (j0j {
*ptemp*pt
ptemp++
pt++
}
}
delete []is
delete []js
delete []m
return 0
}
6 计算结果
外方位元素：
Xs39795452297
Ys27476462210
Zs7572685927
φ 0003987
ω 0002114
κ 0067578
旋转矩阵：
0997709 0067534 0003987
0067526 0997715 0002114
0004121 0001840 0999990
单位权中误差： 0000007
7 心体会建议
挑战会激情成功便更满足感次作业便次挑战总体说空间方交会知识学时点懵书加深理解弄懂编程方面面问题较学期类似作业时没攻破编程题便Excel解决事没继续研究实遗憾直遇次作业样问题悔初没继续学解决实总结起较难点：编程实现矩阵运算尤矩阵求逆搜集资料时网搜索矩阵求逆代码实现找合适考虑运程序问题断调试断改进程难度实现文导入数计算文输出结果次作业然较波折收获解决算拖延久问题然方效仿弄懂试着步步理出说次提高
西南交通学作业
文档香网(httpswwwxiangdangnet)户传
《香当网》用户分享的内容，不代表《香当网》观点或立场，请自行判断内容的真实性和可靠性！
该内容是文档的文本内容，更好的格式请下载文档
作业报告
空间方交会
专 业： 测绘工程
班 级： 2008级X班
姓 名：
指导教师：
2010 年X月X日
1 作业务 3
2 作业思想 3
3 作业条件数 3
4 作业程 3
5 源程序 4
6 计算结果 17
7心体会建议 17
1 作业务
计算似垂直摄影情况方交会解利摄影测量空间方交会方法获取相片6外方位元素限差01
2作业思想
利摄影测量空间方交会方法求解该方法基思想利少三直面控制点坐标A（XAYAZA）B（XBYBZB）C（XCYCZC）影应三点影坐标a（xaya）b（xbyb）c（xcyc）根线方程反求该相片外方位元素XSYSZSφωκ
3作业条件数
已知摄影机距f15324mm四点点坐标相应面坐标列入表：
表1
点号
点坐标
面坐标
x（mm）
y(mm)
X(m)
Y(m)
Z(m)
1
8615
6899
3658941
2527332
219517
2
5340
8221
3763108
3132451
72869
3
1478
7663
3910097
2493498
238650
4
1046
6443
4042654
3031981
75731
4作业程
4．1 获取已知数
相片例尺1m1：10000方位元素f15324mmx0y0获取控制点面测量坐标XtYtZt
4．2 量测控制点点坐标：
次作业中已知见表1
4．3 确定未知数初始值：
似垂直摄影情况胶原素初始值0φ0 ω0 κ00线元素中ZS0Hmf15324mXS0YS0取值四控制点坐标均值：
XS0 3843700
YS0 8910662
4．4 计算旋转矩阵R：
利胶原素似值计算方余弦值组成R阵
4．5 逐点计算点坐标似值：
利未知数似值线方程式计算控制点点坐标似值（x）（y）
4．6 组成误差方程：
逐点计算误差方程式系数常数项
4．7 组成法方程式：
计算法方程系数矩阵ATA常数项ATL
4．8 求解外方位元素：
根法方程式X（AtA）1 ATL解求外方位元素改正数相应似值求外方位元素新似值
4．9 求解外方位元素：
求外方位元素改正数规定限差（01）较限差计算终止否新似值重复第4448步骤计算知道满足求止
5 源程序
#include
#include
#include
const double PRECISION1e5
typedef double DOUBLE[5]
int InputData(int &Num DOUBLE *&Datadouble &mdouble &f)
int Resection(const int &Numconst DOUBLE *&Dataconst double &mconst double &f)
int InverseMatrix(double *matrixconst int &row)
int main(int argc char* argv[])
{
DOUBLE *DataNULL
int Num
double f(0)m(0)
if(InputData(NumDatamf))
{
if (DataNULL)
{
delete []Data
}
return 1
}
if(Resection(NumDatamf))
{
if (DataNULL)
{
delete []Data
}
return 1
}
if (DataNULL)
{
delete []Data
}
printf(解算完毕＼n)
do{
printf(计算结果保存＼结果txt＼文件中＼n
请选择操作（输入P开结果数R开原始数退出程序）：)
fflush(stdin) 刷新输入流
char ordergetchar()
if (＇P＇order || ＇p＇order)
{
system(结果txt)
}
else if (＇R＇order || ＇r＇order)
{
system(datatxt)
}
else
break
system(cls)
}while(1)
system(PAUSE)
return 0
}
**********************************************
*函数名：InputData
*函数介绍：文件（datatxt）中读取数
*文件格式：
*点数 m(未知写作0)
* 方位元素(f x0 y0)
*编号 x y X Y Z
*实例：
4 0
15324 0 0
1 8615 6899 3658941 2527332 219517
2 5340 8221 3763108 3132451 72869
3 1478 7663 3910097 2493498 238650
4 1046 6443 4042654 3031981 75731
*参数：(inout)Num(点数)
*(inout)Data(存放数)mfx0y0
*返回值：int 0成功1文件开失败2控制点
*数足3文件格式错误
**********************************************
int InputData(int &Num DOUBLE *&Datadouble &mdouble &f)
{
double x0y0
FILE *fp_input
if ((fp_inputfopen(datatxtr)))
{
return 1
}
fscanf(fp_inputdlf&Num&m)
if (Num<4)
{
return 2
}
fscanf(fp_inputlflflf&f&x0&y0)
f1000
if (m<0 || f<0)
{
return 3
}
Datanew DOUBLE[Num]
double *temp new double[Num1]
double scale0
int i
for (i0i
读取数忽略编号
if(fscanf(fp_input*dlflflflflf
&Data[i][0]&Data[i][1]&Data[i][2]
&Data[i][3]&Data[i][4])5)
{
return 3
}
单位换算成m
Data[i][0]10000
Data[i][1]10000
}
果m未知算值
if (0m)
{
for (i0i
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]20
}
mscale(Num1)
}
fclose(fp_input)
delete []temp
return 0
}
**********************************************
*函数名：MatrixMul
*函数介绍：求两矩阵积
*参数：Jz1(第矩阵)row(第矩阵行数)
*Jz2(第二矩阵)row(第二矩阵列数)com(第
*矩阵列数)(out)JgJz(存放结果矩阵)
*返回值：void
**********************************************
void MatrixMul(double *Jz1const int &rowdouble *Jz2
const int &lineconst int &comdouble *JgJz)
{
for (int i0i
for (int j0j
double temp0
for (int k0k
temp+*(Jz1+i*com+k)*(*(Jz2+k*line+j))
}
*(JgJz+i*line+j)temp
}
}
}
**********************************************
*函数名：OutPut
*函数介绍：结果txt文件输出数
*参数：Q协数阵m精度m0单位权中误差6外
*方位元素旋转矩阵
*返回值：int0成功1失败
**********************************************
int OutPut(const double *&Qconst double *&mconst double &m0
const double &Xsconst double &Ysconst double &Zs
const double &Phiconst double &Omega
const double &Kappaconst double *R)
{
FILE *fp_out
if ((fp_outfopen(结果txtw)))
{
return 1
}
FILE *fp_input
if ((fp_inputfopen(datatxtr)))
{
return 1
}
fprintf(fp_out**************************************
**************************************
**************************************
*********************************＼n)
fprintf(fp_out＼n空间方交会程序(C＼＼C++)＼n测绘班＼n
学号：20080729＼n姓名：陈闻亚＼n＼n)
fprintf(fp_out**************************************
**************************************
**************************************
*********************************＼n)
fprintf(fp_out已知数：＼n＼n已知点数：)
int numdouble tempxy
fscanf(fp_inputdlf&num&temp)
fprintf(fp_outd＼nnum)
fprintf(fp_out摄影例尺（0表示值位置）：)
fprintf(fp_out100lf＼ntemp)
fprintf(fp_out方位元素(f x0 y0)：)
fscanf(fp_inputlflflf&temp&x&y)
fprintf(fp_out10lf＼t10lf＼t10lf＼ntempxy)
for (int i0i
double temp[5]
fscanf(fp_input*dlflflflflf
&temp[0]&temp[1]&temp[2]&temp[3]&temp[4])
fprintf(fp_out3d＼t10lf＼t10lf＼t10lf＼t10lf＼t10lf＼n
i+1temp[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＼tXs10lf＼nXs)
fprintf(fp_out＼tYs10lf＼nYs)
fprintf(fp_out＼tZs10lf＼nZs)
fprintf(fp_out＼tPhi10lf＼nPhi)
fprintf(fp_out＼tOmega10lf＼nOmega)
fprintf(fp_out＼tKappa10lf＼n＼nKappa)
fprintf(fp_out旋转矩阵：＼n)
for (i0i<3i++)
{
fprintf(fp_out＼t)
for (int j0j<3j++)
{
fprintf(fp_out10lf＼t*(R+i*3+j))
}
fprintf(fp_out＼n)
}
fprintf(fp_out＼n单位权中误差：10lf＼n＼nm0)
fprintf(fp_out协数阵：＼n)
for (i0i<6i++)
{
fprintf(fp_out＼t)
for (int j0j<6j++)
{
fprintf(fp_out20lf＼t*(Q+i*6+j))
}
fprintf(fp_out＼n)
}
fprintf(fp_out＼n外方位元素精度：)
for (i0i<6i++)
{
fprintf(fp_out10lf＼tm[i])
}
fprintf(fp_out＼n)
fprintf(fp_out**************************************
**************************************
**************************************
*********************************＼n)
fclose(fp_out)
return 0
}
**********************************************
*函数名：Resection
*函数介绍：计算
*参数：Num(点数)Data(数)mf(焦距)x0y0
*返回值：int0成功失败
**********************************************
int Resection(const int &Numconst DOUBLE *&Dataconst double &m
const double &f)
{
double Xs0Ys0Zs0
int ij
设置初始值
for (i0i
Xs+Data[i][2]
Ys+Data[i][3]
}
XsNum
YsNum
Zsm*f
double Phi(0)Omega(0)Kappa(0)
double R[3][3]{00}
double *Lnew double[2*Num]
typedef double Double6[6]
Double6 *Anew Double6[2*Num]
double *ATnew double[2*Num*6]
double *ATAnew double[6*6]
double *ATLnew double[6]
double *Xgnew 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 (i0i
double XR[0][0]*(Data[i][2]Xs)+R[1][0]*(Data[i][3]Ys)+
R[2][0]*(Data[i][4]Zs)
double YR[0][1]*(Data[i][2]Xs)+R[1][1]*(Data[i][3]Ys)+
R[2][1]*(Data[i][4]Zs)
double ZR[0][2]*(Data[i][2]Xs)+R[1][2]*(Data[i][3]Ys)+
R[2][2]*(Data[i][4]Zs)
double xxxyyy
xxxf*XZ
yyyf*YZ
常数项
L[2*i]Data[i][0](f*XZ)
L[2*i+1]Data[i][1](f*YZ)
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转置矩阵AT
for (i0i<2*Numi++)
{
for (j0j<6j++)
{
*(AT+j*2*Num+i)A[i][j]
}
}
求ATA
MatrixMul(AT6&A[0][0]62*NumATA)
if(InverseMatrix(ATA6))
return 1
MatrixMul(AT6L12*NumATL)
MatrixMul(ATA6ATL16Xg)
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 *QATA
double *Vnew double[2*Num]
MatrixMul(&A[0][0]2*NumXg16V)
double VTV0
for(i0i<2*Numi++)
{
V[i]L[i]
VTV+V[i]*V[i]
}
double m0sqrt(VTV(2*Num6))
double *mmnew double[6]
for (i0i<6i++)
{
mm[i]sqrt(*(Q+i*6+i))*m0
}
OutPut(Qmmm0XsYsZsPhiOmegaKappa&R[0][0])
delete []L
delete []A
delete []AT
delete []ATA
delete []ATL
delete []Xg
delete []mm
delete []V
return 0
}
void swap(double &adouble &b)
{
double tempa
ab
btemp
}
**********************************************
*函数名：InverseMatrix
*函数介绍：求矩阵逆（高斯约法）
*输入参数：(inout)matrix（矩阵首址）
*（in）row（矩阵阶数）
*输出参数：matrix（原矩阵逆矩阵）
*返回值：int 0成功1失败
*调函数：swap(double&double&)
**********************************************
int InverseMatrix(double *matrixconst int &row)
{
double *mnew double[row*row]
double *ptemp*ptm
int ij
ptempmatrix
for (i0i
for (j0j
*pt*ptemp
ptemp++
pt++
}
}
int k
int *isnew int[row]*jsnew int[row]
for (k0k
double max0
全选元
寻找元素
for (iki
for (jkj
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 (i0i
swap(*(m+k*row+i)*(m+is[k]*row+i))
}
}
列交换
if (js[k]k)
{
for (i0i
swap(*(m+i*row+k)*(m+i*row+js[k]))
}
}
*(m+k*row+k)1(*(m+k*row+k))
for (j0j
if (jk)
{
*(m+k*row+j)**((m+k*row+k))
}
}
for (i0i
if (ik)
{
for (j0j
if(jk)
{
*(m+i*row+j)*(m+i*row+k)**(m+k*row+j)
}
}
}
}
for (i0i
if(ik)
{
*(m+i*row+k)*(*(m+k*row+k))
}
}
}
int r
恢复行列
for (rrow1r>0r)
{
if (js[r]r)
{
for (j0j
swap(*(m+r*row+j)*(m+js[r]*row+j))
}
}
if (is[r]r)
{
for (i0i
swap(*(m+i*row+r)*(m+i*row+is[r]))
}
}
}
ptempmatrix
ptm
for (i0i
for (j0j
*ptemp*pt
ptemp++
pt++
}
}
delete []is
delete []js
delete []m
return 0
}
6 计算结果
外方位元素：
Xs39795452297
Ys27476462210
Zs7572685927
φ 0003987
ω 0002114
κ 0067578
旋转矩阵：
0997709 0067534 0003987
0067526 0997715 0002114
0004121 0001840 0999990
单位权中误差： 0000007
7 心体会建议
挑战会激情成功便更满足感次作业便次挑战总体说空间方交会知识学时点懵书加深理解弄懂编程方面面问题较学期类似作业时没攻破编程题便Excel解决事没继续研究实遗憾直遇次作业样问题悔初没继续学解决实总结起较难点：编程实现矩阵运算尤矩阵求逆搜集资料时网搜索矩阵求逆代码实现找合适考虑运程序问题断调试断改进程难度实现文导入数计算文输出结果次作业然较波折收获解决算拖延久问题然方效仿弄懂试着步步理出说次提高
西南交通学作业
文档香网(httpswwwxiangdangnet)户传
《香当网》用户分享的内容，不代表《香当网》观点或立场，请自行判断内容的真实性和可靠性！
该内容是文档的文本内容，更好的格式请下载文档
