中国像素地图最优路线分析

系统简介

本文对最优路线问题进行分析研究。问题1将中国像素地图计算出0,1矩阵，利用穷举法对问题进行求解，从图片的左上角开始寻找满足数值的点并赋值。问题2是在不考虑实际地表的复杂情况下，求解出A、B、C、D四点不经过境外区域和红色区域的彼此之间的最短路线。利用的是运筹最优化来进行分析求解。在两点之间寻找所要经过的点，在图中画出每个点的范围，用MATLAB穷举范围内的每个点计算距离和，求解出距离和的最小值。

关键词 最优路线 穷举法 运筹最优化 MATLAB

引言

由于交通网络的迅速发展，车辆增加的同时，车辆出行量提高，人们对路程长短、能源消耗也有所要求。运输、出行，都会对路线进行规划，考虑地理位置、成本效益等因素，提高出行便捷水平。在实际应用中，根据位置约束条件，计划出最优路线，提高出行速度，降低能源消耗。
本文就最优路线问题进行分析，研究的主要对象为A、B、C、D四个地点，在不考虑实际地理因素的条件下，不能经过吉林、天津、山东、陕西、青海、江西、广西的区域，计算出彼此之间的最短路线。研究两点之间最短路线时，考虑大致经过的路线，穷举中间点的范围，建立约束规划模型，得出最优的规划模型后，计算出每个最短路线的长度，并对模型结果进行比较和评价。
本文分析了以上内容，依据穷举法及约束规划理论知识，研究最优路线的约束规划模型，同时给出该模型的评价，最终得到合理的最优方案。

1.问题重述和分析

1.1问题重述
下图是在上生成的中国像素地图. 每个小正方形边长为6px小正方形之间没有间距(即间距为0px)
图 1.1
问题 1. 根据图或图, 试计算上面中国像素地图中共有多少个这样的小正方形? 要求写出这样的矩阵, 元素为0或1. 如果小方块在国境内, 即用代表, 否则用表示. 并且要求矩阵中没有列向量为零向量.
问题2. 现在有四个地点. 见图中四个黑色小方块. 彼此之间运输货物, 但是运输途中不能经过吉林、天津、山东、陕西、青海、江西、广西这七个省份/直辖市/自治区(下图中红色区域). 试绘出彼此之间的最短路线, 且分别求出经过多少个小正方形?
这里的所谓最短路线不考虑实际地表的复杂情况（公路、铁路、河流、桥梁、地球曲率等等），仅以所给的平面像素图为准。

图 1.2
1.2问题分析
问题1，用MATLAB进行图片处理，在水平方向和竖直方向上每隔6个点读取一次颜色的数值，将6×6的小方格变成一个像素点，根据图片上小方块的RGB的值，确定小方块属于境内区域还是境外区域，最后可以计算出0-1矩阵A。
问题2，以B到C的路线为例（如图1.3所示），利用运筹学最优化来处理，预先估计两点之间按照4条直线走，中间点分别为x,y,z。在图中画出x,y,z的一个范围，穷举范围内的每个点计算距离和，求距离和的最小值。

图1.3

2.基本假设与符号说明

2.1基本假设
（1）假设所有6×6的小正方形从图片的左上角开始计算；
（2）对于所在中国境内的小正方形面积大于30则认定属于中国境内；
（3）假设人工选择的点范围合理，则所计算出的最短路线长度就是题目中所要求的最短路线；
（4）在实际行走过程趋于离散状态，当前目标线路就近原则；
（5）假设路线总长度在原像素点的长度除以6，则是小方块数。
2.2符号说明
num 中国境内的小方块数
经过路线的转向点
路线长度
P 境外区域、国内可行驶区域、国内不可行驶区域上点的集合
PLP 计算两点所经过的小方块数

3.模型的建立与求解

3.1问题1
3.1.1问题求解
由题目可知，将给定的中国像素地图，分为6×6的小方块，计算出在中国境内有多少个这样的小方块，并写出境外小方块为0、境内小方块为1的矩阵。
首先，用MATLAB读入图片后，根据图片变量确定矩阵。通过函数对图片的进行读取和显示，因为RGB的值大小是0到255，根据图片上每个点的红、绿、蓝三个分量的值，设置阈值，筛选出境外小方块，赋值为0。对于境内区域的确定，由于边境部分并不是完整的6×6小方块（如图3.1所示），

图 3.1
所以利用for循环寻找在境内小方块，计算每个小方块的面积，假设若面积大于30则属于境内小方块。利用穷举法，筛选出面积大于30的小方块，则判定为境内小方块，并赋值为1。在穷举筛选境内小方块的时候，如果符合条件，则小方块数num=num+1。可以得出题目中所要求解出的矩阵A（158×213如图3.2所示）以及境内小方块数14401个。

