PNN实现重力匹配——MATLAB复现论文

论文《一种模式识别神经网络重力匹配算法》2007中提出模式识别的方法进行重力匹配。通过概率神经网络处理重力测量值序列，结合惯性导航系统的轨迹与重力背景图，在重力背景图中找出一组与重力测量值最为相似的重力序列，以此作为重力匹配的结果。

论文匹配方法原理如下图所示：

通过在重力背景图搜索区域中平移惯性导航系统给出的轨迹，得到对应的重力值序列，作为模式k的特征向量。详情见论文。

在测试过程中，假设重力测量值误差非常小（几乎为0）时，才可能得到更精确的重力匹配结果。并且论文中的节点偏置系数（如下图所示）对模式识别的结果影响很大。若该系数设置不当，可能无法得到正确的重力匹配结果，在使用过程中需要斟酌。

MATLAB代码如下所示：

% 《一种模式识别神经网络重力匹配算法》
% 利用概率神经网络进行重力匹配

clear all;
clc;
%% 训练数据
load x2;
load y2;
load z2;

M=3;    %搜索区域左右间隔
N=3;     %上下间隔
load Storedata_caiyang;   %真实轨迹
load INS_caiyang;  %INS轨迹
load g_m;   %重力测量值

% 相对位置关系
L=length(INS_caiyang);      % 采样点个数
del_lamda=zeros(L-1,1);     % 相邻两个采样点之间的经度相对距离
del_fai=zeros(L-1,1);           % 纬度相对距离
for i=1:L-1
    del_lamda(i)=abs(INS_caiyang(L,1)-INS_caiyang(L-i,1));   %经度
    del_fai(i)=abs(INS_caiyang(L,2)-INS_caiyang(L-i,2));    %纬度
end

center_point=zeros(1,2);     
center_point=INS_caiyang(L,1:2);   % 以INS轨迹最后一个点为搜索区域的中心
center_num=find_num(center_point);   %在重力背景图中的坐标

% 定义第k个模式类
S = (2*M+1)*(2*N+1);     %共S个模式类
P = zeros(L,S);
Bk_num=zeros(1,2);
g_k = zeros(L,1);             %重力值序列
for k=1:S
    R1=floor(k/(2*N+1));    
    i_num(k) = center_num(1)-M+R1;    % 搜索区域中，从上到下从左到右第k个网格点的横坐标
    j_num(k) = center_num(2)+N-k+R1*(2*N+1)+1;    % 纵坐标
    lamda(k) = x2(1, i_num(k));   % 经纬度坐标
    fai(k) = y2(j_num(k), 1);
    g_k(1) = z2(j_num(k),i_num(k));    % 对应重力背景图值
    for j=1:L-1 
        Bk_lamda(j) = lamda(k) - del_lamda(j);   % 根据相对位置计算k模式类轨迹上各点的经度
        Bk_fai(j) = fai(k) - del_fai(j);    % 纬度
        point(1) = Bk_lamda(j);
        point(2) = Bk_fai(j);
        Bk_num= find_num(point);   % 对应坐标
        % 模式k对应的特征向量（重力值）
        g_k(j+1) = z2(Bk_num(2),Bk_num(1));
    end
    % 输入矩阵P，包含重力值序列
    P(:,k) = g_k(:);
end
% 模式指示矩阵T
for i=1:S
    T(1,i) = i;
end
for i=1:L
    g_m_new(i, 1) = g_m(L-(i-1), 1);   % 重力测量时间顺序：从后向前
end
%% 神经网络参数
% 重力粗糙程度
start_point=zeros(1,2);     % 搜索区域中心点
start_point=INS_caiyang(1,1:2);
start_num=find_num(start_point);   %在重力背景图中的坐标
X0 = center_num(1) - start_num(1);  %轨迹跨越网格数
Y0 = center_num(2) - start_num(2);
X = 2*M+X0;
Y = 2*N+Y0;
left_end_point(1) = start_num(1) - M;    %搜索区域左上角顶点坐标
left_end_point(2) = start_num(2) - Y0 - N;
% 经度方向
Q = 0;    %粗糙度
for i=1:(X-1)
    for j=1:Y
        H1 = z2(left_end_point(2) + j -1,left_end_point(1) + i - 1);
        H2 = z2(left_end_point(2) + j -1,left_end_point(1) + (i+1) - 1);
        Q = Q +  (H1 - H2)^2;
    end
end
Q_lamda = sqrt(1/((X-1)*Y)*Q);
% 纬度方向
Q = 0;    %粗糙度
for i=1:X
    for j=1:(Y-1)
        H1 = z2(left_end_point(2) + j -1,left_end_point(1) + i - 1);
        H2 = z2( left_end_point(2) + (j+1) -1,left_end_point(1) + i - 1);
        Q = Q +  (H1 - H2)^2;
    end
end
Q_fai= sqrt(1/(X*(Y-1))*Q);
sigma_z = (Q_lamda + Q_fai)/2;   %平均粗糙度

%% 建立PNN
p = P;
T2 = ind2vec(T);     % 生成单位矩阵
t = T2;
spread = 0.5;   %暂时用不到
net = newpnn_test(P, T2, spread);
%% 利用样本训练和测试PNN
x = g_m_new;
y = net(x);      %利用重力测量值测试神经网络
Yc = vec2ind(y);

plot(INS_caiyang(:,1),INS_caiyang(:,2),'k');
hold on;
plot(Storedata_caiyang(:,1),Storedata_caiyang(:,2),'r');
hold on;

