2021華為杯E題思路+demo代碼

2021 年中國研究生數學建模競賽 E 題參考思路
交流群：912166339，非伸手黨群
信號干擾下的超寬帶（UWB）精確定位問題
一、背景
UWB（Ultra-Wideband）技術也被稱之為“超寬帶”，又稱之為脈沖無線電技術，這是一
種無需任何載波，通過發送納秒級脈沖而完成資料傳輸的短距離范圍內無線通信技術，并且信
號傳輸程序中的功耗僅僅有幾十μW，UWB 因其獨有的特點，使其在軍事、物聯網等各個領域
都有著廣闊的應用，其中，基于 UWB 的定位技術具備實時的室內外精確跟蹤能力，定位精度
高，可達到厘米級甚至毫米級定位，UWB 在室內精確的定位將會對衛星導航起到一個極好的
補充作用，可在軍事及民用領域有廣泛應用，比如：電力、醫療、化工行業、隧道施工、危險
區域管控等，UWB 更多應用場景請參見[4—6]，
UWB 的定位技術有多種方法，本文僅考慮基于飛行時間（Time of Flight, TOF）的測距原
理，它是 UWB 定位法中最常見的定位方法之一，TOF 測距技術屬于雙向測距技術，其通過計
算信號在兩個模塊的飛行時間，再乘以光速求出兩個模塊之間的距離，這個距離肯定有不同程
度的誤差，但其精度已經比較高，
在室內定位的應用中，UWB技術可以實作厘米級的定位精度（一般指2維平面定位）,并
具有良好的抗多徑干擾和衰弱的性能以及具有較強的穿透能力，但由于室內環境復雜多變
UWB 通信信號極易受到遮擋，雖然UWB技術具有穿透能力，但仍然會產生誤差，在較強干
擾時，資料會發生例外波動（通常是時間延時），基本無法完成室內定位，甚至會造成嚴重事
故，因此，信號干擾下的超寬帶（UWB）精確定位問題成為亟待解決的問題，
二、問題描述
為解決信號干擾下的超寬帶（UWB）精確定位問題，我們通過實際場景實測，采集到一
定數量的資料，即利用 UWB 的定位技術（TOF），采集到錨點（
anchor）與靶點（Tag）之間
的距離，希望通過數學建模（或演算法）方法 ，無論信號是否干擾，都可以給出目標物（靶點）
的精確定位（
3 維坐標），
三、實驗場景和資料采集
如圖所示，在 5000mm5000mm3000mm 的測驗環境中，分別在 4 個角落 A0，A1，A2，
A3 放置 UWB 錨點（
anchor），錨點向所有方向發送信號，Tag 是 UWB 標簽（靶點），即需
要定位的目標（只在測驗環境范圍內），Tag 接收到 4 個 UWB 錨點（
anchor）的信號（無論
信號是否干擾，Tag 一般都可以接收到信號），利用 TOF 技術，分別解算出對應的 4 個距離數
據，
實驗在實驗場景 1 中采集了 Tag 在 324 個不同位置，在信號無干擾和信號干擾下的 UWB
資料，即每個位置各測驗（采集）2 次，一次信號無干擾，另一次信號有干擾（錨點與靶點間
有遮擋），注意：每次采集資料時，由于 Tag 在同一位置會停留一會兒時間，而錨點與 Tag 之
間每 0.2—0.3 秒之間就會發送、接收信號一次，所以在同一位置點，UWB 會采集到多組資料
（多組資料都代表同一位置的資訊），組數的多少視 Tag 在同一位置的時間而定，停留的時間
越長，組數就越多，資料見檔案夾“附件 1：UWB 資料集”，
圖 1
實測環境示意圖
實驗場景 1：
靶點（Tag）范圍：5000mm5000mm3000mm
錨點（
anchor）位置（單位：mm）：
A0（
0，0，1300）、
A1（
5000，0，1700）、
A2（
0，5000，1700）、A3（
5000，5000，1300）
四、資料檔案說明
（
1）UWB 資料集
“附件 1：UWB 資料集”有 2 個檔案夾和 1 個檔案，1 個檔案（Tag 坐標資訊.txt）存放
324 個不同位置的編號及 3 維坐標資訊，2 個檔案夾中 1 個存放信號無干擾下（正常）采集的
資料（各檔案名為 x.正常.txt，x 表示對應的位置編號），另 1 個存放信號有干擾下（例外）采
集的資料（各檔案名為 x.例外.txt，x 表示對應的位置編號），
（
2）資料檔案
Tag 在每個位置都采集了 2 個資料檔案（1 個正常，另 1 個例外），共有 648 個資料檔案，
無論正常、例外資料，資料格式都一樣，每個資料檔案開頭第 1 行為采集開始行，無實際意義，
接下來，每 4 行為一組，表示 UWB 采集的一組完整資料（一組資料表示一個樣品），如：
T:144235622:RR:0:0:950:950:118:1910
T:144235622:RR:0:1:2630:2630:118:1910
T:144235622:RR:0:2:5120:5120:118:1910
T:144235622:RR:0:3:5770:5770:118:1910
這 4 行資料的含義分別是:
Tag 標識：時間戳：Range Report 的縮寫：Tag ID：錨點 ID：該錨點的測距值(mm)：測
距值的校驗值：資料序列號：資料編號（每個資料之間用“：”分隔），實際上就是提供了 4
個錨點到靶點（Tag）的距離，即
A0 到靶點距離為：950mm
A1 到靶點距離為：2630mm
A2 到靶點距離為：5120mm
A3 到靶點距離為：5770mm
每個資料檔案都有多組資料，表示在同一位置連續時間內 UWB 自動采集到的多組資料，
五、完成任務研賽更多資料請加群：678535491，或關注公眾號：老哥帶你學術數模
試根據上述資料，完成如下任務：
任務 1：資料預處理（清洗）
無論是信號無干擾下采集資料，或信號干擾下采集資料，Tag 在同一坐標點上都采集多組
資料（見附件 1 中 648 個資料檔案），請用某種方法把每個資料檔案相應數值抓取出來，并轉
換成二維表（矩陣）形式（
txt、Excel 或其他資料格式），每一行代表一組資料（即一個樣品），
然后對這些資料檔案進行預處理（清洗），洗掉掉一些“無用”（例外、缺失、相同或相似）
的資料（樣品），經處理后，“正常資料”所有資料檔案和“例外資料”所有資料檔案最后各
保留多少組（多少個樣品）資料，并重點列出以下 4 個資料檔案，經處理后保留的資料（矩陣
形式）；
“正常資料”檔案夾中：
24.正常.txt、 109.正常.txt
“例外資料”檔案夾中：
1.例外.txt、
100.例外.txt
對于任務一而言，主要是進行資料的預處理作業；首先需要用 matlab 撰寫一套讀取檔案的代
碼，參考如下：
List =dir(‘C:\Users\DELL\Desktop\data*.txt’); %設定路徑
for i=1:1:length(List) %挨個讀入 txt 檔案資料
A{i}=importdata([‘C:\Users\DELL\Desktop\data’,List(i).name]);
End
然后對資料做相關的處理作業即可；可選擇不同的方法將無用資料提取出來，其中例外資料可
以用簡單統計分析、3σ原則或者箱線圖來計算；對于缺失資料和相同及相似資料直接用程式
判斷即可，找到上述例外資料后直接剔出，
3σ原則是建立在正態分布的等精度重復測量基礎上而造成奇異資料的干擾或噪聲難以滿足正
態分布，研賽更多資料請加群：678535491，或關注公眾號：老哥帶你學術數模
正態分布：又叫做高斯分布，特征為中間高兩邊低左右對稱，正態分布特性：
集中性：曲線的最高峰位于正中央，并且位置為均數所在的位置，
對稱性：以均數所在的位置為中心呈左右對稱，并且曲線兩段無限趨近于橫軸，
均勻變動性：正態分布曲線以均數所在的位置為中心均勻向左右兩側下降，
正態分布函式公式如下：
如果一組測量資料中的某個測量值的殘余誤差的絕對值 ν
i>3σ，則該測量值為壞值，應剔除，
通常把等于 ±3σ的誤差作為極限誤差,對于正態分布的隨機誤差，落在 ±3σ以外的概率只
有 0.27%，它在有限次測量中發生的可能性很小，故存在 3σ準則，3σ準則是最常用也是最
簡單的粗大誤差判別準則，它一般應用于測量次數充分多( n ≥30)或當 n>10 做粗略判別時的
情況，
σ代表標準差,μ代表均值
樣本資料服從正態分布的情況下：
數值分布在(μ-σ,μ+σ)中的概率為 0.6826
數值分布在(μ-2σ,μ+2σ)中的概率為 0.9544
數值分布在(μ-3σ,μ+3σ)中的概率為 0.9974
可以認為，Y 的取值幾乎全部集中在（μ-3σ,μ+3σ)區間內，超出這個范圍的可能性僅占不
到 0.3%，
任務 2: 定位模型
利用任務 1 處理后的資料，分別對“正常資料”和“例外資料”，設計合適的數學模型（或
演算法），估計（或預測）出 Tag 的精確位置，并說明你所建立的定位模型（或演算法）的有效性；
同時請利用你的定位模型（或演算法）分別對附件 2 中提供的前 5 組（信號無干擾）資料和后 5
組（信號有干擾）資料進行精確定位（
3 維坐標）；
注意：（1）定位模型必須體現實驗場景資訊；
（
2）請同時給出定位模型的 3 維（
x，y，z）精度、2 維（
x，y）精度以及 1 維的各自精度，
根據常識可知，若已知四個點的坐標，并且知道這四個點到另一個點的距離，此時能在空間上
確定該點的坐標；但題目中給了四個已知點，但存在大量的時間戳，因此存在多余解的問題，
對于多余解的問題則可以采用最小二乘的思想進行計算即可得到最終的點位坐標；而對于例外
資料，由于錨點和靶點間存在遮擋，因此此時計算的位置資料將存在偏差；此時可以先按照前
文正常資料計算的思路來時，將會發現計算得到的四組資料存在明顯的偏差；因此可以首先將
可能認為是例外的點剔除（如統計所有樣本資料將方差較大的值洗掉），然后利用其它正常點
來計算坐標即可，
另外在測量領域精度一定要用 3σ理念來表示，比較專業，
任務 3：不同場景應用
我們的訓練資料僅采集于同一實驗場景（實驗場景 1），但定位模型應該能夠在不同實際
場景上使用，我們希望你所建立的定位模型能夠應用于不同場景，附件 3 中 10 組資料采集于
下面實驗場景 2（前 5 組資料信號無干擾，后 5 組資料信號有干擾），請分別用上述建立的定
位模型，對這 10 組資料進行精確定位（
3 維坐標）；
實驗場景 2：
靶點（Tag）范圍：5000mm3000mm3000mm
錨點（
anchor）位置（單位：mm）：
A0（
0，0，1200）、
A1（
5000，0，1600）、
A2（
0，3000，1600），A3（
5000，3000，1200）
直接用任務 2 構建的模型進行計算即可
任務 4: 分類模型
上述定位模型是在已知信號有、無干擾的條件下建立的，但 UWB 在采集資料時并不知道
信號有無干擾，所以判斷信號有無干擾是 UWB 精確定位問題的重點和難點，利用任務 1 處理
后的資料，建立數學模型(或演算法)，以便區分哪些資料是在信號無干擾下采集的資料，哪些數
據是在信號干擾下采集的資料？并說明你所建立的分類模型(或演算法)的有效性；同時請用你所
建立的分類模型（或演算法）判斷附件 4 中提供的 10 組資料（這 10 組資料同樣采集于實驗場景
1）是來自信號無干擾或信號干擾下采集的？
這里是一道典型的分類問題，由于實作并不知道那些是有信號干擾的，那些是沒有信號干擾的，
因此可以采用非監督聚類的思想進行分類作業，非監督聚類模型的主要有 K-means 聚類、高斯
混合聚類和層次聚類等，可以分別試試，類別寫兩類，這樣最終就能將哪些資料是例外資料而
那些是正常資料進行劃分；對于分類演算法有效性判斷可以利用附件 1：UWB 資料集進行驗證，
任何分類演算法都是又誤差的，這時候需要把誤差進行分析描述，如果誤差較小則將其應用于附
件 4 中提供的 10 組資料即可，
K-means 演算法代碼
clear all;
close all;
clc;
% 第一組資料
mu1=[0 0 ]; %均值(是需要生成的資料的均值)
S1=[.1 0 ;0 .1]; %協方差(需要生成的資料的自相關矩陣（相關系數矩陣）)
data1=mvnrnd(mu1,S1,100);
%產生高斯分布資料
%第二組資料
mu2=[1.25 1.25 ];
S2=[.1 0 ;0 .1];
data2=mvnrnd(mu2,S2,100);
% 第三組資料
mu3=[-1.25 1.25 ];
S3=[.1 0 ;0 .1];
data3=mvnrnd(mu3,S3,100);
% 顯示資料
plot(data1(:,1),data1(:,2),‘b+’);
hold on;%不覆寫原圖，要關閉則使用 hold off；
plot(data2(:,1),data2(:,2),‘r+’);
plot(data3(:,1),data3(:,2),‘g+’);
grid on;%顯示表格
% 三類資料合成一個不帶標號的資料類
data=[data1;data2;data3];
N=3;%設定聚類數目
[m,n]=size(data);%表示矩陣 data 大小，m 行 n 列
pattern=zeros(m,n+1);%生成 0 矩陣
center=zeros(N,n);%初始化聚類中心
pattern(:,1:n)=data(:😅;
for x=1:N
center(x,:)=data( randi(300,1)😅;%第一次隨機產生聚類中心
end
while 1 %回圈迭代每次的聚類簇；
distence=zeros(1,N);%最小距離矩陣
num=zeros(1,N);%聚類簇數矩陣
new_center=zeros(N,n);%聚類中心矩陣
for x=1:m
for y=1:N
distence(y)=norm(data(x,:)-center(y,:));%計算到每個類的距離
end
[~, temp]=min(distence);%求最小的距離
pattern(x,n+1)=temp;%劃分所有物件點到最近的聚類中心；標記為 1,2,3；
end
k=0;
for y=1:N
for x=1:m
if pattern(x,n+1)y
new_center(y,:)=new_center(y,:)+pattern(x,1:n);
num(y)=num(y)+1;
end
end
new_center(y,:)=new_center(y,:)/num(y);%求均值，即新的聚類中心；
if norm(new_center(y,:)-center(y,:))<0.1%檢查集群中心是否已收斂，如果是則終止，
k=k+1;
End
end
if kN
break;
else
center=new_center;
end
end
[m, n]=size(pattern);
%最后顯示聚類后的資料
figure;
hold on;
for i=1:m
if pattern(i,n)1
plot(pattern(i,1),pattern(i,2),‘r*’);
plot(center(1,1),center(1,2),‘ko’);%用小圓圈標記中心點；
elseif pattern(i,n)2
plot(pattern(i,1),pattern(i,2),‘g*’);
plot(center(2,1),center(2,2),‘ko’);
elseif pattern(i,n)3
plot(pattern(i,1),pattern(i,2),‘b*’);
plot(center(3,1),center(3,2),‘ko’);
elseif pattern(i,n)4
plot(pattern(i,1),pattern(i,2),‘y*’);
plot(center(4,1),center(4,2),‘ko’);
else
plot(pattern(i,1),pattern(i,2),‘m*’);
plot(center(4,1),center(4,2),‘ko’);
end
end
grid on;
高斯混合聚類代碼
%% 初始化作業空間
clc
clear all
close all
%% 載入資料
load fisheriris
%% 二維資料
% 葉子長度和葉子寬度散點圖(真實標記)
figure;
speciesNum = grp2idx(species);
gscatter(meas(:,3),meas(:,4),speciesNum,[‘r’,‘g’,‘b’])
xlabel(‘葉子長度’)
ylabel(‘葉子寬度’)
title(‘真實標記’)
set(gca,‘FontSize’,12)
set(gca,‘FontWeight’,‘bold’)
% 葉子長度和葉子寬度散點圖（無標記）
figure;
scatter(meas(:,3),meas(:,4),150,’.’)
xlabel(‘葉子長度’)
ylabel(‘葉子寬度’)
title(‘無標記’)
set(gca,‘FontSize’,12)
set(gca,‘FontWeight’,‘bold’)
%% 高斯混合模型聚類
data = [meas(:,3), meas(:,4)];
% 手動初始條件
Mu=[0.25 1.5; 4.0 1.25; 5.5 2.0 ];%平均值設定
Sigma(:,:,1) = [1 1;1 2];%方差設定
Sigma(:,:,2) = [1 1;1 2];
Sigma(:,:,3) = [1 1;1 2];
Pcom=[1/3 1/3 1/3];%混合比例
S = struct(‘mu’,Mu,‘Sigma’,Sigma,‘ComponentPropotion’,Pcom);
GMModel = fitgmdist(data,3,‘Start’,S);
% 分類
T1 = cluster(GMModel,data);
% 標號調整
cen=[mean(data(T11,:));…
mean(data(T12,:));…
mean(data(T13,:))];
dist=sum(cen.^2,2);
[dump,sortind]=sort(dist,‘ascend’);
newT1=zeros(size(T1));
for i =1:3
newT1(T1i)=find(sortindi);
end
% 葉子長度和葉子寬度散點圖(真實標記:實心圓+kmeans 分類:圈)
figure;
gscatter(meas(:,3),meas(:,4),speciesNum,[‘r’,‘g’,‘b’])
hold on
gscatter(data(:,1),data(:,2),newT1,[‘r’,‘g’,‘b’],‘o’,10)
scatter(cen(:,1),cen(:,2),300,‘m*’)
hold off
xlabel(‘葉子長度’)
ylabel(‘葉子寬度’)
title(‘真實標記:實心圓+kmeans 分類:圈’)
set(gca,‘FontSize’,12)
set(gca,‘FontWeight’,‘bold’)
%% 混淆矩陣 ConfusionMatrix
T2ConfMat=confusionmat(speciesNum,newT1)
error23=(speciesNum2)&(newT13);
errDat23=data(error23,:)
error32=(speciesNum3)&(newT1==2);
errDat32=data(error32,:)
%% 高斯模型等高線圖
% 散點圖
figure;
gscatter(meas(:,3),meas(:,4),speciesNum,[‘r’,‘g’,‘b’])
hold on
scatter(cen(:,1),cen(:,2),300,‘m*’)
hold off
xlabel(‘葉子長度’)
ylabel(‘葉子寬度’)
title(‘高斯模型等高線圖’)
set(gca,‘FontSize’,12)
set(gca,‘FontWeight’,‘bold’)
% 疊加等高線
haxis=gca;
xlim = get(haxis,‘XLim’);
ylim = get(haxis,‘YLim’);
dinter=(max([xlim, ylim]) - min([xlim, ylim]))/100;
[Grid1, Grid2] = meshgrid(xlim(1):dinter:xlim(2), ylim(1):dinter:ylim(2));
hold on
GMMpdf=reshape(pdf(GMModel, [Grid1(😃 Grid2(😃]), size(Grid1,1), size(Grid2,2));
contour(Grid1, Grid2, GMMpdf, 30);
hold off
% 混合高斯模型曲面圖
figure;
surf(GMMpdf)
xlabel(‘葉子長度’)
ylabel(‘葉子寬度’)
title(‘高斯混合模型曲面圖’)
view(-3,65)
層次聚類相關代碼
%% 讀取資料
clc;clear all,close all
pre_data=xlsread(‘StockFinance.xlsx’,‘Sheet1’,‘E2:N4735’);
%% 資料歸一化
[rn,cn]=size(pre_data);
step1_data=zeros(rn,cn);
for k=1:cn
%基于均值方差的離群點資料歸一化
xm=mean(pre_data(:,k));
xs=std(pre_data(:,k));
for j=1:rn
if (pre_data(j,k))>xm+2xs
step1_data(j,k)=1;
elseif (pre_data(j,k))<xm-2xs
step1_data(j,k)=0;
else
step1_data(j,k)=(pre_data(j,k)-(xm-2xs))/(4xs);
end
end
end
xlswrite(‘norm_data.xlsx’,step1_data);
%% 層次聚類
numClust=3;
dist_h=‘spearman’;
link=‘weighted’;
hidx=clusterdata(step1_data,‘maxclust’,numClust,‘distance’,dist_h,‘linkage’,link);
任務 5: 運動軌跡定位
運動軌跡定位是 UWB 重要應用之一，利用靜態點的定位模型，加上靶點自身運動規律，
希望給出動態靶點的運動軌跡，附件 5 是對動態靶點采集的資料（一段時間內連續采集的多組
資料），請注意，在采集這些資料時，會隨機出現信號干擾，請對這個運動軌跡進行精確定位，
最侄訓出這條運動軌跡圖（資料采集來自實驗場景 1），
這個就是結合前面各小問的解決方案進行的集成，首先題目明確告知了會隨機出現干擾，，也
就意味著在整個時間線內即存在正常資料也存在干擾資料，因此因該先用附件 4 中的聚類模型
將干擾資料進行篩選，然后利用任務 2 構建的基于干擾資料的定位模型進行點位坐標的確定，
再結合正常資料的坐標即可得到最終的運動軌跡資訊，
參考文獻
[1]Kok M, Hol J D, Schon T B. Indoor Positioning Using Ultrawideband and Inertial
[2]Measurements[J]. IEEE Transactions on Vehicular Technology, 2015, 64(4): 1293—1303.
[2] Arias-De-Reyna E. A Cooperative Localization Algorithm for UWB Indoor Sensor Networks[J].
Wireless Personal Communications, 2013, 72(1):85—99.
[3]繆希仁,范建威等，基站例外情況下基于改進極限學習機的超寬帶室內定位方法，傳感技術
學報，2020，33（
10）：1—10.
[4]http://m.elecfans.com/article/1643155.html
[5]http://m.elecfans.com/article/1116045.html
[6]http://www.qianjia.com/zhike/html/2020-02/11_19601.html