图3.2
3.2问题2
3.2.1模型的建立
通过题中的给出的条件来决定约束条件即位置约束，不能经过吉林、天津、山东、陕西、青海、江西、广西，即图中的境内红色区域以及不经过境外区域。人工判断出两点之间大致走的最短路径，经过的点为，根据假设可得线路总长度，确定位置使得路径最优。通过决策变量、目标函数、约束条件来建立相关的约束规划方程，运用MATLAB进行方程求解最优路线。
在建模过程中，两点之间需要寻找的n个点。这寻找难度较大，所以假设所选定决策变量三点合理，与出发点、目的点所连的线路总长度就是最优值。
本题一共有六条线，B-A,B-C,B-D,C-D,C-A,A-D，分别确定三点的取值范围，在所确定的范围内穷举每两点的距离和，找出最短路径。通常情况下，求解两点之间距离时，会想着用距离公式。但本题是在像素图（即网格图）上找路径，线路放大并不是一条完整的直线，每个小方块趋于离散状态，所以路线就近原则选取小方块，线路就是小方块的对角线。
题目中需要计算的是两点之间所经过的小方块数,用函数PLP表示计算。先假设在直线方程中两点坐标为，通过斜率，寻找像素点，判断计算所经过的小方块数。当时，变化一个单位，则变化小于一个单位，若，，，对四舍五入取整，反之也是如此计算；当时，变化一个单位，则变化少于一个单位，若，，，对四舍五入取整，反之也是如此计算。对于求出的计算出数值，则可以得出小方块的数量，其中P表示三个区域的集合。
对于位置约束的红色区域，用集合P表示。本题中一共有三个区域，境外区域、境内可行驶区域以及境内不可行驶区域，根据RGB的值，决定图中每个小方块所在的区域，将这三个区域分别装进三个集合里，每个集合里的元素赋值为30000,1,10000。如果线路所经过的点值超过10000则判定为经过非法区域，不符合要求，需要重新寻找新的线路。
根据上述内容建立约束规划模型,求解最短路径。
3.2.2模型的求解
在B-A线路中，点B(368,512)，点A(1201,295)，经过第一点的范围为；经过第二点的范围为；经过第三点的范围为。

最终线路在原像素地图中小方块数为925,穷举出的三点坐标分别为(352,420)(738,417)(1045,215)，线路如下图3.3所示

图3.3
在B-C线路中，点B(368,512)，点C(738,529)，经过第一点的范围为；经过第二点的范围为；经过第三点的范围为。

最终线路在原像素地图中小方块数为445，穷举出的三点坐标分别为(355,445)(479,413)(592,453)，线路如下图3.4所示

图3.4
在B-D线路中，点B(368,512)，点D(919,783)，经过第一点的范围为；经过第二点的范围为；经过第三点的范围为。

最终线路在原像素地图中小方块数为886，穷举出的三点坐标分别为(340,517)(346,627)(985,665)，线路如下图3.5所示

图3.5

3.2.3模型的对比
在线路C-D,B-D中，可以选择江西的东西两面，在求解的过程中选择经过的是江西东面，因此需要比较在这两个方向上哪个路径最短。
如下图3.13所示，C-D从江西西面经过，所经过的第一二两点的范围不变，经过第三点的范围为。

最终线路在原像素地图中小方块数为392，穷举出的三点坐标分别为(684,539)(677,624)(874,828)，线路如下图3.9所示。

图3.9
在B-D线路中，从江西西面经过时，至少需要找到四个点，才可以不经过红色区域符合题目要求。经过第一点的范围为；经过第二点的范围为；经过第三点的范围为；经过第四点的范围为。

最终线路在原像素地图中小方块数为699，穷举出的四点坐标分别为(346,524)(346,605)(790,744)(874,828)。线路如下图3.10所示

图3.10
综上所述，B-A,B-C,B-D,C-D,C-A,A-D所经过的小正方形个数分别为155,75,117,92,127。
3.3模型的评价
该模型是求解最短路径，则利用运筹最优化来处理。确定所要经过的大致线路，确定三点范围，利用穷举法，在这范围内寻找最优点连成最短路径。优点是可以逐步优化，在一定条件下求解出最优值；缺点则是计算速度慢，在这范围内一个一个点的测试是否符合条件；同时所算出的最优解并非真正意义上的最优值，可能需要经过n个点，才可以得出最短路。故此有限穷举点范围，具有一定的局限性。

