sphere_inverse.m

这是一个用MATLAB实现的磁法球体相关性自动反演的小程序

%  现在生成最终需要拟合的标准的数据，即作为在实际工作中的采集的数据，
%  在本程序中，磁化率作为已知的参数使用
[dainihe_X,dainihe_Y,dainihe_data]=datacreate(0,0,150,100,10);
%  对以上函数的调用的说明：
%  0，0，150是待拟合或者说是待寻找的目标模型的三维坐标
%  100是一条侧线观测点的数据个数
%  10是观测点的间距

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  现在生成初始的任意八个比较的模型的数据，
%  对于参数的意义同上述一样
%  生成数据放在相关性的计算和比较的内部循环中

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  现在计算以上八个数据与采集数据的相关性
[data_length,data_width]=size(dainihe_data);
data_total=data_length*data_width;
p_max=0; %相关性的最大值
k1=0;
for k=-1:2:1
    k1=k1+1;
    i1=0;
 for i=-1:2:1
     i1=i1+1;
     j1=0;
    for j=-1:2:1
        j1=j1+1;
 % 生成比较的数据

    [dainihe_Xi,dainihe_Yi,dainihe_datai]=datacreate(0+400*i,0+400*j,150+50*k,100,10);
    X_now=0+400*i;
    Y_now=0+400*j; 
    Z_now=150+50*k;
 %  下面进行比较
    phrase1=0;                                  
    phrase2=0;
    phrase3=0;
    for n=1:data_total
          
        phrase1=phrase1+dainihe_datai(n)*dainihe_data(n);
        phrase2=phrase2+dainihe_datai(n)*dainihe_datai(n);
        phrase3=phrase3+dainihe_data(n)*dainihe_data(n);
        
       end

    p(i1)=phrase1/sqrt(phrase2*phrase3);
%  下面求八个相关值的最大值，并保存相应的场源的参数备用；
        if  p_max<p(i1)
        p_max=p(i1);
        X_next=X_now;
        Y_next=Y_now;
        Z_next=Z_now;
        end
    end
 end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  下面用以上的比较作为初值进行后续的比较
break_flag=0;
%   break_flag 是得到满意的结果后的退出标志
for iii=1:150
%  现在设定最多进行150次比较，以防止进入无限的死循环
   k1=0;
for k=-1:2:1
    k1=k1+1;
    i1=0;
 for i=-1:2:1
     i1=i1+1;
     j1=0;
    for j=-1:2:1
        j1=j1+1;
% 生成比较的数据

    [dainihe_Xi,dainihe_Yi,dainihe_datai]=datacreate(X_next+10*i,Y_next+10*j,Z_next+10*k,100,10);
    X_now=X_next+10*i;
    Y_now=Y_next+10*j; 
    Z_now=Z_next+10*k;
% 下面进行比较
    phrase1=0;                                  
    phrase2=0;
    phrase3=0;
    for n=1:data_total
          
        phrase1=phrase1+dainihe_datai(n)*dainihe_data(n);
        phrase2=phrase2+dainihe_datai(n)*dainihe_datai(n);
        phrase3=phrase3+dainihe_data(n)*dainihe_data(n);
        
       end

    p(i1)=phrase1/sqrt(phrase2*phrase3);
%  下面求八个相关值的最大值，并保存相应的场源的参数备用；
        if  p_max<p(i1)
        p_max=p(i1);
        X_next=X_now;
        Y_next=Y_now;
        Z_next=Z_now;
        eventual_datai=dainihe_datai;
        end
%   如果某次比较的相关系数的值大于等于0.9999，则推出循环停止搜索，将次场源的参数作为最后的结果
                  if p_max>=0.9999

                         break;  
                         break_flag=1;     

                  end
       end
                                   if  break_flag==1
                                   break;
                                   end


     end
                                   if  break_flag==1
                                   break;
                                   end
   end    
                                   if  break_flag==1
                                   break;
                                   end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  下面完成最终的结果的可视化显示
%   根据最初的讨论，下面我们仅仅绘制一条测线的数据作为比较

figure;
plot(eventual_datai(10,:));
%  以上取与X 轴平行的测线
hold on
plot(dainihe_data(10,:),'r');
hold off
legend('运算结果','原始结果');