2014全国大学生数学建模a题

附录

求经纬度，matlab代码： function E=diedai(E0) R=1737 u=398610 a=R+50+7.5 c=a-R-15

b=sqrt(a^2-c^2) e=sqrt(1-(b/a)^2) M0=3;E0=1

for k=1:200000

M1=E0-e*sin(E0)

E0=E0+(pi-M1)./(1-e*sin(E0)) M0=M1

if abs(M1-M0)<=0.0000001 break; end end E=E0

T=2*pi*sqrt((a^3)/u)

d=2*atan(sqrt((1+e)/(1-e))*tan(E/2)) w=((1/27.32166)*360)/(3600*24)

s=19.51+180*atan(cos(44.12/pi)*tan(d))-w*T/2 z=asin(sin(44.12/pi)*sin(d))

%以预着陆点（19.51，44.12）为标准远月点的经纬度表示%

function E=diedai(E0) M0=0;E0=0.7 for k=1:20000

M1=E0-0.0237*sin(E0)

E0=E0+(3.14-M1)./(1-0.0237*sin(E0)) M0=M1

if abs(M1-M0)<=0.00001 break; end end E=E0 %迭代

螺旋式搜索算法：

X=imread('附件3 距2400m处的数字高程图.tif'); R=10; %扫描步长

k=tan(pi/6); %斜率控制，倾斜角度不能超过pi/6

15

[ Xe ] = roughAvoid( X, R, k ); %趋避障函数，可得矩阵Xe imwrite(Xe,'Xe.bmp','bmp'); %将Xe矩阵存储成Xe.bmp文件

function [ C ] = isFlat( Y, k )

% 判断Y矩阵是否平坦，如果平坦返回1，如果不平坦，返回0 %

[m,n]=size(Y); C=ones(size(Y)); % Yl=Y(:,1);

% Yr=flipud(Y(:,m)); for i=1:m

j=m-i+1;

kc=(Y(j,m)-Y(i,1))/sqrt((j-i)^2+m^2); if abs(kc)>k

C=zeros(size(Y)); return; end end

for i=1:n

j=n-i+1;

kc=(Y(n,j)-Y(1,i))/sqrt((j-i)^2+n^2); if abs(kc)>k

C=zeros(size(Y)); return; end end end

function [ Xe ] = roughAvoid( X, R, k )

% 粗避障,对矩阵X的元素按边长为2R+1的区域进行避障搜索

% 先在X中以某个点(m.,n)为中心取出(2R+1)*(2R+1)矩阵X(m-R:m+R,n-R:n+R) % 求斜率 % X矩阵

% R控制搜索区域边长 % k最大斜率控制

[mX, nX]=size(X); Xe=ones(mX, nX);

for m=R+1:R:mX-R %以步长为R进行纵向扫描

16

for n=R+1:R:nX-R %以步长为R进行横向扫描

[ Xe(m-R:m+R,n-R:n+R) ] = isFlat( X(m-R:m+R,n-R:n+R), k );

end end end

阀值分割程序：

I = imread('附件3 距2400m处的数字高程图.tif'); figure, imshow(I), title('原图'); C=histc(I,0:255); n=sum(C'); N=sum(n); t=n/N;

[p,threshold]=min(t(120:150)); threshold=threshold+120; tt=find(I>threshold); I(tt)=0;

tt=find(I>=threshold); I(tt)=255; figure;

imshow(I);

title('灰度阀值分割图');

仿真图代码： n=0;s=0;w=0;e=0; for i=1:1149

if Xe(1150-i,1150)==1 w=w+1; else

continue; end end

for i=1:1150

if Xe(1150+i,1150)==1 e=e+1; else

continue; end end

for i=1:1149

if Xe(1150,1150-i)==1

17

s=s+1; else

continue; end end

for i=1:1149

if Xe(1150,1150+i)==1 n=n+1; else

continue; end end

for i=1:1149

if Xe(1150,1150+i)==1 n=n+1; else

continue; end end

18