附录
程序1：

clc;   clear all;  close all;
data=imread('图片2.png');        % 读入图片
figure;
imshow(data);
s=data(:,:,1);
[row0 col0]=size(s);
A=ones(row0,col0);
A(find(data(:,:,1)>=0 & data(:,:,1)<55 & data(:,:,2)>111 &data(:,:,2)<122 & data(:,:,3)>150 & data(:,:,3)<186))=0;  % 设置阈值，筛选境外的像素
A(303:333,1224:1252)=0;         % 字母A位置为境外，手动去除
figure;imagesc(A);
B=zeros(ceil(row0/6),ceil(col0/6));         % 图像矩阵初始化（所求的A）
num=0;
for flag1=1:ceil(row0/6)-1              % 境内小方块寻找，6像素边长
    for flag2=1:ceil(col0/6)-1
        tem=A((flag1-1)*6+1:flag1*6,(flag2-1)*6+1:flag2*6);
        if(sum(sum(tem))>30)            % 小方块内国内像素大于30认为是国内
            B(flag1,flag2)=1;
            num=num+1;
        end
    end
end
for flag1=ceil(row0/6)                  % 边缘部分判断
    for flag2=1:ceil(col0/6)
        tem=A((flag1-1)*6+1:row0,(flag2-1)*6+1:flag2*6);
        if(sum(sum(tem))>30)
        B(flag1,flag2)=1;
        num=num+1;
        end
    end
end
figure(3);imagesc(B)
num
mark=[];
for flag=1:ceil(col0/6)         %去掉B中0列向量
    flag
    if(sum(B(:,flag))==0)
        mark=[mark flag];
    end
end
[row2 col2]=find(s>100 & s<135 & data(:,:,2)>118 &&data(:,:,2)<140);        
row_H=ceil(row2/6);         
col_H=ceil(col2/6);         
LOC1=[row_H col_H];
LOC=sort(LOC1);     
for flag_H=length(col_H):-1:2        
    if(LOC(flag_H,:)==LOC(flag_H-1,:))
        LOC(flag_H,:)=[];
    end
end
 
LOC_H=[];
for flag=1:length(LOC)
    tem=LOC1-kron(LOC(flag,:),ones(length(LOC1),1));
    tem=sum(tem,2);
    [row_tem col_tem]=find(tem==0 );
    length(col_tem);
    if(length(col_tem)>10)                     
        LOC_H=[LOC_H ;LOC(flag,:)];     
    end
end
d_col=length(find(mark<60));            %   去掉0列向量，得到最终结果
LOC_H(:,2)=LOC_H(:,2)-d_col
B(:,mark)=[];
figure;imagesc(B);

