[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'; 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']); 更高清
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']); 更高清
因为颜色比较混乱 可考虑 黑白灰 重新设定
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']); 更高清
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']); 更高清
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是反转色
%% 第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米以下的点位数
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
覆盖率
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))
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
%区域标记 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范围好像错了 所以导致实际的方格数量过少