k = Yc;
R1=floor(k/(2*N+1));
i_num(k) = center_num(1)-M+R1;    %搜索区域中，从上到下从左到右第k个网格点的经度坐标
j_num(k) = center_num(2)+N-k+R1*(2*N+1)+1;    %纬度
lamda(k) = x2(1, i_num(k));   %经纬度坐标
fai(k) = y2(j_num(k), 1); 
Bk(1,1) = lamda(k);
Bk(1,2) = fai(k);
for j=1:L-1
    Bk_lamda(j) = lamda(k) - del_lamda(j);
    Bk_fai(j) = fai(k) - del_fai(j);
    point(1) = Bk_lamda(j);
    point(2) = Bk_fai(j);
    Bk_num= find_num(point);
    Bk(j,1:2) = point(1:2); 
end
plot(Bk(:,1),Bk(:,2),'b');
hold on;

其中newpnn_test函数是根据MATLAB自带的newpnn修改得到的，主要修改了节点偏置部分的定义。MATLAB代码如下所示：

function out1 = newpnn_test(varargin)
%NEWPNN Design a probabilistic neural network.
%
%  Probabilistic neural networks are a kind of radial
%  basis network suitable for classification problems.
%
%  <a href="matlab:doc newpnn">newpnn</a>(P,T,SPREAD) takes an RxQ input matrix P and an SxQ target matrix
%  T, a radial basis function SPREAD and returns a new probabilistic
%  neural network.
%
%  If SPREAD is near zero the network will act as a nearest
%  neighbor classifier.  As SPREAD becomes larger the designed
%  network will take into account several nearby design vectors.
%
%  Here a classification problem is defined with a set of
%  inputs P and class indices Tc.  A PNN is designed to fit this data.
%
%    P = [1 2 3 4 5 6 7];
%    Tc = [1 2 3 2 2 3 1];
%    T = <a href="matlab:doc ind2vec">ind2vec</a>(Tc)
%    net = <a href="matlab:doc newpnn">newpnn</a>(P,T);
%    Y = net(P)
%    Yc = <a href="matlab:doc vec2ind">vec2ind</a>(Y)
%
%  See also SIM, IND2VEC, VEC2IND, NEWRB, NEWRBE, NEWGRNN.

% Mark Beale, 11-31-97
% Copyright 1992-2014 The MathWorks, Inc.

%% =======================================================
%  BOILERPLATE_START
%  This code is the same for all Network Functions.

  persistent INFO;
  if isempty(INFO), INFO = get_info; end
  if (nargin > 0) && ischar(varargin{1}) ...
      && ~strcmpi(varargin{1},'hardlim') && ~strcmpi(varargin{1},'hardlims')
    code = varargin{1};
    switch code
      case 'info',
        out1 = INFO;
      case 'check_param'
        err = check_param(varargin{2});
        if ~isempty(err), nnerr.throw('Args',err); end
        out1 = err;
      case 'create'
        if nargin < 2, error(message('nnet:Args:NotEnough')); end
        param = varargin{2};
        err = nntest.param(INFO.parameters,param);
        if ~isempty(err), nnerr.throw('Args',err); end
        out1 = create_network(param);
        out1.name = INFO.name;
      otherwise,
        % Quick info field access
        try
          out1 = eval(['INFO.' code]);
        catch %#ok<CTCH>
          nnerr.throw(['Unrecognized argument: ''' code ''''])
        end
    end
  else
    [args,param] = nnparam.extract_param(varargin,INFO.defaultParam);
    [param,err] = INFO.overrideStructure(param,args);
    if ~isempty(err), nnerr.throw('Args',err,'Parameters'); end
    net = create_network(param);
    net.name = INFO.name;
    out1 = init(net);
  end
end

function v = fcnversion
  v = 7;
end

%  BOILERPLATE_END
%% =======================================================

function info = get_info
  info = nnfcnNetwork(mfilename,'Probabilistic Neural Network',fcnversion, ...
    [ ...
    nnetParamInfo('inputs','Input Data','nntype.data',{},...
    'Input data.'), ...
    nnetParamInfo('targets','Target Data','nntype.data',{},...
    'Target output data.'), ...
    nnetParamInfo('spread','Radial basis spread','nntype.strict_pos_scalar',0.1,...
    'Distance from radial basis center to 0.5 output.'), ...
    ]);
end

function err = check_param(param)
  err = '';
end

function net = create_network(param)

% Data
  p = param.inputs;
  t = param.targets;
  if iscell(p), p = cell2mat(p); end
  if iscell(t), t = cell2mat(t); end

  % Dimensions
[R,Q] = size(p);
[S,Q] = size(t);

% Architecture
net = network(1,2,[1;0],[1;0],[0 0;1 0],[0 1]);     %构建自定义神经网络

% Simulation
net.inputs{1}.size = R;   %输入维数R
% 构造RBF
net.inputWeights{1,1}.weightFcn = 'dist';   %输入与第一层连接权重，dist欧几里得函数
net.layers{1}.netInputFcn = 'netprod';    %输入函数：计算输入*权重求和的结果
net.layers{1}.transferFcn = 'radbas';       %传递函数（激活函数）
  
net.layers{1}.size = Q;
net.layers{2}.size = S;
net.layers{2}.transferFcn = 'compet';   %竞争函数
net.outputs{2}.exampleOutput = t;    

% Weight and Bias Values
% net.b{1} = zeros(Q,1)+sqrt(-log(.5))/param.spread;   %节点偏置
zita = 0.3;  % 节点偏置系数
sigma_z = 5.79;  % 粗糙度
sigma = zita * sigma_z;
net.b{1} = zeros(Q,1)+sigma;
net.iw{1,1} = p';    %RBF输入权重（节点数据中心）
net.lw{2,1} = t;     %层之间权重
end

测试重力匹配结果如下图所示，其中红色为真实轨迹，黑色为惯性导航系统指示轨迹，蓝色为PNN匹配的轨迹：