移动通信基站布局规划分析

随着社会的快速发展，我们的生活越来越离不开手机，手机已成为我们生活的必需品。大家有没有遇到这样的问题？就是在高楼大厦林立的地方，我们会发现手机信号特别差，通话断断续续并伴有杂音，视频卡成狗，地下车库更是盲区，为什么会出现这样的问题呢？

其实手机信号的传输是一种无线电的发射与接收过程，传播路径主要是以直线传播为主，若在传播过程中遇到阻碍（即与高楼大厦墙壁发生多次反射、绕射等物理现象。其中，绕射是指在电磁波传播路径上，当电波被尺寸较大(与波长相比)的障碍物遮挡时，电磁波改变传播方向的现象），就会大大削弱其信号强度。那么，该如何解决呢？

在现实生活中，运营商主要通过安装基站的方式增强手机信号。那么问题来了，如何最优选择基站的位置将成为重点考虑的问题？

目前，比较有代表性的就是射线跟踪模型。射线跟踪是一种被广泛用于移动通信中的预测无线电波传播特性的技术，由于移动通信中使用的超高频微波和光同属电磁波，有一定近似性，按光学方法辨认出多路径信道中收、发射机间所有主要的传播路径。一旦这些传播路径被辨认后，就可根据电波传播理论来计算每条传播路径信号的幅度、相位、延迟和极化，然后结合天线方向图和系统带宽就可得到到达接收点的所有传播路径的相干合成结果。

对于城市基站的二维模型，建筑群可被划分为一定的“块”，建筑物(即下图中带有灰色阴影的多边形)则被定义为“多边形”，多边形的“边”代表建筑物的表面，多边形的“顶点”则代表了建筑物的拐角。这种简化了的市区平面图大致反映出城市的主体结构，利用它进行射线跟踪，可以得到较为准确的路径损耗。

在多边形的顶点上仅能产生绕射，而在多边形的边上仅能产生反射，这些多次的反射、绕射及其组合便是收、发射机间的传播路径。二维射线跟踪模型可以通过以下两种规律分别确定反射传播路径和绕射传播路径：

(1)反射传播路径，如下图(a)所示，产生反射时入射角 等于反射角 ;

(2)绕射传播路径，如下图(b)所示，不论入射线以任意角度入射到建筑物顶点上，绕射射线都会以任意出射角向没有建筑物覆盖的区域传播。

现在我们要求出，在发射机Tx 即坐标为(500, 200)、接收机Rx 即坐标为(250, 350)之间，通过多次反射与绕射的组合，从而形成的主要传播路径。

拿到这样的题目，我们通过“枚举法”，即发射机Tx，从0°逐步加上步长（比如0.1°，甚至可以更小，步长越小精度越高，但是程序运行时间更长），从这个角度发出一条射线，然后考虑是反射还是绕射情况，直到最后的射线与接收机Rx非常接近（可以利用点到直线的距离公式求得），此次循环才能结束，进入下一个循环，最终当Tx角度超过360°，程序才能结束。

下面，我们通过Matlab编程，运行主程序main.m文件，由于该程序涉及到大量矩阵求逆运算等过程，计算量大，在笔者电脑上面运行该程序大约需要15分钟（取决于电脑配置），执行过程中请大家耐心等待哈。比如，点击运行程序后，可以先去洗个澡，回来后就能看到运行结果啦

最后，我们得到运行结果为：

详细代码如下：

（1）主程序main.m详细代码

clc;
clear;
close all;


