基于matlab的坐标转换 下载本文

内容发布更新时间 : 2024/12/23 4:25:53星期一 下面是文章的全部内容请认真阅读。

河南理工大学本科毕业论文

参考文献

[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