[MCM] 2017研究生数学建模竞赛A题绘图

最新推荐文章于 2022-08-19 11:26:27 发布

clemente620

最新推荐文章于 2022-08-19 11:26:27 发布

阅读量1.1k

点赞数 2

本文链接：https://blog.csdn.net/weixin_38130913/article/details/104471379

版权

clc,clear all;
filename = 'a.csv'; 
%%数据预处理
z = csvread(filename); %读取附件1 区域高程数据
y0=1:1:2913; %按38.2米为1个单位进行简化绘图
x0=1:1:2775;
[xx0,yy0]=meshgrid(x0,y0);
xx=xx0';%依题意 行列对应,进行转置
yy=yy0';

index=find(z<3000); %高程大于3000的矩阵索引 这里可以把小于改成大于,分成两个图来分析
z(index)=nan; %高程小于3000赋值为缺值NaN

zmin = floor(min(z(:))); %高程最小值向下取整
zmax = ceil(max(z(:)));  %高程最大值向上取整
zinc = (zmax - zmin) / 3; % 分母代表颜色层次数
layers = zmin:zinc:zmax; %每个等高线的具体高程值

figure
contourf(xx,yy,z,layers); %等高线按照layers分层绘图
colorbar('Ticks',[0,3000,3623,4245,4800]); %具体数值标记 可以参考layers
%% 标记重点区域的中心坐标(坐标以38.2米为一个单位,简化复杂度) 
point =[793.194,2350.79;
    1727.75,2217.28;
    2575.92,2007.85;
    1929.32,1596.86;
    1515.71,1246.07;
    2272.25,575.916;
    2450.26,1277.49;];

[m,n]=size(point); % 返回m行数值
%axis equal
r=10000/38.2; % 需要巡查的重点区域半径
d=2*r; % rectangle的边长
str = {'A','B','C','D','E','F','G'}; %重点区域编号
hold on;
for i =1:m  %重点区域标记m(行)次循环
    plot(point(i,1),point(i,2),'r.','MarkerSize',20); % 中心点标记
    rectangle('Position',[point(i,1)-r,point(i,2)-r,d,d],'Curvature',[1,1],'edgecolor','r','LineWidth',2); %重点区域 巡查范围
    text(point(i,1)+30,point(i,2),str(i),'Color',[1,0,0],'FontSize',20,'FontWeight','bold'); %重点区域编号
end
hold off;