figure
building={ [5.54892205638474e+002    1.53797468354430e+002;
            4.03648424543947e+002    1.53797468354430e+002;
            4.03648424543947e+002    1.00632911392405e+002;
            5.54892205638474e+002    1.00632911392405e+002;
            5.54892205638474e+002    1.53797468354430e+002];

            [5.54892205638474e+002    1.75738396624473e+002;
            5.54892205638474e+002    1.90928270042194e+002;
            5.29684908789386e+002    1.90928270042194e+002;
            5.29684908789386e+002    1.75738396624473e+002;
            5.54892205638474e+002    1.75738396624473e+002];

            [5.54892205638474e+002    1.95147679324895e+002;
            5.54892205638474e+002    2.59282700421941e+002;
            5.36981757877280e+002    2.59282700421941e+002;
            5.36981757877280e+002    1.95147679324895e+002;
            5.54892205638474e+002    1.95147679324895e+002];

            [4.42786069651741e+002    1.75738396624473e+002;
            4.42786069651741e+002    2.59282700421941e+002;
            4.03648424543947e+002    2.59282700421941e+002;
            4.03648424543947e+002    1.75738396624473e+002;
            4.42786069651741e+002    1.75738396624473e+002];

            [5.54892205638474e+002    2.81223628691983e+002;
            5.54892205638474e+002    3.00632911392405e+002;
            4.99170812603648e+002    3.00632911392405e+002;      %%%%%%淇敼
            4.99170812603648e+002    2.81223628691983e+002;
            5.54892205638474e+002    2.81223628691983e+002];

            [4.89220563847430e+002    2.81223628691983e+002;
            4.89220563847430e+002    2.95569620253165e+002;
            4.56053067993366e+002    2.95569620253165e+002;
            4.56053067993366e+002    2.81223628691983e+002;
            4.89220563847430e+002    2.81223628691983e+002];

            [4.32172470978441e+002    2.81223628691983e+002;
            4.32172470978441e+002    3.39451476793249e+002;
            4.03648424543947e+002    3.39451476793249e+002;
            4.03648424543947e+002    2.81223628691983e+002;
            4.32172470978441e+002    2.81223628691983e+002];

            [5.54892205638474e+002    3.09071729957806e+002;
            5.54892205638474e+002    3.39451476793249e+002;
            4.47429519071310e+002    3.39451476793249e+002;
            4.47429519071310e+002    3.09071729957806e+002;
            5.54892205638474e+002    3.09071729957806e+002];

            [4.75953565505804e+002    3.62236286919831e+002;
            4.75953565505804e+002    3.79113924050633e+002;
            4.71310116086236e+002    3.79113924050633e+002;
            4.71310116086236e+002    3.62236286919831e+002;
            4.75953565505804e+002    3.62236286919831e+002];

            [4.62686567164179e+002    3.62236286919831e+002;
            4.62686567164179e+002    4.18776371308017e+002;
            4.01658374792703e+002    4.18776371308017e+002;
            4.01658374792703e+002    3.62236286919831e+002;
            4.62686567164179e+002    3.62236286919831e+002];

            [4.99834162520730e+002    3.88396624472574e+002;
            4.99834162520730e+002    4.18776371308017e+002;
            4.80597014925373e+002    4.18776371308017e+002;
            4.80597014925373e+002    3.88396624472574e+002;
            4.99834162520730e+002    3.88396624472574e+002];

            [3.79104477611940e+002    4.18776371308017e+002;
            2.94195688225539e+002    4.18776371308017e+002;
            2.94195688225539e+002    3.62236286919831e+002;
            3.79104477611940e+002    3.62236286919831e+002;
            3.79104477611940e+002    4.18776371308017e+002];

            [3.79104477611940e+002    2.81223628691983e+002;
            3.79104477611940e+002    3.39451476793249e+002;
            2.12603648424544e+002    3.39451476793249e+002;
            2.12603648424544e+002    2.81223628691983e+002;
            3.79104477611940e+002    2.81223628691983e+002];

            [3.79104477611940e+002    1.53797468354430e+002;
            3.51243781094527e+002    1.53797468354430e+002;
            3.51243781094527e+002    1.00632911392405e+002;
            3.79104477611940e+002    1.00632911392405e+002;
            3.79104477611940e+002    1.53797468354430e+002];

            [3.79104477611940e+002    1.75738396624473e+002;
            3.79104477611940e+002    2.59282700421941e+002;
            3.32669983416252e+002    2.59282700421941e+002;
            3.32669983416252e+002    2.16244725738397e+002;
            3.57213930348259e+002    2.16244725738397e+002;
            3.57213930348259e+002    1.75738396624473e+002;
            3.79104477611940e+002    1.75738396624473e+002];

            [3.38640132669983e+002    1.53797468354430e+002;
            2.41791044776119e+002    1.53797468354430e+002;
            2.41791044776119e+002    1.25949367088608e+002;
            3.02819237147595e+002    1.25949367088608e+002;
            3.02819237147595e+002    1.00632911392405e+002;
            3.27363184079602e+002    1.00632911392405e+002;
            3.27363184079602e+002    1.25949367088608e+002;
            3.38640132669983e+002    1.25949367088608e+002;
            3.38640132669983e+002    1.53797468354430e+002];

            [2.90215588723051e+002    1.75738396624473e+002;
            2.90215588723051e+002    2.08649789029536e+002;
            2.65008291873964e+002    2.08649789029536e+002;
            2.65008291873964e+002    2.59282700421941e+002;
            2.12603648424544e+002    2.59282700421941e+002;
            2.12603648424544e+002    2.33966244725738e+002;
            2.43117744610282e+002    2.33966244725738e+002;
            2.43117744610282e+002    2.01054852320675e+002;
            2.12603648424544e+002    2.01054852320675e+002;
            2.12603648424544e+002    1.75738396624473e+002;
            2.90215588723051e+002    1.75738396624473e+002];

            };