程序2：
%% 
data1=data(:,:,1);
data2=data(:,:,2);
data3=data(:,:,3); 
% line([Xb,x10,x20,x30,Xa],[Yb,y10,y20,y30,Ya])
X=zeros(6,5);
Y=zeros(6,5);
TL00=zeros(1,6);
figure(1)
for kk=1:6
    if kk==1
        % % B-A  3点穷举
        Xa=1201;Ya=295;
        Xb=368;Yb=512;
        X1=345:352;Y1=420:441;
        X2=738:746;Y2=417:417;
        X3=1045:1060;Y3=215:222;
    elseif kk==2
        % % B-C  3点穷举
        Xa=738;Ya=529;
        Xb=368;Yb=512;
        X1=349:355;Y1=445:454;
        X2=479:487;Y2=413:421;
        X3=592:602;Y3=453:463;
    elseif kk==3
        % % B-D  3点穷举
        Xb=368;Yb=512;
        Xa=919;Ya=783;
        X1=340:346;Y1=517:524;
        X2=340:346;Y2=627:633;
        X3=985:992;Y3=658:665;
    elseif kk==4
        % % C-D  3点穷举
        Xb=738;Yb=529;
        Xa=919;Ya=783;
        X1=677:684;Y1=534:540;
        X2=677:684;Y2=624:630;
        X3=871:899;Y3=828:840;
    elseif kk==5
        % %  C-A  3点穷举
        Xa=1201;Ya=295;
        Xb=738;Yb=529;
        X1=693:698;Y1=486:491;
        X2=910:910;Y2=308:323;
        X3=1051:1060;Y3=222:232;
    elseif kk==6
        
        % % A-D  3点穷举
        Xa=919;Ya=783;
        Xb=1201;Yb=295;
        X1=1019:1031;Y1=219:225;
        X2=830:845;Y2=548:551;
        X3=985:992;Y3=658:665;
    else
        disp('!!!!!')
    end
    TL0=90000;
    % 图片阈值处理
    PIC=data1;
    PIC(find((data1>1)&(data1<47)))=30000; %非中国
    PIC(find((data1>47)&(data1<100)))=1; %能行驶区域
    PIC(find((data1>100)&(data1<200)))=10000; %不能行驶区域
    PIC(find(data1==0))=1;
    for i1=1:length(X1)
        for j1=1:length(Y1)
            for i2=1:length(X2)
                for j2=1:length(Y2)
                    for i3=1:length(X3)
                        for j3=1:length(Y3)
                            Tx1=X1(i1);
                            Ty1=Y1(j1);
                            Tx2=X2(i2);
                            Ty2=Y2(j2);
                            Tx3=X3(i3);
                            Ty3=Y3(j3);
                            TL1=PLP(Xb,Yb,Tx1,Ty1,PIC);
                            TL2=PLP(Tx1,Ty1,Tx2,Ty2,PIC);
                            TL3=PLP(Tx2,Ty2,Tx3,Ty3,PIC);
                            TL4=PLP(Tx3,Ty3,Xa,Ya,PIC);
                            TL=TL1+TL2+TL3+TL4;
                            if TL<TL0
                                x10=Tx1;
                                x20=Tx2;
                                x30=Tx3;
                                y10=Ty1;
                                y20=Ty2;
                                y30=Ty3;
                                TL0=TL;
                            else
                            end
                            
                        end
                    end
                end
            end
            
        end
       
    end
    kk
    PIC=data1;
    PIC(find((data1>1)&(data1<47)))=300; %非中国
    PIC(find((data1>47)&(data1<100)))=1; %能行驶区域
    PIC(find((data1>100)&(data1<200)))=100; %不能行驶区域
    PIC(find(data1==0))=200;
   subplot(2,3,kk)
    imagesc(PIC)
    hold on
    line([Xb,x10,x20,x30,Xa],[Yb,y10,y20,y30,Ya])
    
    X(kk,:)=[Xb,x10,x20,x30,Xa];
    Y(kk,:)=[Yb,y10,y20,y30,Ya];
    TL00(1,kk)=TL0;
end
 
figure(2)
 imagesc(PIC)
 hold on
 for jjj=1:6
line(X(jjj,:),Y(jjj,:))
 end

程序3：

function A=PLP(X,Y,XX,YY,Map)
       k=(YY-Y)/(XX-X);
       
       if abs(k)<1  % 则 x变化一个单位 y变化小于一个单位 
           if X>XX
               TTx=XX:X;
               TTy0=YY+(TTx-XX)*k;
               TTy=round(TTy0);
           else
               TTx=X:XX;
               TTy0=Y+(TTx-X)*k;
               TTy=round(TTy0);  
           end
       else
           if Y>YY
               TTy=YY:Y;
               TTx0=XX+(TTy-YY)*1/k;
               TTx=round(TTx0);
           else
               TTy=Y:YY;
               TTx0=X+(TTy-Y)*1/k;
               TTx=round(TTx0);

结论

本文主要运用到数学建模的知识，利用穷举法和约束规划知识建立模型.
建立的模型主要回答了以下问题：
（1）问题1，矩阵A是158×213的0.1矩阵（如图3.2所示），在中国境内的小方块数14401个。
（2）问题2，B-A,B-C,B-D,C-D,C-A,A-D所经过的小正方形个数分别为155,75,117,92,127。
本文中，在假设成立的情况下，利用穷举法，再结合约束规划，求出最优路径。通过比较数值，得出合理的数据，增加了模型的准确性。
同时，本文建立的模型也存在一些缺陷。模型的可行性是在模型的假设基础上成立的，由于每两点之间的路径，选取的是三点和四点转向，算出的小方块数，并非真正意义上的最小值。其次是穷举法计算速度慢，时间复杂度为。

致谢

转眼间，在宿迁学院四年的学习生活即将结束。在这四年里，从初入学校的不习惯到现在即将毕业时的不舍，很庆幸这四年里遇到了很多帮助过我的老师同学。.在此，向一切给过我支持和帮助的老师和同学致以深深的谢意.
感谢我的导师周克元老师，在周老师精心指导下完成本篇论文，论文框架的构建、相关资料的供给、以及后续的修改等各阶段都提出的宝贵建议。
未来的挫折还有很多，但我永远不会忘记我在宿院学习到的精神品质，这些都是我一生难得的财富。
最后再次感谢论文写作期间，给予我帮忙和支持的每一位老师、同学。