基于虚拟几何环境的无人机地空信道预测方法附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)));endfor 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)));endh = 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);endfigure, 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]);
⛳️ 运行结果