N=17;


for i=1:N
    plot(building{i}(:,1),building{i}(:,2));
    fill(building{i}(:,1),building{i}(:,2),[0.5 0.5 0.5],'FaceAlpha',0.8);
    hold on
end
axis equal
xlabel('X(m)');
ylabel('Y(m)');
%grid on


TX=[500;200];
plot(TX(1),TX(2),'*r');
h_text1 = text(TX(1)-6,TX(2)-10,['Tx'],'color','r');
set(h_text1,'fontsize',14)    %设置字体大小


RX=[250;350];
plot(RX(1),RX(2),'or');     
h_text2 = text(RX(1)-5,RX(2)+10,['Rx'],'color','r');
set(h_text2,'fontsize',14)    %设置字体大小    


[aa]=calculate_point(TX(1),TX(2),building);
num1=size(aa);
num=num1(1);
for i=1:num
    [bb]=calculate_point(aa(i,1),aa(i,2),building);
    num2=size(bb);
    number=num2(1);
    for j=1:number
         for thea=0:0.1:2*pi
            x1=bb(j,1);y1=bb(j,2);x2=x1+cos(thea);y2=y1+sin(thea); 
            node=calculate_node(x1,y1,building);
            if x2<=max(building{node}(:,1))&&x2 >=min(building{node}(:,1)) && y2<=max(building{node}(:,2))&&y2 >=min(building{node}(:,2))
                continue
            end
               [insect1,insect2]=insection_point(x1,y1,x2,y2,building);
               if insect1(1)==0&&insect1(2)==0
                   continue
               end
               value1=calculate(insect1(1),insect1(2),RX(1),RX(2),building);
               zhi=([RX(1),RX(2)]-[insect1(1),insect1(2)])*([insect2(1),insect2(2)]-[insect1(1),insect1(2)])'/norm(([RX(1),RX(2)]-[insect1(1),insect1(2)]))/norm([insect2(1),insect2(2)]-[insect1(1),insect1(2)]);
               if zhi >=0.95&&zhi<=1&&value1==1
                   plot([TX(1),aa(i,1)],[TX(2),aa(i,2)],'k');
                   plot([aa(i,1),bb(j,1)],[aa(i,2),bb(j,2)],'k');
                   plot([bb(j,1),insect1(1)],[bb(j,2),insect1(2)],'k');
                   plot([insect1(1),RX(1)],[insect1(2),RX(2)],'k');
               end
        end
    end     
end


