实验二 RSSI测距定位技术性能程序设计

1、RSSI的测量值由对数路径损耗模型产生，为减小波动造成的误差，其值可由多次测量取平均值来得到。

2、对数路径损耗模型中的参考距离路径损耗和路径损耗因子可通过参考点相互之间的测量值估计。

3、完成理想情况下（参考距离路径损耗和路径损耗因子已知）与实际情况下的RMSE曲线对比图，横坐标为噪声方差，纵坐标为RMSE。

1、main主函数写入RSSI定位算法的仿真

  (1)第一段代码为每个BS的RSSI测量值

代码定义了四个固定节点[0,0] ;[500,0];BS3=[500,500]; BS4=[0,500]; 和一个移动节点MS=[100，500]，四个已知节点构成的坐标矩阵，仿真假设已知参考功率电平(pd0)为0db，路径损耗指数(n)为3。和用于进行RSSI测量的时间间隔(tt)。使用具有不同标准差的路径损失模型生成RSSI测量值，并对一定数量的锚节点的测量值进行平均。

图1 每个BS的RSSI测量值代码

（2）嵌套循环代码

用于模拟加性高斯噪声的不同标准差(std_var)的RSSI定位算法。外部循环迭代std_var的值，内部循环使用具有这些标准差的路径损耗模型生成RSSI测量值。

对于每组RSSI测量值，代码使用已知的参考功率电平和路径损耗指数(理想情况)或根据RSSI测量值(实际情况)估计这些参数来估计移动节点和固定节点之间的距离。

然后，代码使用到达时间(TOA)三角测量计算移动节点的估计位置，并计算移动节点的真实位置与其估计位置之间的均方根误差(RMSE)，在多次迭代中取平均值。

最后，代码为std_var的每个值存储理想情况和实际情况下的RMSE值。

图2 存储RMSE值

   （3）plot绘图代码

代码生成了RSSI测量(std_var)中均方根误差(RMSE)与加性高斯噪声标准偏差(std_var)的图，用于RSSI定位算法的理想和实际情况。该图显示了两条线，每一条代表一种情况，用标记表示std_var和RMSE的值。x轴表示噪声的标准差(以db为单位)，y轴表示RMSE。xlabel和ylabel函数标记坐标轴，legend函数创建一个图例，表示理想的和实际的情况。

图3 绘制结果图代码

2、parameter_est.m

（1）该函数基于来自固定节点(BS)的RSSI测量。

  以固定节点的坐标(A)和RSSI测量中加性高斯噪声的标准差(sigma)作为输入。首先计算每对固定节点之间的距离，并使用已知值pd0和n的路径损耗模型生成RSSI测量值。

  然后构建RSSI测量值与固定节点之间距离的线性方程组(Ax = B)

其中a为距离矩阵，x = [pd0, n]为未知参数向量，B为RSSI测量值向量。

使用最小二乘回归求解pd0和n的参数估计值。

将这些估计值分别返回为pd0_est和n_est。

图4 parameter_est函数

3、TOALLOP函数

  此函数基于一组参考节点(BBS)的到达时间(TOA)测量来估计移动节点的二维位置。

  以参考节点的坐标(A)、TOA测量值(p)和用于三角测量的参考节点索引(j)作为输入。

  该函数首先使用TOA测量值计算移动节点与每个参考节点之间的距离。然后，构建了一个将距离平方与参考节点坐标相关的线性方程组。该函数使用最小二乘回归来求解移动节点的未知位置，给定参考节点的平方距离和坐标。

  最后，该函数返回移动节点的估计二维位置，作为x和y坐标(theta)的矢量。

图5 theta函数

4、运行结果

图6 运行结果

代码

main

