河南理工大学本科毕业论文
参考文献
[1] 刘亚静,毛善君,郭达志,等.基于VC++的坐标系统转换设计与实现[J].湖南科技大学学报（自然科学版），2006，9（3）：61-64. [2]宁津生.现代大地测量参考系统[J].测绘学报，2002，（5）：7-11. [3]王沫然.MATLAB与科学计算[M].北京:电子工业出版社,2004.2.
[4]姚连璧，周小平.基于MATLAB的控制网平差程序设计[M].上海:同济大学出版社,2006.4.
[5]边少锋，柴洪洲，金际航，等.大地坐标系与大地基准[M].北京:国防工业出版社,2005.5. [6]孔祥元，郭际明，刘宗泉，等.大地测量学基础[M].武汉:武汉大学出版社,2006.1. [7]李开友. 基于MATLAB的工程运算可视化系统的设计与实现:(硕士学位论文). 昆明理工大学,2006.
[8]董钧祥，杨德宏，等.测量坐标转换模型及其应用[J].昆明理工大学学报(理工版)，2006，4(3):1-4.
[9]刘大杰,施一民，等.全球定位系统（GPS）的原理与数据处理[M].上海:同济大学出版社,1996.
[10]孔祥元,梅是义，等.控制测量学(下)[M].武汉:武汉大学出版社,2002. [11]施一民.现代大地控制测量[M].上海:同济大学出版社,2003.6.
[12]陈清礼,胡家华.大地坐标转换程序[J].江汉石油学院学报，1998，3(1):45-49. [13]范新云.利用GPS测定地方坐标系转换的四参数法[J].海洋测绘，2005，7(4):35-37. [14]Bowring B. R Transformation from Spatial to Geographic Coordinates. Survey Review. No.181,1987.
[15]Paul,M. K. A Note on Computation of Geographic Coordinates. Bulletin Geodesique. No.108,1978.
[16]李岳.坐标系统的设计与实现：[D].武汉：中国地质大学，2010.
[17]徐生望，周建国，王超，等.坐标转换模型问题研究[J].西部探矿工，2009，（2）：162-165.
[18]黄海礼，官云兰，韩子太，等.基于MATLAB的坐标转换实践[J].测绘科学，2012，27（1）：20-22.
[19]郑航行.坐标转换研究：[D].郑州：解放军信息工程大学，2007.
[20]刘大杰，陶本藻，等.实用测量数据处理方法[M].北京：测绘出版社，2000. [21]姚吉利.三维坐标转换参数直接计算的严密公式[J].测绘通报，2006，（5）：7-10.
31
河南理工大学本科毕业论文
附录
1.大地坐标转空间直角的matlab计算代码 %大地坐标转换成空间直角坐标函数dd2kj.m
% -------------------------------------------------------------------- function XYZ=dd2kj(a,Bl,BB) format long g global XYZ;
B=BB(:,1);L=BB(:,2);H=BB(:,3); %将输入数据分配给变量 b=a*(1-Bl);
e=sqrt(a^2-b^2)/a; B=d2r(B); L=d2r(L);
N=a./sqrt(1-sin(B).^2.*e^2); X=(N+H).*cos(B).*cos(L); Y=(N+H).*cos(B).*sin(L); Z=(N.*(1-e^2)+H).*sin(B); XYZ=[X Y Z];
% ------------------------------------------------------------------------- % ------------------------------------------------------------------------
2.空间直角坐标转为大地坐标matlab计算代码 %空间直角坐标转换为大地坐标转换函数kj2dd.m
% -------------------------------------------------------------------- function BLH=kj2dd(a,Bl,XYZ) format long g global BLH BLH=[];
X=XYZ(1,1);Y=XYZ(1,2);Z=XYZ(1,3);%将输入数据分配给变量 b=a*(1-Bl);
e=sqrt(a^2-b^2)/a;
L=atan(Y/X);%求大地经度 if L<0
L=L+pi; end
B(1)=atan((Z/sqrt(X^2+Y^2))*(1+e^2)); for i=1:4%迭代求大地纬度
B(i+1)=atan((1/sqrt(X^2+Y^2))*(Z+(a*e^2*tan(B(i))/(sqrt(1+(1-e^2)*tan(B(i))))))); end
Bn=B(i+1);
H=sqrt(X^2+Y^2)*cos(Bn)+Z*sin(Bn)-a/sqrt(1-sin(Bn)^2*e^2)*(1-sin(Bn)^2*e^2);%求大地高
32
河南理工大学本科毕业论文
Bn=r2d(Bn); L=r2d(L);
BLH=[BLH;Bn L H]; end
% -------------------------------------------------------------------- % --------------------------------------------------------------------
3.高斯正反算matlab计算代码
正算：
%高斯正算[x,y]=BLtoxy(B,L,L0,ellipse) function [x,y]=BLtoxy(B,L,L0,ellipse) switch(ellipse)
case 0 %WGS-84椭球 a=6378137;
b=6356752.3142;
case 1 T椭球 a=6378245;
b=6356863.0187730473;
case 2 ?椭球 a=6378140;
b=6356755.2881575287; end
format long %定义数据类型为双精度型
l=(L-L0)*60*60/206265 ; %把角度转为弧度 B=pi*B/180;
%求取第一二偏心率以及卯酉圈、和各个参数的数值 e1=sqrt(a^2-b^2)/a; e2=sqrt(a^2-b^2)/b;
N=a/sqrt(1-e1^2*sin(B)^2); t=tan(B); n=e2*cos(B);
a0=1+3/4*e1^2+45/64*e1^4+350/512*e1^6+11025/16384*e1^8;
a2=-(3/4*e1^2+60/64*e1^4+60/64*e1^4+525/512*e1^6+17640/16384*e1^8)/2; a4=(15/64*e1^4+210/512*e1^6+8820/16384*e1^8)/4; a6=-(35/512*e1^6+2520/16384*e1^8)/6; a8=(315/16384*e1^8)/8;
X=a*(1-e1^2)*(a0*B+a2*sin(2*B)+a4*sin(4*B)+a6*sin(6*B)+a8*sin(8*B));
x=X+1/2*N*cos(B)^2*l^2+1/24*N*t*(5-t^2+9*n^2+4*n^4)*cos(B)^4*l^4+1/720*N*t*(61-58*t^2+t^4)*cos(B)^6*l^6
y=N*cos(B)*l+1/6*N*(1-t^2+n^2)*cos(B)^3*l^3+1/120*N*(5-18*t^2+t^4+14*n^2-58*n^2*t^2)*cos(B)^5*l^5; y=y+500000 反算：
%高斯反算[B,L]=xytoBL(x,y,L0,ellipse)
33
河南理工大学本科毕业论文
function [B,L]=xytoBL(x,y,L0,ellipse) format long
'ellipse值：选择WGS-84椭球输入0，选择54椭球输入1，选择80椭球输入2' switch(ellipse)
case 0 %WGS-84椭球 a=6378137;
b=6356752.3142;
case 1 T椭球 a=6378245;
b=6356863.0187730473;
case 2 ?椭球 a=6378140;
b=6356755.2881575287; end
e1=sqrt(a^2-b^2)/a; e2=sqrt(a^2-b^2)/b; y=y-500000;
A0=1+3/4*e1^2+45/64*e1^4+350/512*e1^6+11025/16384*e1^8+43659/65536*e1^10; B0=X/a/(1-e1^2)/A0;
k0=1/2*(3/4*e1^2+45/64*e1^4+350/512*e1^6+11025/16384*e1^8); k2=-1/3*(63/64*e1^4+1108/512*e1^6+58239/16384*e1^8); k4=1/3*(604/512*e1^6+68484/16384*e1^8); k6=-1/3*(26328/16384*e1^8); si=sin(B0)^2;
Bf=B0+sin(2*B0)*(k0+si*(k2+si*(k4+k6*si))); tf=tan(Bf);
ankef=e2*cos(Bf); c=a/sqrt(1+ankef^2); Nf=c/Vf;
Mf=a*(1-e1^2)/(1-e1^2*sin(Bf))^1.5; an=ankef^2;
B=Bf-1/2*Vf^2*tf*(y^2/Nf^2-1/12*(5+3*tf^2+an-9*an*tf^2)*y^4/Nf^4+1/360*(61+90*tf^2+45*tf^4)*y^6/Nf^6);
l=y/Nf/cos(Bf)-y^3*(1+2*tf^2+an)/6/Nf^3/cos(Bf)+y^5*(5+28*tf^2+24*tf^4+6*an+8*an*tf^2)/120/Nf^5/cos(Bf); B=B*180/pi; l=l*180/pi; L=L0+l;
4.空间直角与空间直角之间坐标转换matlab计算代码 %七参数计算函数burse.m
% -------------------------------------------------------------------- function y = burse(X) global dX
X1=X(:,1);Y1=X(:,2);Z1=X(:,3); X2=X(:,4);Y2=X(:,5);Z2=X(:,6); N=size(X1); n=N(1);
B=zeros(3*n,7);
34
河南理工大学本科毕业论文
for i=1:n
%V(3*i:3*(i+1),1)=[VX2(i) VY2(i) Vy3(i)];
B((3*i-2):3*i,1:7)=-[1 0 0 X1(i) 0 -Z1(i) Y1(i) 0 1 0 Y1(i) Z1(i) 0 -X1(i) 0 0 1 Z1(i) -Y1(i) X1(i) 0]; %dX=[dX0 dX0 dX0 a1 a2 a3 a4]';
L(3*i-2:3*i,1)=[X2(i) Y2(i) Z2(i)]'; end
P=eye(3*n);
dX=-inv(B'*P*B)*B'*P*L
% -------------------------------------------------------------------- % --------------------------------------------------------------------
35