for thea=0:0.1:2*pi
    x1=TX(1);y1=TX(2);x2=x1+cos(thea);y2=y1+sin(thea);     
    [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
    if insect1(1)==0&&insect1(2)==0
        continue
    end
    x1=insect1(1);y1=insect1(2);
    x2=insect2(1);y2=insect2(2);    
    [x,y]=line_to_point(x1,y1,x2,y2,building);
    if x==0&&y==0
        continue
    end
    [x1,y1]=point_in(x,y,insect3(1),insect3(2),a,b,c,d);
    if x1==0&&y1==0
        continue
    end
   [bb]=calculate_point(x,y,building);
    num2=size(bb);
    number=num2(1);
    for j=1:number     
        value1=calculate(bb(j,1),bb(j,2),RX(1),RX(2),building);
        if value1==1
            plot([TX(1),x1],[TX(2),y1],'r');
            plot([x1,x],[y1,y],'r');
            plot([x,bb(j,1)],[y,bb(j,2)],'r');
            plot([bb(j,1),RX(1)],[bb(j,2),RX(2)],'r');
        end
    end  
end        


[aa]=calculate_point(TX(1),TX(2),building);
num1=size(aa);
num=num1(1);
for i=1:num
    for thea=0:0.1:2*pi
        x1=aa(i,1);y1=aa(i,2);x2=x1+cos(thea);y2=y1+sin(thea);     
        node=calculate_node(x1,y1,building);
        if x2<=max(building{node}(:,1))&&x2 >=min(building{node}(:,1)) && y2<=max(building{node}(:,2))&&y2 >=min(building{node}(:,2))
            continue
        end
        [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
        if insect1(1)==0&&insect1(2)==0
            continue
        end
        x1=insect1(1);y1=insect1(2);
        x2=insect2(1);y2=insect2(2);    
        [x,y]=line_to_point(x1,y1,x2,y2,building);
        if x==0&&y==0
            continue
        end
        [x1,y1]=point_in(x,y,insect3(1),insect3(2),a,b,c,d);
        if x1==0&&y1==0
            continue
        end  
        value1=calculate(x,y,RX(1),RX(2),building);
        if value1==1
            plot([TX(1),aa(i,1)],[TX(2),aa(i,2)],'r');
            plot([aa(i,1),x1],[aa(i,2),y1],'r');
            plot([x1,x],[y1,y],'r');
            plot([x,RX(1)],[y,RX(2)],'r');
        end
    end  
end

[aa]=calculate_point(TX(1),TX(2),building);
num1=size(aa);
num=num1(1);
for i=1:num
    [bb]=calculate_point(aa(i,1),aa(i,2),building); 
    num2=size(bb);
    number=num2(1);
    for j=1:number
        value1=calculate(bb(j,1),bb(j,2),RX(1),RX(2),building);
        if value1==1
            plot([TX(1),aa(i,1)],[TX(2),aa(i,2)],'g');
            plot([aa(i,1),bb(j,1)],[aa(i,2),bb(j,2)],'g');
            plot([bb(j,1),RX(1)],[bb(j,2),RX(2)],'g');
        end
    end
end      


[aa]=calculate_point(TX(1),TX(2),building);
num1=size(aa);
num=num1(1);
for i=1:num  
    for thea=0:0.1:2*pi
         for k=1:4
            x1=aa(i,1);y1=aa(i,2);
            x2=x1+cos(thea);y2=y1+sin(thea);   
            node=calculate_node(x1,y1,building);
            if x2<=max(building{node}(:,1)) && x2 >=min(building{node}(:,1)) && y2<=max(building{node}(:,2))&&y2 >=min(building{node}(:,2))
                continue
            end
            for j=1:k
               [insect1,insect2]=insection_point(x1,y1,x2,y2,building);
               if insect1(1)==0&&insect1(2)==0
                   break
               end
               xm{j}=insect1;ym{j}=insect2;
               x1=insect1(1);y1=insect1(2);
               x2=insect2(1);y2=insect2(2); 
            end
            value1=calculate(insect1(1),insect1(2),RX(1),RX(2),building);
            zhi=([RX(1),RX(2)]-[insect1(1),insect1(2)])*([insect2(1),insect2(2)]-[insect1(1),insect1(2)])'/norm(([RX(1),RX(2)]-[insect1(1),insect1(2)]))/norm([insect2(1),insect2(2)]-[insect1(1),insect1(2)]);
            if zhi >=0.95&&zhi<=1&&value1==1
               plot([TX(1),aa(i,1)],[TX(2),aa(i,2)],'k');
               xx=aa(i,1);yy=aa(i,2);
               for m=1:j
                  plot([xx,xm{m}(1)],[yy,xm{m}(2)],'k');
                  xx=xm{m}(1);yy=xm{m}(2);
               end
               plot([xx,RX(1)],[yy,RX(2)],'k');
           end            
         end
    end
end


 for thea=0:0.1:2*pi
    for thea1=0:0.1:2*pi
         for k=1:3

             x1=TX(1);y1=TX(2);x2=x1+cos(thea);y2=y1+sin(thea);     
            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2);    
            [x,y]=line_to_point(x1,y1,x2,y2,building);
            if x==0&&y==0
                continue
            end
            [x1,y1]=point_in(x,y,insect3(1),insect3(2),a,b,c,d);
            if x1==0&&y1==0
                continue
            end
            aa1=x1;bb1=y1;
            aa2=x;bb2=y;
            x2=x+cos(thea1);y2=y+sin(thea1);   
            node=calculate_node(x,y,building);
            if x2<=max(building{node}(:,1)) && x2 >=min(building{node}(:,1)) && y2<=max(building{node}(:,2))&&y2 >=min(building{node}(:,2))
                continue
            end
            for j=1:k
               [insect1,insect2]=insection_point(x,y,x2,y2,building);
               if insect1(1)==0&&insect1(2)==0
                   break
               end
               xm{j}=insect1;ym{j}=insect2;
               x=insect1(1);y=insect1(2);
               x2=insect2(1);y2=insect2(2); 
            end
            value1=calculate(insect1(1),insect1(2),RX(1),RX(2),building);
            zhi=([RX(1),RX(2)]-[insect1(1),insect1(2)])*([insect2(1),insect2(2)]-[insect1(1),insect1(2)])'/norm(([RX(1),RX(2)]-[insect1(1),insect1(2)]))/norm([insect2(1),insect2(2)]-[insect1(1),insect1(2)]);
            if zhi >=0.95&&zhi<=1&&value1==1
               plot([TX(1),aa1],[TX(2),bb1],'k');
               plot([aa1,aa2],[bb1,bb2],'k');
               xx=aa2;yy=bb2;
               for m=1:j
                  plot([xx,xm{m}(1)],[yy,xm{m}(2)],'k');
                  xx=xm{m}(1);yy=xm{m}(2);
               end
               plot([xx,RX(1)],[yy,RX(2)],'k');
           end            
         end
    end    
 end           


 for thea=0:0.1:2*pi
    for thea1=0:0.1:2*pi
         for k=1:2

             x1=TX(1);y1=TX(2);x2=x1+cos(thea);y2=y1+sin(thea);     
            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2);  
            aa3=x1;bb3=y1;

            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2); 

            [x,y]=line_to_point(x1,y1,x2,y2,building);
            if x==0&&y==0
                continue
            end
            [x1,y1]=point_in(x,y,insect3(1),insect3(2),a,b,c,d);
            if x1==0&&y1==0
                continue
            end
            aa1=x1;bb1=y1;
            aa2=x;bb2=y;
            x2=x+cos(thea1);y2=y+sin(thea1);   
            node=calculate_node(x,y,building);
            if x2<=max(building{node}(:,1)) && x2 >=min(building{node}(:,1)) && y2<=max(building{node}(:,2))&&y2 >=min(building{node}(:,2))
                continue
            end
            for j=1:k
               [insect1,insect2]=insection_point(x,y,x2,y2,building);
               if insect1(1)==0&&insect1(2)==0
                   break
               end
               xm{j}=insect1;ym{j}=insect2;
               x=insect1(1);y=insect1(2);
               x2=insect2(1);y2=insect2(2); 
            end
            value1=calculate(insect1(1),insect1(2),RX(1),RX(2),building);
            zhi=([RX(1),RX(2)]-[insect1(1),insect1(2)])*([insect2(1),insect2(2)]-[insect1(1),insect1(2)])'/norm(([RX(1),RX(2)]-[insect1(1),insect1(2)]))/norm([insect2(1),insect2(2)]-[insect1(1),insect1(2)]);
            if zhi >=0.95&&zhi<=1&&value1==1
               plot([TX(1),aa3],[TX(2),bb3],'k');
               plot([aa3,aa1],[bb3,bb1],'k');
               plot([aa1,aa2],[bb1,bb2],'k');
               xx=aa2;yy=bb2;
               for m=1:j
                  plot([xx,xm{m}(1)],[yy,xm{m}(2)],'k');
                  xx=xm{m}(1);yy=xm{m}(2);
               end
               plot([xx,RX(1)],[yy,RX(2)],'k');
           end            
         end
    end    
 end 
 
 for thea=0:0.1:2*pi
    for thea1=0:0.1:2*pi
         for k=1:1

             x1=TX(1);y1=TX(2);x2=x1+cos(thea);y2=y1+sin(thea);     
            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2);  
            aa4=x1;bb4=y1;

            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2); 
            aa3=x1;bb3=y1;  

            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2); 

            [x,y]=line_to_point(x1,y1,x2,y2,building);
            if x==0&&y==0
                continue
            end
            [x1,y1]=point_in(x,y,insect3(1),insect3(2),a,b,c,d);
            if x1==0&&y1==0
                continue
            end
            aa1=x1;bb1=y1;
            aa2=x;bb2=y;
            x2=x+cos(thea1);y2=y+sin(thea1);   
            node=calculate_node(x,y,building);
            if x2<=max(building{node}(:,1)) && x2 >=min(building{node}(:,1)) && y2<=max(building{node}(:,2))&&y2 >=min(building{node}(:,2))
                continue
            end
            for j=1:k
               [insect1,insect2]=insection_point(x,y,x2,y2,building);
               if insect1(1)==0&&insect1(2)==0
                   break
               end
               xm{j}=insect1;ym{j}=insect2;
               x=insect1(1);y=insect1(2);
               x2=insect2(1);y2=insect2(2); 
            end
            value1=calculate(insect1(1),insect1(2),RX(1),RX(2),building);
            zhi=([RX(1),RX(2)]-[insect1(1),insect1(2)])*([insect2(1),insect2(2)]-[insect1(1),insect1(2)])'/norm(([RX(1),RX(2)]-[insect1(1),insect1(2)]))/norm([insect2(1),insect2(2)]-[insect1(1),insect1(2)]);
            if zhi >=0.95&&zhi<=1&&value1==1
               plot([TX(1),aa4],[TX(2),bb4],'k');
               plot([aa4,aa3],[bb4,bb3],'k');
               plot([aa3,aa1],[bb3,bb1],'k');
               plot([aa1,aa2],[bb1,bb2],'k');
               xx=aa2;yy=bb2;
               for m=1:j
                  plot([xx,xm{m}(1)],[yy,xm{m}(2)],'k');
                  xx=xm{m}(1);yy=xm{m}(2);
               end
               plot([xx,RX(1)],[yy,RX(2)],'k');
           end            
         end
    end    
 end         


 for thea=0:0.1:2*pi
    for thea1=0:0.1:2*pi
             x1=TX(1);y1=TX(2);x2=x1+cos(thea);y2=y1+sin(thea);     
            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2);  
            aa5=x1;bb5=y1;

            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2);  
            aa4=x1;bb4=y1;

            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2); 
            aa3=x1;bb3=y1;

            [insect1,insect2,insect3,a,b,c,d]=insection_point(x1,y1,x2,y2,building);
            if insect1(1)==0&&insect1(2)==0
                continue
            end
            x1=insect1(1);y1=insect1(2);
            x2=insect2(1);y2=insect2(2); 

            [x,y]=line_to_point(x1,y1,x2,y2,building);
            if x==0&&y==0
                continue
            end
            [x1,y1]=point_in(x,y,insect3(1),insect3(2),a,b,c,d);
            if x1==0&&y1==0
                continue
            end
            aa2=x1;bb2=y1;
            aa1=x;bb1=y;
            value1=calculate(x,y,RX(1),RX(2),building);
            if value1==1
                plot([TX(1),aa5],[TX(2),bb5],'r');
                plot([aa5,aa4],[bb5,bb4],'r');
                plot([aa4,aa3],[bb4,bb3],'r');
                plot([aa3,aa2],[bb3,bb2],'r');
                plot([aa2,aa1],[bb2,bb1],'r');
                plot([aa1,RX(1)],[bb1,RX(2)],'r');
            end
    end      
 end

