?作者簡(jiǎn)介：熱愛(ài)科研的Matlab仿真開(kāi)發(fā)者，修心和技術(shù)同步精進(jìn)，matlab項(xiàng)目合作可私信。

??個(gè)人主頁(yè)：Matlab科研工作室

??個(gè)人信條：格物致知。

更多Matlab完整代碼及仿真定制內(nèi)容點(diǎn)擊??

智能優(yōu)化算法 ? ? ? 神經(jīng)網(wǎng)絡(luò)預(yù)測(cè) ? ? ? 雷達(dá)通信 ? ? ?無(wú)線傳感器 ? ? ? ?電力系統(tǒng)

信號(hào)處理 ? ? ? ? ? ? ?圖像處理 ? ? ? ? ? ? ? 路徑規(guī)劃 ? ? ? 元胞自動(dòng)機(jī) ? ? ? ?無(wú)人機(jī)

?? 內(nèi)容介紹

無(wú)人機(jī)通信技術(shù)是近年來(lái)備受關(guān)注的領(lǐng)域，其在軍事、民用等領(lǐng)域都有廣泛的應(yīng)用。而在無(wú)人機(jī)通信中，地面終端與低空無(wú)人機(jī)之間的信道增益預(yù)測(cè)、無(wú)線中繼規(guī)劃和源定位等問(wèn)題是需要解決的。本文提出了一種基于虛擬幾何環(huán)境的無(wú)人機(jī)地空信道預(yù)測(cè)方法，旨在通過(guò)構(gòu)建一張地空無(wú)人機(jī)信道增益預(yù)測(cè)地圖，為無(wú)人機(jī)通信提供可靠的基礎(chǔ)。

然而，由于信道預(yù)測(cè)地圖數(shù)據(jù)點(diǎn)的維度較高，每個(gè)數(shù)據(jù)點(diǎn)都是6維的，其中發(fā)送器和接收器各有三個(gè)空間自由度，因此構(gòu)建信道預(yù)測(cè)地圖所需的測(cè)量樣本不足，同時(shí)將大型信道預(yù)測(cè)地圖傳輸或共享給移動(dòng)決策者也是一項(xiàng)昂貴的任務(wù)。傳統(tǒng)的方法，如k最近鄰算法（KNN）和Kriging算法，可能因數(shù)據(jù)不足而失敗。因此，本文提出了一種基于環(huán)境幾何語(yǔ)義的傳播特性的無(wú)人機(jī)地空信道預(yù)測(cè)方法，即通過(guò)重構(gòu)虛擬幾何環(huán)境來(lái)構(gòu)建信道預(yù)測(cè)地圖。

具體而言，本文開(kāi)發(fā)了一個(gè)嵌入環(huán)境感知的多類虛擬障礙模型和多度信道模型，制定了一個(gè)最小二乘問(wèn)題來(lái)學(xué)習(xí)虛擬障礙地圖和模型參數(shù)。本文研究了最小二乘問(wèn)題的部分?jǐn)M凸性，并基于此開(kāi)發(fā)了一種高效的無(wú)人機(jī)地空信道預(yù)測(cè)算法。此外，本文采用數(shù)據(jù)驅(qū)動(dòng)方法構(gòu)建了一個(gè)剩余陰影地圖，以進(jìn)一步提高構(gòu)建的信道預(yù)測(cè)地圖的細(xì)節(jié)。數(shù)值結(jié)果表明，與Kriging基線相比，本文提出的方法顯著提高了預(yù)測(cè)精度，并將所需測(cè)量次數(shù)減少了一半以上。當(dāng)構(gòu)建的信道預(yù)測(cè)地圖應(yīng)用于基于接收信號(hào)強(qiáng)度（RSS）的定位時(shí)，在密集城市環(huán)境下實(shí)現(xiàn)了亞20米的精度，取得了顯著的性能提升。

總之，本文提出的基于虛擬幾何環(huán)境的無(wú)人機(jī)地空信道預(yù)測(cè)方法在無(wú)人機(jī)通信中具有重要的應(yīng)用價(jià)值。其通過(guò)重構(gòu)虛擬幾何環(huán)境，構(gòu)建了一張地空無(wú)人機(jī)信道增益預(yù)測(cè)地圖，為無(wú)人機(jī)通信提供了可靠的基礎(chǔ)，同時(shí)也為無(wú)人機(jī)通信領(lǐng)域的研究提供了新的思路和方法。

?? 部分代碼

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]);

?? 運(yùn)行結(jié)果

?? 參考文獻(xiàn)