%存图用命令 saveas(gcf,['e:\','test','.png']);  更高清

2017年研究生数模竞赛A题 3000高程以上绘图

clc,clear all;
% 程序运行计时开始
t0 = clock;

filename = 'a.csv'; 
%% 数据预处理
z = csvread(filename); %读取附件1 区域高程数据
y0=1:1:2913; %按38.2米为1个单位进行简化绘图
x0=1:1:2775;
[xx0,yy0]=meshgrid(x0,y0);
xx=xx0';%依题意 行列对应,进行转置
yy=yy0';

index=find(z>=3000); %高程大于3000的矩阵索引 这里可以把小于改成大于,分成两个图来分析
z(index)=nan; %高程大于3000赋值为缺损

zmin = floor(min(z(:))); %高程最小值向下取整
zmax = ceil(max(z(:)));  %高程最大值向上取整
zinc = (zmax - zmin) / 3; % 分母代表颜色层次数
layers = zmin:zinc:zmax; %每个等高线的具体高程值

figure
contourf(xx,yy,z,layers); %等高线按照layers分层绘图
colorbar('Ticks',[1200,1745,2372,3000,3623,4245,4800]); %具体数值标记 可以参考layers
%% 标记重点区域的中心坐标(坐标以38.2米为一个单位,简化复杂度) 
point =[793.194,2350.79;
    1727.75,2217.28;
    2575.92,2007.85;
    1929.32,1596.86;
    1515.71,1246.07;
    2272.25,575.916;
    2450.26,1277.49;];

[m,n]=size(point); % 返回m行数值
%axis equal
r=10000/38.2; % 需要巡查的重点区域半径
d=2*r; % rectangle的边长
str = {'A','B','C','D','E','F','G'}; %重点区域编号
%% 绘图
hold on;
for i =1:m  %重点区域标记m(行)次循环
    plot(point(i,1),point(i,2),'r.','MarkerSize',20); % 中心点标记
    rectangle('Position',[point(i,1)-r,point(i,2)-r,d,d],'Curvature',[1,1],'edgecolor','r','LineWidth',2); %重点区域 巡查范围
    text(point(i,1)+30,point(i,2),str(i),'Color',[1,0,0],'FontSize',20,'FontWeight','bold'); %重点区域编号
end
hold off;
Time_Cost = etime(clock,t0); %程序执行时间
disp(['程序执行时间:' num2str(Time_Cost),'秒']);
%存图用命令 saveas(gcf,['e:\','test','.png']);  更高清

2017年研究生数模竞赛A题 3000高程以下绘图

因为颜色比较混乱可考虑黑白灰重新设定

colormap(flipud(gray)); %flipud是反转！

更改颜色体系

clc,clear all;
% 程序运行计时开始
t0 = clock;

filename = 'a.csv'; 
%% 数据预处理
z = csvread(filename); %读取附件1 区域高程数据
y0=1:1:2913; %按38.2米为1个单位进行简化绘图
x0=1:1:2775;
[xx0,yy0]=meshgrid(x0,y0);
xx=xx0';%依题意 行列对应,进行转置
yy=yy0';

index=find(z>=3000); %高程大于3000的矩阵索引 这里可以把小于改成大于,分成两个图来分析
z(index)=nan; %高程大于3000赋值为缺损

zmin = floor(min(z(:))); %高程最小值向下取整
zmax = ceil(max(z(:)));  %高程最大值向上取整
zinc = (zmax - zmin) / 3; % 分母代表颜色层次数
layers = zmin:zinc:zmax; %每个等高线的具体高程值

figure
contourf(xx,yy,z,layers); %等高线按照layers分层绘图
%colorbar('Ticks',[1200,1745,2372,3000,3623,4245,4800]); %具体数值标记 可以参考layers
colormap(gray);
colorbar;
%% 标记重点区域的中心坐标(坐标以38.2米为一个单位,简化复杂度) 
point =[793.194,2350.79;
    1727.75,2217.28;
    2575.92,2007.85;
    1929.32,1596.86;
    1515.71,1246.07;
    2272.25,575.916;
    2450.26,1277.49;];

[m,n]=size(point); % 返回m行数值
%axis equal
r=10000/38.2; % 需要巡查的重点区域半径
d=2*r; % rectangle的边长
str = {'A','B','C','D','E','F','G'}; %重点区域编号
%% 绘图
hold on;
for i =1:m  %重点区域标记m(行)次循环
    plot(point(i,1),point(i,2),'r.','MarkerSize',20); % 中心点标记
    rectangle('Position',[point(i,1)-r,point(i,2)-r,d,d],'Curvature',[1,1],'edgecolor','r','LineWidth',2); %重点区域 巡查范围
    text(point(i,1)+30,point(i,2),str(i),'Color',[1,0,0],'FontSize',20,'FontWeight','bold'); %重点区域编号
end
hold off;
Time_Cost = etime(clock,t0); %程序执行时间
disp(['程序执行时间:' num2str(Time_Cost),'秒']);
%存图用命令 saveas(gcf,['e:\','test','.png']);  更高清

小于3000 灰色体系

clc,clear all;
% 程序运行计时开始
t0 = clock;

filename = 'a.csv'; 
%% 数据预处理
z = csvread(filename); %读取附件1 区域高程数据
y0=1:1:2913; %按38.2米为1个单位进行简化绘图
x0=1:1:2775;
[xx0,yy0]=meshgrid(x0,y0);
xx=xx0';%依题意 行列对应,进行转置
yy=yy0';

index=find(z<3000); %高程大于3000的矩阵索引 这里可以把小于改成大于,分成两个图来分析
z(index)=nan; %高程大于3000赋值为缺损

zmin = floor(min(z(:))); %高程最小值向下取整
zmax = ceil(max(z(:)));  %高程最大值向上取整
zinc = (zmax - zmin) / 3; % 分母代表颜色层次数
layers = zmin:zinc:zmax; %每个等高线的具体高程值

figure
contourf(xx,yy,z,layers); %等高线按照layers分层绘图
%colorbar('Ticks',[1200,1745,2372,3000,3623,4245,4800]); %具体数值标记 可以参考layers
colormap(flipud(gray));
colorbar;
%% 标记重点区域的中心坐标(坐标以38.2米为一个单位,简化复杂度) 
point =[793.194,2350.79;
    1727.75,2217.28;
    2575.92,2007.85;
    1929.32,1596.86;
    1515.71,1246.07;
    2272.25,575.916;
    2450.26,1277.49;];

[m,n]=size(point); % 返回m行数值
%axis equal
r=10000/38.2; % 需要巡查的重点区域半径
d=2*r; % rectangle的边长
str = {'A','B','C','D','E','F','G'}; %重点区域编号
%% 绘图
hold on;
for i =1:m  %重点区域标记m(行)次循环
    plot(point(i,1),point(i,2),'r.','MarkerSize',20); % 中心点标记
    rectangle('Position',[point(i,1)-r,point(i,2)-r,d,d],'Curvature',[1,1],'edgecolor','r','LineWidth',2); %重点区域 巡查范围
    text(point(i,1)+30,point(i,2),str(i),'Color',[1,0,0],'FontSize',20,'FontWeight','bold'); %重点区域编号
end
hold off;
Time_Cost = etime(clock,t0); %程序执行时间
disp(['程序执行时间:' num2str(Time_Cost),'秒']);
%存图用命令 saveas(gcf,['e:\','test','.png']);  更高清

高于3000 灰色体系

clc,clear all;
% 程序运行计时开始
t0 = clock;

filename = 'a.csv'; 
%% 数据预处理
z = csvread(filename); %读取附件1 区域高程数据
y0=1:1:2913; %按38.2米为1个单位进行简化绘图
x0=1:1:2775;
[xx0,yy0]=meshgrid(x0,y0);
xx=xx0';%依题意 行列对应,进行转置
yy=yy0';
% 手工设定完成颜色的分界 0-1500白色 
% 1500-3000浅灰色 3000-4200深灰 4200以上黑色
layers = [0,1500,3000,4200,5000]; 

figure
contourf(xx,yy,z,layers); %等高线按照layers分层绘图
%colorbar('Ticks',[1200,1745,2372,3000,3623,4245,4800]); %具体数值标记 可以参考layers
colormap(flipud(gray)); %高程由低到高 呈现白到黑 flipud是反转色
colorbar;
%% 标记重点区域的中心坐标(坐标以38.2米为一个单位,简化复杂度) 
point =[793.194,2350.79; 
    1727.75,2217.28;
    2575.92,2007.85;
    1929.32,1596.86;
    1515.71,1246.07;
    2272.25,575.916;
    2450.26,1277.49;];

[m,n]=size(point); % 返回m行数值
%axis equal
r=10000/38.2; % 需要巡查的重点区域半径
d=2*r; % rectangle的边长
str = {'A','B','C','D','E','F','G'}; %重点区域编号
%% 绘图
hold on;
for i =1:m  %重点区域标记m(行)次循环
    plot(point(i,1),point(i,2),'r.','MarkerSize',20); % 中心点标记
    rectangle('Position',[point(i,1)-r,point(i,2)-r,d,d],'Curvature',[1,1],'edgecolor','r','LineWidth',2); %重点区域 巡查范围
    text(point(i,1)+30,point(i,2),str(i),'Color',[1,0,0],'FontSize',20,'FontWeight','bold'); %重点区域编号
end
hold off;
Time_Cost = etime(clock,t0); %程序执行时间
disp(['程序执行时间:' num2str(Time_Cost),'秒']);
%存图用命令 saveas(gcf,['e:\','test','.png']);  更高清

不作过滤全高程绘图

3000以下是优先巡查区域 ;5000是无人机最大飞行高度 4200是假设恒定的飞行高度(优劣待定,会影响无人机扫描半径) 整个区域最大高程4867,最低1118

layers = [0,1500,3000,4200,5000]; % 手工设定完成颜色的分界 0-1500白色 1500-3000浅灰色 3000-4200深灰 4200以上黑色

colormap(flipud(gray)); %高程由低到高呈现白到黑 flipud是反转色

选定圆内数据选择性绘图可算圆内平均高程

clear all

filename = 'C:\matlab_code\a.csv'; 
%% 数据预处理
z = csvread(filename); %读取附件1 区域高程数据
y0=1:1:2913; %按38.2米为1个单位进行简化绘图
x0=1:1:2775;
[xx0,yy0]=meshgrid(x0,y0);
xx=xx0';%依题意 行列对应,进行转置
yy=yy0';

%% 高程数据的行列下标
row=2775;
col=2913;

%圆半径
r=10000/38.2; 
%圆内元素值及索引放置矩阵
index=[];

ZZ=zeros(2775,2913)*nan;

in=1; %索引循环重置
center=[793.194,2350.79]; % A圆心位置行列
for i=1:row
     for j=1:col
          coordinate=[i j];
          if norm(center-coordinate)<=r %norm为范数函数，默认2-范数，用来求距离center点位小于 半径r的所有点i,j
             ZZ(i,j)=z(i,j); %符合条件的元素值,赋值到新的ZZ,ZZ其余位置元素为NaN
             index(in,:)=coordinate; %符合条件的元素索引
             in=in+1;
          end
     end
end
meanA=nanmean(ZZ(:)) %算A区域内平均高程 由于ZZ是累加的,单独计算每个区域内平均高程要单独进行 nanmean是不包括nan计算平均值

in=1; % 索引循环重置
center=[1727.75,2217.28];% B圆心位置行列
for i=1:row
     for j=1:col
          coordinate=[i j];
          if norm(center-coordinate)<=r %norm为范数函数，默认2-范数，用来求距离center点位小于 半径r的所有点i,j
             ZZ(i,j)=z(i,j); %符合条件的元素值,赋值到新的ZZ,ZZ其余位置元素为NaN
             index(in,:)=coordinate; %符合条件的元素索引
             in=in+1;
          end
     end
end


in=1; % 索引循环重置
center=[2575.92,2007.85];% C圆心位置行列
for i=1:row
     for j=1:col
          coordinate=[i j];
          if norm(center-coordinate)<=r %norm为范数函数，默认2-范数，用来求距离center点位小于 半径r的所有点i,j
             ZZ(i,j)=z(i,j); %符合条件的元素值,赋值到新的ZZ,ZZ其余位置元素为NaN
             index(in,:)=coordinate; %符合条件的元素索引
             in=in+1;
          end
     end
end

in=1; % 索引循环重置
center=[1929.32,1596.86]; % D圆心位置行列
for i=1:row
     for j=1:col
          coordinate=[i j];
          if norm(center-coordinate)<=r %norm为范数函数，默认2-范数，用来求距离center点位小于 半径r的所有点i,j
             ZZ(i,j)=z(i,j); %符合条件的元素值,赋值到新的ZZ,ZZ其余位置元素为NaN
             index(in,:)=coordinate; %符合条件的元素索引
             in=in+1;
          end
     end
end

in=1; % 索引循环重置
center=[1515.71,1246.07]; % E圆心位置行列
for i=1:row
     for j=1:col
          coordinate=[i j];
          if norm(center-coordinate)<=r %norm为范数函数，默认2-范数，用来求距离center点位小于 半径r的所有点i,j
             ZZ(i,j)=z(i,j); %符合条件的元素值,赋值到新的ZZ,ZZ其余位置元素为NaN
             index(in,:)=coordinate; %符合条件的元素索引
             in=in+1;
          end
     end
end

in=1; % 索引循环重置
center=[2272.25,575.916]; % F圆心位置行列
for i=1:row
     for j=1:col
          coordinate=[i j];
          if norm(center-coordinate)<=r %norm为范数函数，默认2-范数，用来求距离center点位小于 半径r的所有点i,j
             ZZ(i,j)=z(i,j); %符合条件的元素值,赋值到新的ZZ,ZZ其余位置元素为NaN
             index(in,:)=coordinate; %符合条件的元素索引
             in=in+1;
          end
     end
end

in=1; % 索引循环重置
center=[2450.26,1277.49]; % G圆心位置行列
for i=1:row
     for j=1:col
          coordinate=[i j];
          if norm(center-coordinate)<=r %norm为范数函数，默认2-范数，用来求距离center点位小于 半径r的所有点i,j
             ZZ(i,j)=z(i,j); %符合条件的元素值,赋值到新的ZZ,ZZ其余位置元素为NaN
             index(in,:)=coordinate; %符合条件的元素索引
             in=in+1;
          end
     end
end
meanAll=nanmean(ZZ(:))

%% 绘图
layers = [0,1500,3000,4200,5000]; 
figure
contourf(xx,yy,ZZ,layers); %等高线按照layers分层绘图
colormap(flipud(gray)); %高程由低到高 呈现白到黑 flipud是反转色

指定圆内数据进行绘图

第四问附件2 问题3,4中地面终端的地理位置

%% 第9章 蚁群算法及MATLAB实现——TSP问题
% 程序9-1
 
%% 数据准备
% 清空环境变量
clear all
clc
 
% 程序运行计时开始
t0 = clock;
 
% 导入数据
citys = xlsread('terminal.xlsx','B2:C73');
 
%% 计算城市间距离
n = size(citys,1); %城市数
D = zeros(n,n);
for i = 1:n
    for j = 1:n
        if i ~= j
            D(i,j) = sqrt(sum( ( citys(i,:) - citys(j,:) ).^2 ) ); %两点之间的距离
        else
            D(i,j) = 1e-4;          %为保证启发函数的分母不为0,设定的对角矩阵修正值为一个较小正值
        end
    end
end
 
%% 初始化参数
m = 72;                     % 蚂蚁数量最好接近城市数量的1.5倍
alpha = 1;                  % 信息素重要程度因子[1,4]最好
beta = 5;                   % 启发函数重要程度因子 5最好
vol = 0.2;                  % 信息素挥发(volatilization)因子 
Q = 10;                     % 常系数
Heu_F = 1./D;               % 启发函数(heuristic function)
Tau = ones(n,n);            % 信息素矩阵
Table = zeros(n,n);         % 路径记录表
iter = 1;                   % 迭代次数初值
iter_max = 300;             % 最大迭代次数 [100,500]
Route_best = zeros(iter_max,n);     % 各代最佳路径
Length_best = zeros(iter_max,1);    % 各代最佳路径的长度
Length_ave = zeros(iter_max,1);     % 各代路径的平均长度
Limit_iter = 0;                     % 程序收敛时迭代次数
 
%% 迭代寻找最佳路径
while iter <= iter_max
    % 随机产生各个蚂蚁的起点城市
    start = zeros(m,1);
    for i = 1:m
        temp = randperm(n); %randperm函数打乱顺序 1-n随机排序
        start = temp(1);
    end
    Table(:,1) = start; %路径记录表
    % 构建解空间
    citys_index = 1:n;
    % 逐个蚂蚁路径选择
    for i =1:m
        % 逐个城市路径选择
        for j = 2:n
            tabu = Table(i,1:(j - 1));                  % 已访问的城市集合（禁忌表）
            allow_index = ~ismember(citys_index,tabu);  % 1.ismember函数判断一个变量中的元素是否在另一个变量中出现,返回0-1矩阵；
            allow = citys_index(allow_index);           % 待访问的城市集合
            P = allow;
            % 计算城市间转移概率
            for k = 1:length(allow)
                P(k) = Tau(tabu(end),allow(k))^alpha *  Heu_F(tabu(end),allow(k))^beta;  %线路选择概率的分子
            end
            P = P / sum(P); %线路选择概率的分母
            % 轮盘赌法选择下一个访问城市
            Pc = cumsum(P);             %cumsum函数用于求变量中累加元素的和,如A=[1,2,3,4,5],那么cumsum(A)=[1,3,6,10,15]
            target_index = find(Pc >= rand);
            target = allow(target_index(1));
            Table(i,j) = target;
        end
    end
    % 计算各个蚂蚁的路径距离
    Length = zeros(m,1);
    for i = 1:m
        Route = Table(i,:);
        for j = 1:(n - 1)
            Length(i) = Length(i) + D(Route(j),Route(j + 1));
        end
        Length(i) = Length(i) + D(Route(n),Route(1)); %最后回到起点
    end
    % 计算最短路径距离及平均距离
    if iter == 1
        [min_Length, min_index] = min(Length);
        Length_best(iter) = min_Length;
        Length_ave(iter) = mean(Length);
        Route_best(iter,:) = Table(min_index,:);
        Limit_iter = 1;
    else
        [min_Length,min_index] = min(Length);
        Length_best(iter) = min(Length_best(iter - 1),min_Length);
        Length_ave(iter) = mean(Length);
        if Length_best(iter) == min_Length
            Route_best(iter,:) = Table(min_index,:);
            Limit_iter = iter;
        else
            Route_best(iter,:) = Route_best((iter - 1),:);
        end
    end
    % 更新信息素
    Delta_Tau = zeros(n,n); %信息素增量
    % 逐个蚂蚁计算
    for i = 1:m
        % 逐个城市计算
        for j = 1:(n - 1)
            Delta_Tau(Table(i,j),Table(i,j+1)) = Delta_Tau(Table(i,j),Table(i,j+1)) + Q/Length(i); %Q为常数
        end
        Delta_Tau(Table(i,n),Table(i,1)) = Delta_Tau(Table(i,n),Table(i,1)) + Q/Length(i); %最后回到第一个城市的最终值
    end
    Tau = (1 - vol) * Tau + Delta_Tau; %信息素挥发损失的剩下部分+新的增量
    % 迭代次数加1,清空路径记录表
    iter = iter + 1;
    Table = zeros(m,n);
end
 
%% 结果显示
[Shortest_Length,index] = min(Length_best);
Shortest_Route = Route_best(index,:);
Time_Cost = etime(clock,t0);
disp(['最短距离:' num2str(Shortest_Length)]);
disp(['最短路径:' num2str( [Shortest_Route Shortest_Route(1)] )]);
disp(['收敛迭代次数:' num2str(Limit_iter)]);
disp(['程序执行时间:' num2str(Time_Cost),'秒']);
 
%% 绘图
set(gca,'LineWidth',1.5); %边框加粗,美观
figure(1)
%  假设各个城市的X坐标为zuobiao_X，Y坐标为zuobiao_Y，zuobiao_X(i)表示第i个城市的横坐标，一共有n个城市，那么，采用以下循环语句进行画图：
%  for i=1:n-1
%  plot([zuobiao_X(i),zuobiao_X(i+1)],[zuobiao_Y(i),zuobiao_Y(i+1)],'-r');
%  hold on;
%  end
plot([ citys(Shortest_Route,1);citys(Shortest_Route(1),1) ], [ citys(Shortest_Route,2);citys(Shortest_Route(1),2) ], 'k.-','Markersize',20);
set(gca,'LineWidth',1.5); %边框加粗,美观
grid on;
for i = 1:size(citys,1)
    text(citys(i,1),citys(i,2),['  ' num2str(i)]);
end
text(citys(Shortest_Route(1),1),citys(Shortest_Route(1),2),'    起点');
text(citys(Shortest_Route(end),1),citys(Shortest_Route(end),2),'    终点')
xlabel('城市位置横坐标');
ylabel('城市位置纵坐标');
title(['ACA最优化路径(最短距离:' num2str(Shortest_Length) ')']);
figure(2);
plot(1:iter_max,Length_best,'b');
set(gca,'LineWidth',1.5); %边框加粗,美观
legend('最短距离');
xlabel('迭代次数');
ylabel('距离');
title('算法收敛轨迹');
set(gca,'LineWidth',1.5);  %边框加粗,美观
axis equal

蚁群算法整体路径规划

为了保持x轴从0开始命令行输入 axis equal

%%
clc;
clear all;
%区域标记 

zz=xlsread('区域高程数据.xlsx');
zz=zz';
x=(0:1:size(zz,2)-1)*38.2/1000; % 在 x 轴上取点
y=(0:1:size(zz,1)-1)*38.2/1000; % 在 y 轴上取点
layers = [0,1500,3000,4200,5000]; 
[xx,yy]=meshgrid(x, y);
[c,h]=contourf(xx, yy, zz,layers);
colormap(flipud(gray)); %高程由低到高 呈现白到黑 flipud是反转色

P=[30.3 89.8;
    66.0 84.7; 
    98.4 76.7; 
    73.7 61.0; 
    57.9 47.6; 
    86.8 22.0; 
    93.6 48.8];
%%
set(h,'ShowText','off','LevelList',[0,3000,4150])%设定等高线的值
hold on
plot(P(:,1),P(:,2),'s','MarkerEdgeColor','k','MarkerFaceColor','m','MarkerSize',10); %重点点位
plot(110,0,'s','MarkerEdgeColor','k','MarkerFaceColor','g','MarkerSize',10); %基地

for i=1:size(P,1)
    rectangle('Position',[P(i,1)-10, P(i,2)-10, 2*10, 2*10],'Curvature',[1,1],'EdgeColor','r','LineWidth',3);   %重点范围范围
    axis equal;
end


%% 找出每个 P 点周围的方格
% zav=sum(zz(1==(zz<=3000)))/sum(sum((zz<=3000))); 
zav=3000;
cover=2*(4150-zav)*tan(30/180*pi)/1000;   %无人机覆盖半径
W1=Wfind(P(1,1),P(1,2),10,cover);  % R=10 储存可扫描区域的点位  A区有sum(W1(:,3))个点位
W2=Wfind(P(2,1),P(2,2),10,cover); 
W3=Wfind(P(3,1),P(3,2),10,cover); 
W4=Wfind(P(4,1),P(4,2),10,cover); 
W5=Wfind(P(5,1),P(5,2),10,cover); 
W6=Wfind(P(6,1),P(6,2),10,cover); % 没有符合的点
W7=Wfind(P(7,1),P(7,2),10,cover); %有一个点3000以下

%% 判断方格是否符合条件
W1f=W1(find((W1(:,4)>0.5)==1),:);
W2f=W2(find((W2(:,4)>0.5)==1),:);
% W2f=W2f(((W2f(:,1)-110).^2+(W2f(:,2)-60).^2).^0.5<52,:); 圆圈约束条件
W3f=W3(find((W3(:,4)>0.5)==1),:);
W4f=W4(find((W4(:,4)>0.5)==1),:);
W5f=W5(find((W5(:,4)>0.7)==1),:);

Wxy=[W1f(:,[1,2])
    W2f(:,[1,2])
    W3f(:,[1,2])
    W4f(:,[1,2])
    W5f(:,[1,2])];  %筛选后的坐标点位
plot(Wxy(:,1),Wxy(:,2),'rs'); 
save('Wxy.mat','Wxy');
save('.\Wf.mat','W2f','W3f','W4f','W5f');


%计算覆盖率
sum(W1f(:,4))/sum(W1(:,3)) % A区
sum(W2f(:,4))/sum(W2(:,3)) % B区  分母是3000以下区域的面积 分子是有效扫描区域
sum(W3f(:,4))/sum(W3(:,3)) % C区
sum(W4f(:,4))/sum(W4(:,3)) % D区
sum(W5f(:,4))/sum(W5(:,3)) % E区
%%求总的覆盖率
sum_all=( sum(W1f(:,4))+sum(W2f(:,4))+sum(W3f(:,4))+sum(W4f(:,4))+sum(W5f(:,4))+0+0) /...
(sum(W1(:,3))+sum(W2(:,3))+sum(W3(:,3))+sum(W4(:,3))+sum(W5(:,3))) %分子是扫描有效点位 分母是符合3000米以下的点位数

main.m

function W= Wfind( x,y,R,cover )%W=[Wx,Wy,Wf],小网格的 x，y 坐标，R=10 半径，无人机覆盖半径

zz=xlsread('区域高程数据.xlsx');
zz=zz';
W=[];
n=ceil(R/cover); %ceil向上取整　n为重点区直径与无人机覆盖直径之比（大圈包含几个小圈）
for i=-n:1:n % 横坐标方向移动
    for j=-n:1:n %划分小网格的范围  纵坐标方向移动
            num=0;
            fnum=0;
            xn=x+i*cover; %小网格的坐标 
            yn=y+j*cover;
%小网格在地图中的可探测点 
zj=max(1,ceil((xn-cover/2)*1000/38.2)):min(size(zz,2)-1,floor((xn+cover/2)*1000/38.2));
        zi=ceil((yn-cover/2)*1000/38.2):floor((yn+cover/2)*1000/38.2); 
        for jj=zj
            for ii=zi 
                zx=38.2*(jj-1)/1000; 
                zy=38.2*(ii-1)/1000;
                 if zz(ii,jj)<=3000&&((zx-x)^2+(zy-y)^2)^0.5<=10 
                    num=num+1;

dn=[-(zz(ii,jj+1)-zz(ii,jj-1))/2/38.2,-(zz(ii+1,jj)-zz(ii-1,jj))/2/38.2,1];%法向量

                    dl=[xn-zx,yn-zy,(4150-zz(ii,jj))/1000]; 
                    if dn*dl'>=0&&atan(dl(3)/(dl(1)^2+dl(2)^2)^0.5)>=60/180*pi%判断不会被遮挡和在可视角范围内
                        fnum=fnum+1;
                    end
                 end
            end
        end
    if num>0

    W=[W; xn,yn,num/((zj(end)-zj(1)+1)*(zi(end)-zi(1)+1)),fnum/((zj(end)-zj(1)+1)*(zi(end)-zi(1)+1))]; %返回坐标和覆盖点数比例并累加
    end

end

end
end

Wfind.m

覆盖率

A 0.6937

B 0.8937

C 0.6726

D 0.7067

E 0.7561

总覆盖率 0.7706

次生灾害巡查分析

飞机对全局的 4000m 以下点进行巡查；

%%
clc,clear all;
%区域标记 clc ,clear all;
zz=xlsread('区域高程数据.xlsx');
zz=zz';
x=(0:1:size(zz,2)-1)*38.2/1000; % 在 x 轴上取点
y=(0:1:size(zz,1)-1)*38.2/1000; % 在 y 轴上取点
layers = [0,1500,3000,4200,5000]; 
[xx,yy]=meshgrid(x, y);
[c,h]=contourf(xx, yy, zz,layers);
colormap(flipud(gray)); %高程由低到高 呈现白到黑 flipud是反转色

set(h,'ShowText','off','LevelList',[0,4000,4150])%设定等高线的值 
hold on
%%
%找出每个 P 点周围的方格
% zav=sum(zz(1==(zz<=4000)))/sum(sum((zz<=4000))); 

zav=2500;
cover=2*(4150-zav)*tan(30/180*pi)/1000; 
R=max(x(end)/2,y(end)/2); 
W=Wfind(x(end)/2,y(end)/2,R,cover); 
save('W.mat','W');
%%

%判断方格是否符合条件 load W.mat

Wf=W(find((W(:,4)>0.3)==1),:);  %储存获取的方块坐标 
% Wf=Wf(~((Wf(:,1)<40 & Wf(:,2)<70)),:);

save('Wf.mat','Wf'); hold on
plot(Wf(:,1),Wf(:,2),'rs'); 
sum(Wf(:,3))/sum(W(find((W(:,3)>=1)==1),3))

main.m

function W= Wfind( x,y,R,cover )%W=[Wx,Wy,Wf],小网格的 x，y 坐标，重点区域半径, 扫描半径

zz=xlsread('区域高程数据.xlsx');
zz=zz';

W=[]; 
nx=ceil((size(zz,2)-1)*38.2/1000/cover/2); 
ny=ceil((size(zz,1)-1)*38.2/1000/cover/2); 
k=0;

for i=-nx:1:nx 
    k=k+1;
    for j=-ny:1:ny%划分小网格的范围 
        num=0;
        fnum=0;
        xn=x+i*cover; %小网格的坐标 
        yn=y+j*cover;
zj=max(2,ceil((xn-cover/2)*1000/38.2)):min(size(zz,2)-2,floor((xn+cover/2)* 1000/38.2));%小网格在地图中的可探测点
zi=max(2,ceil((yn-cover/2)*1000/38.2)):min(size(zz,1)-2,floor((yn+cover/2)* 1000/38.2));

        for jj=zj 
            for ii=zi
            zx=38.2*(jj-1)/1000; 
            zy=38.2*(ii-1)/1000; 
            if zz(ii,jj)<=4000
                num=num+1;
dn=[-(zz(ii,jj+1)-zz(ii,jj-1))/2/38.2,-(zz(ii+1,jj)-zz(ii-1,jj))/2/38.2,1];%法向量
            dl=[xn-zx,yn-zy,(4150-zz(ii,jj))/1000]; 
                if dn*dl'>=0&&atan(dl(3)/(dl(1)^2+dl(2)^2)^0.5)>=60/180*pi %判断不会被遮挡和在可视角范围内

                    fnum=fnum+1;
                end
            end
         end
        end
        if num>0

W=[W;xn,yn,num/((zj(end)-zj(1)+1)*(zi(end)-zi(1)+1)),fnum/((zj(end)-zj(1)+1)*(zi(end)-zi(1)+1))];

        end
    end

end
end

Wfind.m

%区域标记 
clc ,clear all;

zz=xlsread('区域高程数据.xlsx');
zz=zz';
x=(0:1:size(zz,2)-1)*38.2/1000; % 在 x 轴上取点
y=(0:1:size(zz,1)-1)*38.2/1000; % 在 y 轴上取点 
layers = [0,1500,3000,4200,5000]; 
[xx,yy]=meshgrid(x, y);
[c,h]=contourf(xx, yy, zz,layers);
colormap(flipud(gray)); %高程由低到高 呈现白到黑 flipud是反转色

set(h,'ShowText','off','LevelList',[0,3000,7000])%设定等高线的值 
axis equal
hold on

%%

%找出每个 P 点周围的方格
Lc=2*500*1.732/1000;

R=max(x(end)/2,y(end)/2);

W=Wfind(x(end)/2,y(end)/2,R,Lc); save('W.mat','W');

%%

%判断方格是否符合条件 
load W.mat
Wf=W(find((W(:,4)>0.9)==1),:);
Wxy=Wf(:,[1,2]);
Wxy=Wxy(2:end,:); 
save('Wxy.mat','Wxy'); 
plot(Wxy(:,1),Wxy(:,2),'rs'); 
sum(Wf(:,4))/sum(W(:,3))

第二问生命探测区域

这里的W.mat范围好像错了所以导致实际的方格数量过少

clemente620

关注

2
点赞
踩
9

收藏

觉得还不错? 一键收藏
0
评论
[MCM] 2017研究生数学建模竞赛A题绘图

clc,clear all;filename = 'a.csv'; %%数据预处理z = csvread(filename); %读取附件1 区域高程数据y0=1:1:2913; %按38.2米为1个单位进行简化绘图x0=1:1:2775;[xx0,yy0]=meshgrid(x0,y0);xx=xx0';%依题意行列对应,进行转置yy=yy0'...
复制链接

扫一扫