（2）calculate.m详细代码

function value=calculate(x2,y2,x3,y3,build)
    N=17;
    k=1;mm=0;
    value=0;
    while k >=1
        row_column1=size(build{k});
        row1=row_column1(1);
        for p=1:row1-1
            X1=[x3,y3];Y1=[x2,y2];X2=build{k}(p,:);Y2=build{k}(p+1,:);


            if X1(1)==Y1(1)
                X=X1(1);
                k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
                b2=X2(2)-k2*X2(1); 
                Y=k2*X+b2;
            end
            if X2(1)==Y2(1)
                X=X2(1);
                k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
                b1=X1(2)-k1*X1(1);
                Y=k1*X+b1;
            end
            if X1(1)~=Y1(1)&&X2(1)~=Y2(1)
                k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
                k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
                b1=X1(2)-k1*X1(1);
                b2=X2(2)-k2*X2(1);
                if k1==k2
                    continue
                else
                    X=(b2-b1)/(k1-k2);
                    Y=k1*X+b1;
                end
            end


            Q=([X,Y]-X1)*(Y1-[X,Y])';
            W=([X,Y]-X2)*(Y2-[X,Y])';
            if (Q >0.001)&&(W >0.001)
                mm=1;
            end
        end 
        if mm==1
            mm=0;
            k=1;
            break
        end    
        k=k+1;     
        if k >N
            break
        end
    end
    if k >N
        value=1;
        k=1;
    end  
   end

