MSU的邢老师真是个牛人~

最新推荐文章于 2021-01-16 08:09:33 发布

原创最新推荐文章于 2021-01-16 08:09:33 发布 · 876 阅读

0 ·

CC 4.0 BY-SA版权

Matlab 专栏收录该内容

7 篇文章

订阅专栏

本文介绍了一种使用机器人传感器网络进行准确性感知水下扩散过程剖析的方法，并通过最大似然估计实现参数优化。文中提供了MATLAB代码示例，展示了如何模拟扩散过程并估计其参数。

最近在跟邢老师的Accuracy-Aware Aquatic Diffusion Process Proﬁling Using Robotic Sensor Networks, CRB这个东西真是让我很头大，论文里的参考文献是经典的《模式分类》，翻了半天都没找到CRB这个东西，汗一个。

这是部分做最大似然估计的代码，虽然图画得还凑合看，但参数好像有点不对。跟老杜讨论了一会儿，老杜提的一些想法我发现邢老师都已经全集成进这篇论文了，于是老杜一拍大腿：“TMD他好像把这个问题做死了~~”

当时就。。。想笑。。。

效果：

代码：

function [x,y,concentration]=diffusion_grediant_ascent(time,Amount,D)
%模拟扩散过程的Profiling
%参数为:时间 污染物总量 扩散系数
%梯度上升的办法
figure(1);
axis([-50,50,-50,50,0,1]);
%图1绘制随时间变化的浓度情况,注意到初始条件为一冲击函数
for t=0:0.1:time
for i=1:1:101
    for j=1:1:101
      x(i,j)=i-51;
      y(i,j)=j-51;
      distance=sqrt(x(i,j)*x(i,j)+y(i,j)*y(i,j));
      concentration(i,j)=(Amount/(4*pi*D*t))*exp(-(distance*distance)/4*D*time);
    end
end
subplot(2,2,1);
surf(x,y,concentration);
title('扩散过程');
axis([-50,50,-50,50,0,1000]);
drawnow;
end

hold on

pause;

%传感器数目
num_of_sensor=10;
%偏差值
bias=0; 
%sigma
sigma=10;
%在新窗口中选择传感器防止位置
[x,y]=select_sensors(num_of_sensor);
subplot(2,2,2);
axis([-50,50,-50,50]);
hold on
scatter(x,y,'p','MarkerEdgeColor','k','MarkerFaceColor','k');
for select_num=1:num_of_sensor
z(select_num)=concentration(x(select_num)+51,y(select_num)+51);%获取传感器信息
z(select_num)=z(select_num)+bias;%传感器度数存在偏差
z(select_num)=z(select_num)+normrnd(0,sigma);%传感器度数存在噪声
%text(x,y,' ');
end
title('投放传感器');

%fun=@(parameter)((z-parameter(1)/sigma*exp(-1*parameter(2)*((x-parameter(3))^2+(y-parameter(4))^2)))^2);

%最大似然估计优化目标
fun=@(parameter)(((z(1)-bias)/sigma-parameter(1)/sigma*exp(-1*parameter(2)*((x(1)-parameter(3))^2+(y(1)-parameter(4))^2)))^2+((z(2)-bias)/sigma-parameter(1)/sigma*exp(-1*parameter(2)*((x(2)-parameter(3))^2+(y(2)-parameter(4))^2)))^2+((z(3)-bias)/sigma-parameter(1)/sigma*exp(-1*parameter(2)*((x(3)-parameter(3))^2+(y(3)-parameter(4))^2)))^2+((z(4)-bias)/sigma-parameter(1)/sigma*exp(-1*parameter(2)*((x(4)-parameter(3))^2+(y(4)-parameter(4))^2)))^2+((z(5)-bias)/sigma-parameter(1)/sigma*exp(-1*parameter(2)*((x(5)-parameter(3))^2+(y(5)-parameter(4))^2)))^2+((z(6)-bias)/sigma-parameter(1)/sigma*exp(-1*parameter(2)*((x(6)-parameter(3))^2+(y(6)-parameter(4))^2)))^2+((z(7)-bias)/sigma-parameter(1)/sigma*exp(-1*parameter(2)*((x(7)-parameter(3))^2+(y(7)-parameter(4))^2)))^2+((z(8)-bias)/sigma-parameter(1)/sigma*exp(-1*parameter(2)*((x(8)-parameter(3))^2+(y(8)-parameter(4))^2)))^2+((z(9)-bias)/sigma-parameter(1)/sigma*exp(-1*parameter(2)*((x(9)-parameter(3))^2+(y(9)-parameter(4))^2)))^2+((z(10)-bias)/sigma-parameter(1)/sigma*exp(-1*parameter(2)*((x(10)-parameter(3))^2+(y(10)-parameter(4))^2)))^2);
%parameter(1)=aplpha parameter(2)=beta parameter(3)=x0 parameter(4)=y0

%优化时从准确值处开始优化
%alpha=Amount/(4*pi*D*time);
%beta=1/(4*D*time);
%x_real=0;
%y_real=0;

%优化时从[0,0,0,0]开始优化
[answer,fval,exitflag,output]=fminsearch(fun,[0,0,0,0]);
while (exitflag~=1)%直到优化终止条件为找到最小值才停止
[answer,fval,exitflag,output]=fminsearch(fun,[answer(1),answer(2),answer(3),answer(4)])
                 % alpha beta x0 y0                 
end

t_hat=(4*D*answer(2))^-1%对应的时间估计
A_hat=pi*answer(1)*(answer(2))^-1%对应的总量估计
x_hat=answer(3)%对应的当前源点中心位置估计
y_hat=answer(4)%对应的当前源点中心位置估计

%绘制估计的情况
subplot(2,2,3);
for i=1:1:101
    for j=1:1:101
      x(i,j)=i-51;
      y(i,j)=j-51;      
      distance=sqrt((x(i,j)-x_hat)^2+(y(i,j)-y_hat)^2);
      concentration_hat(i,j)=(A_hat/(4*pi*D*t_hat))*exp(-(distance*distance)/4*D*time);%依据扩散模型
    end
end
surf(x,y,concentration_hat);
axis([-50,50,-50,50,0,1000]);
title('估计值');

%绘制估计偏差
subplot(2,2,4);
surf(x,y,concentration_hat-concentration);
axis([-50,50,-50,50,-100,100]);
title('估计值偏差');

%计算CRB



%移动


end

function [x,y]=select_sensors(select_num)
%选择传感器
pic2=figure(2);
axis([-50,50,-50,50]);
title('请选择传感器投放位置:');
for i=1:select_num
    [x(i),y(i)]=ginput(1);
    x(i)=fix(x(i));%取整对应格点
    y(i)=fix(y(i));%取整对应格点
    scatter(x,y,'p','MarkerEdgeColor','k','MarkerFaceColor','k');
    axis([-50,50,-50,50]);
    title('请选择传感器投放位置:');
end
close(pic2)
end