?? 部分理論引用網(wǎng)絡(luò)文獻(xiàn)，若有侵權(quán)聯(lián)系博主刪除

?? ?關(guān)注我領(lǐng)取海量matlab電子書(shū)和數(shù)學(xué)建模資料

?? ?私信完整代碼和數(shù)據(jù)獲取及論文數(shù)模仿真定制

1 各類智能優(yōu)化算法改進(jìn)及應(yīng)用

生產(chǎn)調(diào)度、經(jīng)濟(jì)調(diào)度、裝配線調(diào)度、充電優(yōu)化、車間調(diào)度、發(fā)車優(yōu)化、水庫(kù)調(diào)度、三維裝箱、物流選址、貨位優(yōu)化、公交排班優(yōu)化、充電樁布局優(yōu)化、車間布局優(yōu)化、集裝箱船配載優(yōu)化、水泵組合優(yōu)化、解醫(yī)療資源分配優(yōu)化、設(shè)施布局優(yōu)化、可視域基站和無(wú)人機(jī)選址優(yōu)化

2 機(jī)器學(xué)習(xí)和深度學(xué)習(xí)方面

卷積神經(jīng)網(wǎng)絡(luò)（CNN）、LSTM、支持向量機(jī)（SVM）、最小二乘支持向量機(jī)（LSSVM）、極限學(xué)習(xí)機(jī)（ELM）、核極限學(xué)習(xí)機(jī)（KELM）、BP、RBF、寬度學(xué)習(xí)、DBN、RF、RBF、DELM、XGBOOST、TCN實(shí)現(xiàn)風(fēng)電預(yù)測(cè)、光伏預(yù)測(cè)、電池壽命預(yù)測(cè)、輻射源識(shí)別、交通流預(yù)測(cè)、負(fù)荷預(yù)測(cè)、股價(jià)預(yù)測(cè)、PM2.5濃度預(yù)測(cè)、電池健康狀態(tài)預(yù)測(cè)、水體光學(xué)參數(shù)反演、NLOS信號(hào)識(shí)別、地鐵停車精準(zhǔn)預(yù)測(cè)、變壓器故障診斷

2.圖像處理方面

圖像識(shí)別、圖像分割、圖像檢測(cè)、圖像隱藏、圖像配準(zhǔn)、圖像拼接、圖像融合、圖像增強(qiáng)、圖像壓縮感知

3 路徑規(guī)劃方面

旅行商問(wèn)題（TSP）、車輛路徑問(wèn)題（VRP、MVRP、CVRP、VRPTW等）、無(wú)人機(jī)三維路徑規(guī)劃、無(wú)人機(jī)協(xié)同、無(wú)人機(jī)編隊(duì)、機(jī)器人路徑規(guī)劃、柵格地圖路徑規(guī)劃、多式聯(lián)運(yùn)運(yùn)輸問(wèn)題、車輛協(xié)同無(wú)人機(jī)路徑規(guī)劃、天線線性陣列分布優(yōu)化、車間布局優(yōu)化

4 無(wú)人機(jī)應(yīng)用方面

無(wú)人機(jī)路徑規(guī)劃、無(wú)人機(jī)控制、無(wú)人機(jī)編隊(duì)、無(wú)人機(jī)協(xié)同、無(wú)人機(jī)任務(wù)分配、無(wú)人機(jī)安全通信軌跡在線優(yōu)化

5 無(wú)線傳感器定位及布局方面

傳感器部署優(yōu)化、通信協(xié)議優(yōu)化、路由優(yōu)化、目標(biāo)定位優(yōu)化、Dv-Hop定位優(yōu)化、Leach協(xié)議優(yōu)化、WSN覆蓋優(yōu)化、組播優(yōu)化、RSSI定位優(yōu)化

6 信號(hào)處理方面

信號(hào)識(shí)別、信號(hào)加密、信號(hào)去噪、信號(hào)增強(qiáng)、雷達(dá)信號(hào)處理、信號(hào)水印嵌入提取、肌電信號(hào)、腦電信號(hào)、信號(hào)配時(shí)優(yōu)化

7 電力系統(tǒng)方面

微電網(wǎng)優(yōu)化、無(wú)功優(yōu)化、配電網(wǎng)重構(gòu)、儲(chǔ)能配置

8 元胞自動(dòng)機(jī)方面

交通流 人群疏散 病毒擴(kuò)散 晶體生長(zhǎng)

9 雷達(dá)方面

卡爾曼濾波跟蹤、航跡關(guān)聯(lián)、航跡融合