（3）calculate_node.m详细代码

function node=calculate_node(x1,y1,build)
    N=17;m=0;node=0;
    for k=1:N
        row_column1=size(build{k});
        row1=row_column1(1);
        for p=1:row1-1
            x2=build{k}(p,1);y2=build{k}(p,2);
         if abs(x1-x2)<=0.1&&abs(y1-y2)<=0.1
             node=k;
             m=1;
             break
         end
         if m==1
             break 
         end
        end  
    end
end

（4）calculate_point.m详细代码

function [x1]=calculate_point(x2,y2,build)
    N=17;
    k=1;m=1;mm=0;dd=0;
    node=calculate_node(x2,y2,build);
    for i=1:N
        if i==node
            continue
        end
        row_column=size(build{i});
        row=row_column(1);

        value=findway(x2,y2,build);
        if value==i
            continue
        end

        for j=1:row-1
            x=build{i}(j,1);
            y=build{i}(j,2);
            while k >=1
                row_column1=size(build{k});
                row1=row_column1(1);
                for p=1:row1-1
                    X1=[x,y];Y1=[x2,y2];X2=build{k}(p,:);Y2=build{k}(p+1,:);

                    if (X1(1)==Y1(1)&&X1(2)==Y1(2))  ||  (X2(1)==Y2(1)&&X2(2)==Y2(2))
                        continue
                    end
                    if X1(1)==Y1(1)
                        X=X1(1);
                        k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
                        b2=X2(2)-k2*X2(1); 
                        Y=k2*X+b2;
                    end
                    if X2(1)==Y2(1)
                        X=X2(1);
                        k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
                        b1=X1(2)-k1*X1(1);
                        Y=k1*X+b1;
                    end
                    if X1(1)~=Y1(1)&&X2(1)~=Y2(1)
                        k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
                        k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
                        b1=X1(2)-k1*X1(1);
                        b2=X2(2)-k2*X2(1);
                        if k1==k2
                            continue
                        else
                            X=(b2-b1)/(k1-k2);
                            Y=k1*X+b1;
                        end
                    end

                    Q=([X,Y]-X1)*(Y1-[X,Y])';
                    W=([X,Y]-X2)*(Y2-[X,Y])';
                    if (Q >0.001)&&(W >0.001)
                        mm=1;
                    end
                end 
                if mm==1
                    mm=0;
                    k=1;
                    break
                end    
                k=k+1;     
                if k >N
                    break
                end
            end
            if k >N
                for ii=1:N
                    row_column11=size(build{ii});
                    row11=row_column11(1);
                    for pp=1:row11-1
                        x22=build{ii}(pp,1);y22=build{ii}(pp,2);
                        QQ=([x22,y22]-Y1)*([x,y]-[x22,y22])'/norm([x22,y22]-Y1)/norm([x,y]-[x22,y22]);
                        if QQ==1
                            dd=1;
                            break
                        end
                    end
                    if dd==1
                            break
                    end
                end
                if dd==0
                    x1(m,1)=x;x1(m,2)=y;
                    m=m+1;
                end
                k=1;
                dd=0;
            end          
        end
    end
