【无人机】基于最小二乘法实现无线电地图重构附matlab代码

✅作者简介：热爱科研的Matlab仿真开发者，修心和技术同步精进，matlab项目合作可私信。

🍎个人主页：Matlab科研工作室

🍊个人信条：格物致知。

更多Matlab仿真内容点击👇

智能优化算法       神经网络预测       雷达通信       无线传感器        电力系统

信号处理              图像处理               路径规划       元胞自动机        无人机 

❤️ 内容介绍

在当今数字化时代，我们生活在一个无线连接的世界中。无线通信技术的快速发展使得我们能够随时随地进行无线通信，无论是通过手机、电视还是互联网。然而，随着无线设备的不断增多，无线信号的传播环境也变得越来越复杂。为了更好地理解和利用这个复杂的无线环境，无线电地图学习和重建成为了一个备受关注的研究领域。

无线电地图学习和重建是指利用接收到的无线信号来推断和重建无线信号的传播路径和环境特征。这种技术可以帮助我们更好地了解无线信号在不同环境中的传播规律，从而提高无线通信的性能和可靠性。在过去的几十年中，许多研究者已经提出了各种各样的方法和算法来实现无线电地图学习和重建。

其中一种常用的方法是利用环境语义来辅助无线电地图的学习和重建。环境语义是指利用环境中的物体、结构和地理信息来推断无线信号的传播规律。例如，我们可以利用建筑物的结构和材料来推测无线信号在城市环境中的传播路径，或者利用地形和植被信息来推断无线信号在农村或森林中的传播规律。通过利用环境语义，我们可以更准确地重建无线信号的传播路径和环境特征，从而提高无线通信系统的性能。

在实际应用中，无线电地图学习和重建可以应用于许多领域。例如，在智能交通系统中，我们可以利用无线电地图学习和重建来推断车辆之间的距离和速度，从而实现更智能化的交通管理。在无线传感器网络中，无线电地图学习和重建可以帮助我们更好地了解传感器节点之间的通信质量，从而优化网络拓扑结构和传输协议。此外，无线电地图学习和重建还可以应用于室内定位、无线电频谱管理等领域。

然而，无线电地图学习和重建也面临着许多挑战。首先，由于无线信号的传播路径和环境特征受到许多因素的影响，如多径效应、信号衰减和噪声干扰，因此准确地学习和重建无线电地图是一项复杂的任务。其次，由于环境语义的复杂性和多样性，如何有效地利用环境语义来辅助无线电地图学习和重建也是一个具有挑战性的问题。此外，由于无线信号的频谱资源有限，如何在保证学习和重建准确性的同时降低无线信号采样的成本也是一个重要的研究方向。

总的来说，无线电地图学习和重建是一个具有挑战性但又备受关注的研究领域。通过利用环境语义，我们可以更好地了解和利用复杂的无线环境，从而提高无线通信的性能和可靠性。随着无线通信技术的不断发展和应用需求的不断增加，无线电地图学习和重建将会在未来发挥更加重要的作用。

无线电地图可以用于链路性能预测、无线中继规划和源定位。本文构建了一个空中到地面的无线电地图，用于预测连接地面终端和低空无人机之间的每个链路的信道增益。挑战在于无线电地图的测量样本不足，每个数据点是6维的，因为发射器和接收器都有三个空间自由度。此外，将大尺寸的无线电地图传达或共享给移动决策者也是昂贵的。经典方法，如k最近邻（KNN）和Kriging，在数据不足时可能会失败。本文提出利用环境的几何语义中的传播特性来辅助无线电地图的构建。具体而言，通过重建虚拟几何环境来构建无线电地图。

• 开发了一个嵌入环境感知多度信道模型的多类虚拟障碍物模型。

• 制定了一个最小二乘问题来学习虚拟障碍物地图和模型参数。

本文研究了最小二乘问题的部分凸性，并基于此，开发了一种高效的无线电地图学习算法。此外，采用数据驱动方法构建残余阴影地图，进一步提高了构建的无线电地图的细节。我们的数值结果证实，与Kriging基准相比，所提出的方法显著提高了预测准确性，并将所需测量量减少了一半以上。当将构建的无线电地图应用于基于接收信号强度（RSS）的定位时，在密集城市环境中观察到了显著的性能改进，达到了不到20米的准确度。

🔥核心代码