%RSSI定位算法的仿真
clear all;
clc;
BS1=[0,0];
BS2=[500,0];
BS3=[500,500];
BS4=[0,500];
MS=[100,500];
std_var=[0,2,4,6];%标准偏差
A=[BS1;BS2;BS3;BS4];%A表示每一个已知节点构成坐标矩阵
number=300;
pd0=0;%pd0接收功率=0db
n=3;%n表示路径损耗因子
tt=5;%每5个锚节点测量RSSI值求平均得到RSSI值
%tt=10;
%tt=30;
% 每个BS的RSSI测量值
for j=1:length(std_var)
    error1=0;
    error2=0;
    std_var1=std_var(j);
    for i=1:number%循环次数number800
        r1=A-ones(4,1)*MS;
        r2=(sum(r1.^2,2)).^(1/2);
        for k=1:tt%内部再循环，因为多次测量
            rssi(:,k)=pd0-10*n*log10(r2)-10^(std_var1/10)*randn(4,1);%将db转换为非db,产生对数正态分布得到RSSI
        end
        RSSI1=mean(rssi,2);%求平均得到未知节点到已知节点的RSSI值，估计未知节点的坐标
        % 理想情况,prd0、n已知
        r1=10.^((RSSI1-pd0)/(-10*n));%得到r1距离
        % 实际情况
        [p_est,n_est]=parameter_est(A,std_var1);
        r2=10.^((RSSI1-p_est)/(-10*n_est));%r2维移动台未知节点到已知节点距离
        theta1=TOALLOP(A,r1,1);%输入A，已知r1,参考节点为1得到theta1估计
        theta2=TOALLOP(A,r2,1);%得到theta2估计
        error1=error1+norm(MS-theta1)^2;%均方误差
        error2=error2+norm(MS-theta2)^2;
    end
    RMSE1(j)=(error1/number)^(1/2);%理想情况RSSI标准偏差
    RMSE2(j)=(error2/number)^(1/2);%实际情况RSSI标准偏差
end
% plot
plot(std_var,RMSE1,'-O',std_var,RMSE2,'-s')
xlabel('RSSI标准偏差 (db)');
ylabel('均方根误差RMSE');
legend('理想情况','实际情况');

 parameter_est.m

function [pd0_est,n_est]=parameter_est(A,sigma)%参数估计函数，A是已知节点构成的坐标，sigma是标准偏差
% A是很多BS的坐标
% sigma是RSSI测量的标准差
[m,~]=size(A);%有多少个已知节点参加定位
pd0=0;
n=3;
d=zeros(m,m);%初始化固定节点坐标为0
tt=5;
% 每个BS的RSSI测量的个数
sigma1=10^(sigma/10);%db转化为非db
h1=[];
G1=[];
for i=1:m
    for j=1:m
        if i~=j %两个节点距离不想等
            d(i,j)=norm(A(i,:)-A(j,:));%计算几种不同类型的矩阵范数
            for k=1:tt
                prd(k)=pd0-10*n*log10(d(i,j))-sigma1*randn;%距离为k的每个测量值的prd
            end
            RSSI=mean(prd);%平均值
            d_distance=-10*log10(d(i,j));
            h1=[h1;RSSI];%对应Ax=B的B
            G1=[G1;d_distance];%G1对应Ax=B的A
        end
    end
end
h=h1;
[m1,~]=size(h);
G=[ones(m1,1),G1];%将接收功率的式子转化为矩阵形式，构造G为n*2维矩阵
x=inv(G'*G)*G'*h;%最小二乘法得到位置估计
pd0_est=x(1,1);%pd0估计
n_est=x(2,1);%路径损耗因子估计
end

 TOALLOP.m

function theta=TOALLOP(A,p,j)
%A是BBS的坐标
%P是范围测量
%J是参考BS的索引
[m,~]=size(A);  %size得到A的行列数赋值给[m,~]，~表示占位，就是只要行m的值
k=sum(A.^2,2);%矩阵A每个元素分别平方，得到新矩阵，在行求和，最为矩阵K
k1=k([1:j-1,j+1:m],:); %取出J行
A1=A([1:j-1,j+1:m],:); %取出J行
A2=A1-ones(m-1,1)*A(j,:); %得到D，就是j行与其余行对应值相减
p1=p([1:j-1,j+1:m],:); %取出J行
p2=p(j).^2*ones(m-1,1)-p1.^2-(k(j)*ones(m-1,1)-k1); %得到b,(Rn*Rn-R1*R1-Kn+K1)其中Kn为对应第n个x^2+y^2
theta=1/2*inv(A2'*A2)*A2'*p2; %利用最小二乘解，得位置估计
theta=theta';%转换为（x,y）形式