end

（5）findway.m详细代码

function value=findway(xx,yy,build)
    N=17;
    value=0;
    for i=1:N
        row_column=size(build{i});
        row=row_column(1);
        for j=1:row-1
            x=build{i}(j,1);
            y=build{i}(j,2);
            if x==xx&&y==yy
                value=i;
                break
            end
        end
        if value==i
            break
        end
    end
end

（6）insection_point.m详细代码

function [point1,point2,point3,a,b,c,d]=insection_point(x2,y2,x3,y3,build)
    N=17;min_num=10^10;
    point1(1)=0;point1(2)=0;point2(1)=0;point2(2)=0;point3(1)=0;point3(2)=0;a=0;b=0;c=0;d=0;
    for k=1:N
        row_column1=size(build{k});
        row1=row_column1(1);
        for p=1:row1-1
            X1=[x3,y3];Y1=[x2,y2];X2=build{k}(p,:);Y2=build{k}(p+1,:);
            if X1(1)==Y1(1)
                X=X1(1);
                k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
                b2=X2(2)-k2*X2(1); 
                Y=k2*X+b2;
            end
            if X2(1)==Y2(1)
                X=X2(1);
                k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
                b1=X1(2)-k1*X1(1);
                Y=k1*X+b1;
            end
            if X1(1)~=Y1(1)&&X2(1)~=Y2(1)
                k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
                k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
                b1=X1(2)-k1*X1(1);
                b2=X2(2)-k2*X2(1);
                if k1==k2
                    continue
                else
                    X=(b2-b1)/(k1-k2);
                    Y=k1*X+b1;
                end
            end


            Q=(X1-Y1)*([X,Y]-Y1)';
            W=([X,Y]-X2)*(Y2-[X,Y])';
            if (Q >0.001)&&(W >0.001)
                num1=(x2-X)^2+(y2-Y)^2;
               if min_num >=num1
                   min_num=num1;
                   point1=[X,Y];
                   a=build{k}(p,1);b=build{k}(p,2);
                   c=build{k}(p+1,1);d=build{k}(p+1,2);
                       if a==c
                            point2=[2*a-x2,y2];
                            point3=point2;
                            point2=2.*point1-point2;
                       end
                       if b==d
                            point2=[x2,2*d-y2];
                            point3=point2;
                            point2=2.*point1-point2;
                       end
               end              
            end            
        end  
    end