​set(0,'defaultfigurecolor',[1 1 1])clear, close alladdpath Sub%rng('default')​%% DATA COLLECTION --------------------------------------------------------% Topology parametersK_class = 2;sample_rate = 1;meter_pixel = 9;map_height = 50;    % Drone height for air-to-ground power mapsample_height = 50; % Drone height for learning (Learning stage)neighbour_mat = [0 0; -1 0; 0 -1; 1 0; 0 1] * 39; %neighbour_weight = [0.6; 0.1; 0.1; 0.1; 0.1]; % Sum-to-onetest_positions = round(1:100);noise_guass = 3; % dB stdn_ue = 50; % [8*8 9*10 11*11 13*13 15*16 18*18 21*22 25*26 30*30 35*36]*8n_uav = 50 * length(sample_height);    % 15:50:9 50:200:5 100:400n_mea = n_ue * n_uav;n_uav = round(n_uav / length(sample_height));residual_on = false;residsec_on = false;residthi_on = false;residfou_on = false;residfiv_on = true;segment_on = false;segresid_on = false;plot_on = false;kmeans_on = false;​% DATA = load('radiomap_shanghai100tx_2.5GHz.mat');% DATA = load('radiomap_simulated100tx_noise.mat');% DATA = load('radiomap_shanghai115tx_28GHz.mat');% DATA = load('radiomap_simulated100tx_nonoise');DATA = load('radiomap_simulated100tx_3class.mat');pos_ue = DATA.PosUE;Xdata = zeros(n_mea, 6); Gdata = zeros(n_mea, 1);Ddata = zeros(n_mea, 1);Gtrue = cell(n_ue, 1);cnt = 0;​figure(1)title("Users' locations")hold onpos_ue_all = zeros(size(pos_ue));userids = sort(randperm(length(pos_ue), n_ue));for id = 1:length(userids)    % i must increase    i = userids(id);    % Training data    for j = 1:length(sample_height)        uav_height = sample_height(j);        x = pos_ue(i,1); y = pos_ue(i,2); z = pos_ue(i,3);        I = (DATA.RadioMap(:, 6) == uav_height) ...            & (DATA.RadioMap(:, 1) == x) ...            & (DATA.RadioMap(:, 2) == y)...            & (DATA.RadioMap(:, 3) == z);        Rm2D = DATA.RadioMap(I, :);        Xvec = Rm2D(:, 4); Yvec = Rm2D(:, 5); Zvec = Rm2D(:, end);        [Xmat,Ymat] = meshgrid((min(Xvec):5:max(Xvec)),(min(Yvec):5:max(Yvec)));        Zmat = griddata(Xvec, Yvec, Zvec, Xmat, Ymat);        position = [x-min(Xvec)+1 y-min(Yvec)+1 z];        Xmat = Xmat - min(Xvec) + 1;        Ymat = Ymat - min(Yvec) + 1;        % Data pre-processing        Zmat(isnan(Zmat)) = min(min(Zmat));        Zmat = Zmat + randn(size(Zmat)) * noise_guass;        Xvec = Xmat(:);        Yvec = Ymat(:);        Zvec = Zmat(:);        dist = log10(vecnorm([Xvec Yvec uav_height*ones(length(Xvec),1)] - ...            ones(length(Xvec),3).*position, 2, 2));        position = [ceil(position(1:2)) z];        if ~exist('idx', 'var')            %idx = 1:ceil(length(Xvec)/nDrone):length(Xvec);            idx = randperm(length(Xvec), n_uav);        end        uav_positon = ceil([Xvec(idx) Yvec(idx)]);        uav_num = length(uav_positon);        cnt = cnt + uav_num;        Xdata(cnt-uav_num+1:cnt, :) = [uav_positon ...            uav_height.*ones(uav_num,1) ones(uav_num,3).*position];        Gdata(cnt-uav_num+1:cnt, :) = Zvec(idx);        Ddata(cnt-uav_num+1:cnt, :) = dist(idx);​        figure(1)        plot(position(1), position(2), '.');        text(position(1), position(2), string(i));​        pos_ue_all(i, :) = position;    end    % Test data    uav_height = map_height;    x = pos_ue(i,1); y = pos_ue(i,2); z = pos_ue(i,3);    I = (DATA.RadioMap(:, 6) == uav_height) ...        & (DATA.RadioMap(:, 1) == x) ...        & (DATA.RadioMap(:, 2) == y)...        & (DATA.RadioMap(:, 3) == z);    Rm2D = DATA.RadioMap(I, :);    Xvec = Rm2D(:, 4); Yvec = Rm2D(:, 5); Zvec = Rm2D(:, end);    [Xmat,Ymat] = meshgrid((min(Xvec):5:max(Xvec)),(min(Yvec):5:max(Yvec)));    Zmat = griddata(Xvec, Yvec, Zvec, Xmat, Ymat);    position = [x-min(Xvec)+1 y-min(Yvec)+1 z];    Xmat = Xmat - min(Xvec) + 1;    Ymat = Ymat - min(Yvec) + 1;    % Data pre-processing    Xvec = Xmat(:);    Yvec = Ymat(:);    Zvec = Zmat(:);    Zmat(isnan(Zmat)) = min(Zvec);    Gtrue{i} = Zmat;endlenx = ceil(max(Xvec)/meter_pixel);leny = ceil(max(Yvec)/meter_pixel);n_mea = cnt;Xdata = Xdata(1:cnt, :); Ydata = Gdata(1:cnt, :);Ddata = Ddata(1:cnt, :);maps.BldMapZ = zeros(lenx, leny);maps.BldPosMat = ones(lenx, leny);maps.FolPosMat = zeros(lenx, leny);maps.meterPerPixel = meter_pixel;​% Generate collinear obstacles setobstacles = (1:lenx*leny)';nObst = length(obstacles);cols = cell(nObst, 2);   % collinear measurements of each obstaclecovB = cell(n_mea, 2);for id = 1:n_mea    uav_pos_meter = Xdata(id, 1:3);    uav_pos_pixel = [floor(uav_pos_meter(1:2)/meter_pixel)+1, uav_pos_meter(3)];    ue_pos_meter = Xdata(id, 4:6);    ue_pos_pixel = [floor(ue_pos_meter(1:2)/meter_pixel)+1, ue_pos_meter(3)];    [cov_bld, cov_z] = covBldZ(uav_pos_pixel, ue_pos_pixel, lenx);    cov_bld = cov_bld(cov_z > 0 & cov_z <= map_height);    cov_z = cov_z(cov_z > 0 & cov_z <= map_height);    covB{id, 1} = cov_bld;    covB{id, 2} = cov_z;    for i = 1:length(cov_bld)        j = cov_bld(i);        ib = find(obstacles == j, 1, 'first');        if ~isempty(ib)            if ~isempty(cols{ib, 1})                cols{ib, 1} = [cols{ib, 1} id];                cols{ib, 2} = [cols{ib, 2} cov_z(i)];            else                cols{ib, 1} = id;                cols{ib, 2} = cov_z(i);            end        end    endend​%% LOCAL POLY RECONSTRUCTION ----------------------------------------------MAXLOOP = 9;tolerance = 1e-3;metrics = cell(MAXLOOP, 6);maps.droneHeightMap = map_height;maps.neighbourMat = neighbour_mat;maps.neighbourWeight = neighbour_weight;​% InitializationR.X = Xdata;P = polyfit(Ddata, Ydata, 1);delta = Ydata - (P(1) * Ddata + P(2));w = delta > 0;W = [w w ~w ~w];A = [Ddata ones(n_mea, 1)];A_ = W .* [A A];X = (A_' * A_) \ A_' * Ydata;A = repmat(A, 1, K_class+1);R.Z1 = zeros(n_mea, K_class+1);R.Alpha = zeros(1, K_class+1);R.Beta = zeros(1, K_class+1);for k = 1:K_class+1    R.Alpha(k) = X(1) - (X(1) - X(3)) * ((k-1)/K_class);    R.Beta(k) = X(2) - (X(2) - X(4)) * ((k-1)/K_class);    R.Z1(:, k) = -abs(Ydata-(R.Alpha(k)*Ddata+R.Beta(k)));end​% -- coarse resolution --dsfactor = 4;       % Downsampling factor[small_maps, S, Hs] = downsampleMaps(maps, dsfactor, R);R.Hs = Hs; R.S = S;smap_hat = ones([K_class numel(small_maps.BldMapZ)]) * map_height;​for i = 1:MAXLOOP    [smap_hat, W] = optimizeH(R, small_maps, smap_hat);    A_ = W .* A;    X = (A_' * A_) \ A_' * Ydata;    R.Alpha = X(1:2:end);    R.Beta = X(2:2:end);    for k = 1:K_class+1        R.Z1(:, k) = -abs(Ydata-(R.Alpha(k)*Ddata+R.Beta(k)));    endend​% -- medium resolution --dsfactor2 = 2;  [medium_maps, S, Hs] = downsampleMaps(maps, dsfactor2, R);R.Hs = Hs; R.S = S;mmap_hat = zeros([K_class numel(medium_maps.BldMapZ)]);for k = 1:K_class    upsample_matrix = repelem(reshape(smap_hat(k, :), ...        size(small_maps.BldMapZ)), round(dsfactor/dsfactor2), ...        round(dsfactor/dsfactor2));    mmap_hat(k, :) = vec(upsample_matrix(1:size(medium_maps.BldMapZ, 1),...        1:size(medium_maps.BldMapZ, 2)));end​for i = 1:MAXLOOP    [mmap_hat, W] = optimizeH(R, medium_maps, mmap_hat);    A_ = W .* A;    X = (A_' * A_) \ A_' * Ydata;    R.Alpha = X(1:2:end);    R.Beta = X(2:2:end);    for k = 1:K_class+1        R.Z1(:, k) = -abs(Ydata-(R.Alpha(k)*Ddata+R.Beta(k)));    endend​% -- fine resolution (3 meter) --dsfactor3 = 1;[maps, S, Hs] = downsampleMaps(maps, dsfactor3, R);R.Hs = Hs; R.S = S;map_hat = zeros([K_class numel(maps.BldMapZ)]);for k = 1:K_class    upsample_matrix = repelem(reshape(mmap_hat(k, :), ...        size(medium_maps.BldMapZ)), round(dsfactor2/dsfactor3), ...        round(dsfactor2/dsfactor3));    map_hat(k, :) = vec(upsample_matrix(1:size(maps.BldMapZ, 1),...        1:size(maps.BldMapZ, 2)));end​h  = waitbar(0, 'Estimate building heights');for i = 1:MAXLOOP    waitbar(i / MAXLOOP, h);    [map_hat, W] = optimizeH(R, maps, map_hat);    A_ = W .* A;    X = (A_' * A_) \ A_' * Ydata;    R.Alpha = X(1:2:end);    R.Beta = X(2:2:end);    for k = 1:K_class+1        R.Z1(:, k) = -abs(Ydata-(R.Alpha(k)*Ddata+R.Beta(k)));    endendclose(h);​% hf = showmap(BldMapHat + FolMapHat, Maps.meterPerPixel, 8);% title('Esitmated building and foliage height');​if kmeans_on == true    c_i = D1 ./ abs((X(1)*Ddata+X(2)) - (X(3)*Ddata+X(4)));    obs_clu = zeros(nObst, map_height);    for i = 1:map_height        for j = 1:nObst            obs_clu(j, i) = mean(c_i(cols{j, 1}(ceil(cols{j, 2}) == i)));        end    end    obs_clu(isnan(obs_clu)) = 0;    bld_ind = zeros(lenx, leny);    bld_ind(maps.BldPosMat > 0) = kmeans(obs_clu, 5);end​figure, hold onfor k = 1:K_class+1    plot(Ddata(W(:, 2*k) > 1/(K_class+1)), ...        Ydata(W(:, 2*k) > 1/(K_class+1)), '.', 'Color', rand([1 3]));endfigure, hold onfor k = K_class+1:-1:1    plot(Ddata(W(:, 2*k) > 1/(K_class+1)), ...        Ydata(W(:, 2*k) > 1/(K_class+1)), '.', 'Color', rand([1 3]));endif isfield(DATA, 'BldPosMat')    k = K_class;    map_hat(k, DATA.BldPosMat(round(1:DATA.lenX/lenx:DATA.lenX), ...        round(1:DATA.lenY/leny:DATA.lenY)) < 1) = 0;endheights = reshape(map_hat(K_class, :), size(maps.BldMapZ));for i = 1:length(pos_ue_all)    x = round(pos_ue_all(i, 1) / meter_pixel);    y = round(pos_ue_all(i, 2) / meter_pixel);    if x == 0 || y == 0, continue, end    heights(x, y) = 0;end% heights(heights >= droneHeightMap | heights == droneHeightMap/2) = 0;% indicator = heights >= UserHeight;% heights = medfilt1(heights .* (heights <= DroneHeight)) .* indicator;figure; showmap(heights, meter_pixel);title(sprintf('Height from L model at %d meters', round(map_height)));​uav_num = uav_num * length(sample_height);​%% RADIOMAP AND PERFORMANCE -----------------------------------------------% Channel model C.A1 = X(1); C.B1 = X(2); C.S1 = 0;C.A2 = 0; C.B2 = 0; C.S2 = 0; C.A3 = X(3); C.B3 = X(4); C.S3 = 0;channel_model.C = C;channel_model.los_nlos_trans = 0;channel_model.noise = 0;​metrics_mse = []; metrics_mae = [];XYkr = [vec(meshgrid(1:lenx, 1:leny)'), vec(meshgrid(1:leny, 1:lenx)), ...    ones(lenx * leny, 1) * map_height / meter_pixel];Xann = [vec(meshgrid(1:lenx, 1:leny)'), vec(meshgrid(1:leny, 1:lenx)) ...    * meter_pixel, ones(lenx * leny, 1) * map_height];net = fitnet([32 32 16 16 8 8 4 4]); %fitnet([(4:8).^2 (8:-1:2).^2]);net = train(net, [Xdata Ddata X(1)*Ddata+X(2) X(3)*Ddata+X(4)]', Ydata', ...    'useParallel', 'no', 'showResources', 'no');​j = 1;for i = 1:length(pos_ue)    if pos_ue_all(i, 1) == 0 || pos_ue_all(i, 2) == 0, continue, end    if ~ismember(i, test_positions), j = j + 1; continue, end    Gtru = Gtrue{i}';    % Radiomap reconstruction    Xdata_i = [Xdata((j-1)*uav_num+1:j*uav_num,1:2)/meter_pixel, ...        Xdata((j-1)*uav_num+1:j*uav_num,3)/meter_pixel];    Ydata_i = Ydata((j-1)*uav_num+1:j*uav_num);    idx = randperm(length(Xdata_i), round(length(Xdata_i) / sample_rate));    Xdata_i = Xdata_i(idx, :); Ydata_i = Ydata_i(idx, :);    d_variogram = variogram(Xdata_i(:, 1:2), Ydata_i);    [~, ~, ~, vstruct] = variogramfit(d_variogram.distance, d_variogram.val, ...        [], [], [], 'model', 'exponential', 'plotit', false);    if n_mea > 1e9        idx = randperm(length(Xdata_i), round(length(Xdata_i) / sample_rate));        Gkri = kriging(vstruct,Xdata_i(idx,:),false,Ydata_i(idx),XYkr,false);    else        idx = randperm(length(Xdata), round(length(Xdata) / sample_rate));        Gkri = kriging(vstruct, Xdata(idx,:)/meter_pixel, false, Ydata(idx), ...        [XYkr ones(length(XYkr), 3).*pos_ue_all(i, :)/meter_pixel], false);    end    Xann_upos = ones(length(XYkr), 3).*pos_ue_all(i, :);    Xann_dist = log10(vecnorm(Xann - Xann_upos, 2, 2));    Xann_features = [Xann Xann_upos ...        Xann_dist X(1)*Xann_dist+X(2) X(3)*Xann_dist+X(4)]';    Gann = net(Xann_features);    Gann = reshape(Gann, [lenx leny]);    Gkri = reshape(Gkri, [lenx leny]);    Gknn = powerMapReconKNN(R,Ydata,pos_ue_all(i, :),map_height,maps);    Gbld = powerMapReconBld(R,heights,heights,pos_ue_all(i, :),map_height,maps);    Gbdk = powerMapReconKBld(R, map_hat, pos_ue_all(i, :), map_height, maps);    % Radiomap reconstruction with residual    if residual_on == true    [~, id] = min(abs(sample_height - map_height));    droneHeightRes = sample_height(id)/meter_pixel;    ResX = floor(Xdata_i(:, 1:2)) + 1;    ResX = ResX(Xdata_i(:, 3) == droneHeightRes, :);    ResG = Ydata_i(Xdata_i(:, 3) == droneHeightRes);    Gbld = powerMapReconBld(R,heights,heights,pos_ue_all(i, :),map_height,maps);    Gbld = imgaussfilt(Gbld, 1);    for k = 1:length(ResX)        ResG(k) = ResG(k) - Gbld(ResX(k, 1), ResX(k, 2));    end    d_variogram = variogram(ResX, ResG);    [~, ~, ~, vstruct] = variogramfit(d_variogram.distance, d_variogram.val, ...        [], [], [], 'model', 'exponential', 'plotit', false);    Rkri = kriging(vstruct, ResX(:,1), ResX(:,2), ...        ResG, meshgrid(1:lenx, 1:leny)', meshgrid(1:leny, 1:lenx));    Gbld = Gbld + reshape(Rkri, [lenx leny]);    end    if residsec_on == true    ResX = Xdata / meter_pixel;    P = (X(1)*Ddata + X(2) - Ydata)./(X(1)*Ddata + X(2) - X(3)*Ddata - X(4));    P(P < 0) = 0; P(P > 1) = 1;    ResG = Ydata - w1.*(X(1)*Ddata + X(2)) - w2.*(X(3)*Ddata + X(4));    d_variogram = variogram(Xdata((j-1)*uav_num+1:j*uav_num,1:2)/meter_pixel, ...        Ydata((j-1)*uav_num+1:j*uav_num));    [~, ~, ~, vstruct] = variogramfit(d_variogram.distance, d_variogram.val, ...        [], [], [], 'model', 'exponential', 'plotit', false);    ResRec = kriging(vstruct, ResX, false, ResG, ...    [XYkr ones(length(XYkr), 3).*pos_ue_all(i, :)/meter_pixel], false);    ResRec = reshape(ResRec, [lenx leny]);    Gbld = Gbld + ResRec;    end    if residthi_on == true    Rind = true(length(ResX), 1);    for k = 1:length(ResX)        ResG(k) = ResG(k) - Gbld(ResX(k, 1), ResX(k, 2));        Rind(k) = Gind(ResX(k, 1), ResX(k, 2));    end    d_variogram = variogram(ResX(Rind,:), ResG(Rind));    [~, ~, ~, vstruct] = variogramfit(d_variogram.distance, d_variogram.val, ...        [], [], [], 'model', 'exponential', 'plotit', false);    LOSSeg = kriging(vstruct, ResX(Rind,1), ResX(Rind,2), ...        ResG(Rind), meshgrid(1:lenx, 1:leny)', meshgrid(1:leny, 1:lenx));    d_variogram = variogram(ResX(~Rind,:), ResG(~Rind));    [~, ~, ~, vstruct] = variogramfit(d_variogram.distance, d_variogram.val, ...        [], [], [], 'model', 'exponential', 'plotit', false);    NLOSSeg = kriging(vstruct, ResX(~Rind,1), ResX(~Rind,2), ...        ResG(~Rind), meshgrid(1:lenx, 1:leny)', meshgrid(1:leny, 1:lenx));    Gbld = Gbld + LOSSeg .* Gind + NLOSSeg .* ~Gind;    end    if residfou_on == true    ResY = []; ResR = [];    for id = 1:length(sample_height)        h = sample_height(id)/meter_pixel;        ResX = [floor(Xdata_i(:, 1:2)) + 1 Xdata_i(:, 3)];        ResX = ResX(Xdata_i(:, 3) == h, :);        ResG = Ydata_i(Xdata_i(:, 3) == h);        Gbld = powerMapReconBld(R,heights,heights,pos_ue_all(i, :),h,maps);        Gbld = imgaussfilt(Gbld, 3);        for k = 1:length(ResX)            ResG(k) = ResG(k) - Gbld(ResX(k, 1), ResX(k, 2));        end        x = ResX(:, 1) - pos_ue_all(i, 1)/meter_pixel;        y = ResX(:, 2) - pos_ue_all(i, 2)/meter_pixel;        h = h/meter_pixel - pos_ue_all(i, 3)/meter_pixel;        l = vecnorm([x y], 2, 2);        theta1 = atan2(y,x); theta2 = atan2(h,l);        ResR = [ResR; min(theta1+pi,pi-theta1)/pi*180 ...            min(theta2+pi,pi-theta2)/pi*180 sqrt(l.^2 + h^2)];        ResY = [ResY; ResG];    end    d_variogram = variogram(ResR, ResY);    [~, ~, ~, vstruct] = variogramfit(d_variogram.distance, d_variogram.val, ...        [], [], [], 'model', 'exponential', 'plotit', false);    Rkri = krigingR(vstruct, Xdata((j-1)*uav_num+1:j*uav_num,1:3)/ ...        meter_pixel, ResY, XYkr, pos_ue_all(i, :)/meter_pixel);    Gbld = powerMapReconBld(R,heights,heights,pos_ue_all(i, :),map_height,maps);    Gbld = imgaussfilt(Gbld, 3) + reshape(Rkri, [lenx leny]);    end    if residfiv_on == true    [~, id] = min(abs(sample_height - map_height));    droneHeightRes = sample_height(id)/meter_pixel;    ResX = floor(Xdata_i(:, 1:2)) + 1;    ResX = ResX(Xdata_i(:, 3) == droneHeightRes, :);    ResG = Ydata_i(Xdata_i(:, 3) == droneHeightRes);    Gind = powMapRecBldSof(R,heights,[],pos_ue_all(i, :),map_height,maps,[]);    Gind = imgaussfilt(Gind, 3);    Glos = powMapRecBldSof(R,heights-realmax,[],pos_ue_all(i, :),map_height,maps);    Gnlo = powMapRecBldSof(R,heights+realmax,[],pos_ue_all(i, :),map_height,maps);    Gbld = Gind .* Glos + (1 - Gind) .* Gnlo;    Gbdn = Gbld;    for k = 1:length(ResX)        ResG(k) = ResG(k) - Gbld(min(lenx, ResX(k, 1)), min(leny, ResX(k, 2)));    end    d_variogram = variogram(ResX, ResG);    [~, ~, ~, vstruct] = variogramfit(d_variogram.distance, d_variogram.val, ...        [], [], [], 'model', 'exponential', 'plotit', false);    Rkri = kriging(vstruct, ResX(:,1), ResX(:,2), ...        ResG, meshgrid(1:lenx, 1:leny)', meshgrid(1:leny, 1:lenx));    Gbld = Gbld + reshape(Rkri, [lenx leny]);    end    % Radiomap reconstruction with segmentation    if segment_on == true    d_variogram = variogram(Xdata(w1, :)/meter_pixel, Ydata(w1));    [~, ~, ~, vstruct] = variogramfit(d_variogram.distance, d_variogram.val, ...        [], [], [], 'model', 'exponential', 'plotit', false);    LOSSeg = kriging(vstruct,Xdata(w1, :)/meter_pixel, false, Ydata(w1), ...        [XYkr ones(length(XYkr), 3).*pos_ue_all(i, :)/meter_pixel], false);    d_variogram = variogram(Xdata(w2, :)/meter_pixel, Ydata(w2));    [~, ~, ~, vstruct] = variogramfit(d_variogram.distance, d_variogram.val, ...        [], [], [], 'model', 'exponential', 'plotit', false);    NLOSSeg = kriging(vstruct,Xdata(w2, :)/meter_pixel, false, Ydata(w2), ...        [XYkr ones(length(XYkr), 3).*pos_ue_all(i, :)/meter_pixel], false);    LOSSeg = reshape(LOSSeg, [lenx leny]);    NLOSSeg = reshape(NLOSSeg, [lenx leny]);    Gind = powerMapReconBld(R,heights,heights,pos_ue_all(i, :),map_height,maps,1);    Gbld = LOSSeg .* Gind + NLOSSeg .* ~Gind;    end    % Radiomap reconstruction with residual and segmentation    if segresid_on == true    Xlos = Xdata(w1, :) / meter_pixel;     Ylos = Ydata(w1) - Ddata(w1) * X(1) - X(2);    Xnlos = Xdata(w2, :) / meter_pixel;    Ynlos = Ydata(w2) - Ddata(w2) * X(3) - X(4);    d_variogram = variogram(Xlos, Ylos);    [~, ~, ~, vstruct] = variogramfit(d_variogram.distance, d_variogram.val, ...        [], [], [], 'model', 'exponential', 'plotit', false);    LOSSeg = kriging(vstruct, Xlos, false, Ylos, ...        [XYkr ones(length(XYkr), 3).*pos_ue_all(i, :)/meter_pixel], false);    d_variogram = variogram(Xnlos, Ynlos);    [~, ~, ~, vstruct] = variogramfit(d_variogram.distance, d_variogram.val, ...        [], [], [], 'model', 'exponential', 'plotit', false);    NLOSSeg = kriging(vstruct, Xnlos, false, Ynlos, ...        [XYkr ones(length(XYkr), 3).*pos_ue_all(i, :)/meter_pixel], false);    LOSSeg = reshape(LOSSeg, [lenx leny]);    NLOSSeg = reshape(NLOSSeg, [lenx leny]);    Gind = powerMapReconBld(R,heights,heights,pos_ue_all(i, :),map_height,maps,1);    Gbld = Gbld + LOSSeg .* Gind + NLOSSeg .* ~Gind;    end    % Performance and figures    xm = length(Zmat);    if xm <= lenx        Gbld = Gbld(round(1:lenx/xm:lenx),round(1:lenx/xm:lenx));        Gkri = Gkri(round(1:lenx/xm:lenx),round(1:lenx/xm:lenx));        Gknn = Gknn(round(1:lenx/xm:lenx),round(1:lenx/xm:lenx));        Gann = Gann(round(1:lenx/xm:lenx),round(1:lenx/xm:lenx));        Gbdk = Gbdk(round(1:lenx/xm:lenx),round(1:lenx/xm:lenx));    else        Gtru = Gtru(round(1:xm/lenx:xm),round(1:xm/lenx:xm));    end    if plot_on == true    figure(300 + i);    subplot(2, 3, 1); showmap(Gkri, meter_pixel);    subplot(2, 3, 2); showmap(Gknn, meter_pixel);    subplot(2, 3, 3); showmap(Gann, meter_pixel);    subplot(2, 3, 4); showmap(Gbdk, meter_pixel);    subplot(2, 3, 5); showmap(Gbld, meter_pixel);    subplot(2, 3, 6); showmap(Gtru, meter_pixel);    % figure(400 + i); showmap(reshape(Rkri, [lenX lenY]), 9);    disp_str = strcat('Position ',string(i),': ',string(mean(abs(vec(...        Gtru-Gkri)))),' vs. ',string(mean(abs(vec(Gtru-Gknn)))),' vs. ',...        string(mean(abs(vec(Gtru-Gbld)))),'; ',string(mse(Gtru-Gkri)),...        ' vs. ',string(mse(Gtru-Gknn)),' vs. ',string(mse(Gtru-Gbld)));    disp(disp_str); input('');    end    metrics_mse = [metrics_mse; mse(Gtru-Gkri) mse(Gtru-Gknn) ...        mse(Gtru-Gann) mse(Gtru-Gbdk) mse(Gtru-Gbld)];    metrics_mae = [metrics_mae; mean(abs(vec(Gtru-Gkri))) ...        mean(abs(vec(Gtru-Gknn))) mean(abs(vec(Gtru-Gann))) ...        mean(abs(vec(Gtru-Gbdk))) mean(abs(vec(Gtru-Gbld)))];    j = j + 1;enddisp([mean(metrics_mae, 1) mean(metrics_mse, 1) n_ue]);​​​

❤️ 运行结果

⛄ 参考文献

⛳️ 代码获取关注我

❤️部分理论引用网络文献，若有侵权联系博主删除

❤️ 关注我领取海量matlab电子书和数学建模资料

🍅 仿真咨询

1 各类智能优化算法改进及应用

生产调度、经济调度、装配线调度、充电优化、车间调度、发车优化、水库调度、三维装箱、物流选址、货位优化、公交排班优化、充电桩布局优化、车间布局优化、集装箱船配载优化、水泵组合优化、解医疗资源分配优化、设施布局优化、可视域基站和无人机选址优化

2 机器学习和深度学习方面

卷积神经网络（CNN）、LSTM、支持向量机（SVM）、最小二乘支持向量机（LSSVM）、极限学习机（ELM）、核极限学习机（KELM）、BP、RBF、宽度学习、DBN、RF、RBF、DELM、XGBOOST、TCN实现风电预测、光伏预测、电池寿命预测、辐射源识别、交通流预测、负荷预测、股价预测、PM2.5浓度预测、电池健康状态预测、水体光学参数反演、NLOS信号识别、地铁停车精准预测、变压器故障诊断

2.图像处理方面

图像识别、图像分割、图像检测、图像隐藏、图像配准、图像拼接、图像融合、图像增强、图像压缩感知

3 路径规划方面

旅行商问题（TSP）、车辆路径问题（VRP、MVRP、CVRP、VRPTW等）、无人机三维路径规划、无人机协同、无人机编队、机器人路径规划、栅格地图路径规划、多式联运运输问题、车辆协同无人机路径规划、天线线性阵列分布优化、车间布局优化

4 无人机应用方面

无人机路径规划、无人机控制、无人机编队、无人机协同、无人机任务分配

、无人机安全通信轨迹在线优化

5 无线传感器定位及布局方面

传感器部署优化、通信协议优化、路由优化、目标定位优化、Dv-Hop定位优化、Leach协议优化、WSN覆盖优化、组播优化、RSSI定位优化

6 信号处理方面

信号识别、信号加密、信号去噪、信号增强、雷达信号处理、信号水印嵌入提取、肌电信号、脑电信号、信号配时优化

7 电力系统方面

微电网优化、无功优化、配电网重构、储能配置

8 元胞自动机方面

交通流 人群疏散 病毒扩散 晶体生长 火灾扩散

9 雷达方面

卡尔曼滤波跟踪、航迹关联、航迹融合、状态估计