end

（7）line_to_point.m详细代码

function [x,y]=line_to_point(x1,y1,x2,y2,build)
    x=0;y=0;
    N=17;
    for k=1:N
        row_column1=size(build{k});
        row1=row_column1(1);
        for p=1:row1-1
            x3=build{k}(p,1);y3=build{k}(p,2);
            value1=calculate(x1,y1,x3,y3,build);
            zhi=([x3,y3]-[x1,y1])*([x2,y2]-[x1,y1])'/norm([x3,y3]-[x1,y1])/norm([x2,y2]-[x1,y1]);
            if zhi >=0.99&&zhi<=1&&value1==1
                x=x3;y=y3;
            end
        end
    end
end

（8）point_in.m详细代码

function [value1,value2]=point_in(x2,y2,x3,y3,a,b,c,d)
    value1=0;value2=0;
    X1=[x3,y3];Y1=[x2,y2];X2=[a,b];Y2=[c,d];
    if X1(1)==Y1(1)
        X=X1(1);
        k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
        b2=X2(2)-k2*X2(1); 
        Y=k2*X+b2;
    end
    if X2(1)==Y2(1)
        X=X2(1);
        k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
        b1=X1(2)-k1*X1(1);
        Y=k1*X+b1;
    end
    if X1(1)~=Y1(1)&&X2(1)~=Y2(1)
        k1=(Y1(2)-X1(2))/(Y1(1)-X1(1));
        k2=(Y2(2)-X2(2))/(Y2(1)-X2(1));
        b1=X1(2)-k1*X1(1);
        b2=X2(2)-k2*X2(1);
        if k1==k2
            value1=0;value2=0;
        else
            X=(b2-b1)/(k1-k2);
            Y=k1*X+b1;
        end
    end
    Q=([X,Y]-X1)*(Y1-[X,Y])';
    W=([X,Y]-X2)*(Y2-[X,Y])';
    if (Q >0.001)&&(W >0.001)
       value1=X;value2=Y;
    end
end