一维大地电磁测深模拟退火反演

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%一维大地电磁测深模拟退火法反演
%by xiejingtao
%2016.05.21
%模拟退火法的优点：质量高；初值鲁棒性强；简单、通用、易实现
%缺点：由于要求较高的初始温度，较慢的降温速率、较低的终止温度，以及各温度下足够多次的抽样率，因此优化过程较长。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc
clear all
close all

%初始化参数
u=4*pi*1e-7;
f=[1500.000
1500.000
1000.000
800.000    
560.000
480.000
400.000
320.000
240.000
160.000
120.000
80.0000
60.0000
40.0000
30.0000
20.0000
15.0000
10.0000
7.50000
5.00000
3.75000
2.50000
1.87500
1.25000
.937500
.625000
.468800
.312500
.234400
.156300
.117200
.781300E-01
.585900E-01
.390600E-01
.293000E-01
.195300E-01
.146500E-01
.977000E-02
.733000E-02
.488000E-02
.366000E-02
.244000E-02
.183000E-02
.122000E-02
.920000E-03
.610000E-03
.460000E-03];
f=f‘;
nfreq=length（f）;


%给定初始模型
m0=[300;200;100;100;100];

%给定模型取值范围
area=[01000;01000;01000;01000;01000];

for j=1:5
    while m0（j1）area（j2）
        if m0（j1）            m0（j1）=2*area（j1）-m0（j1）;
        elseif m0（j1）>area（j2）
                m0（j1）=2*area（j2）-m0（j1）;
        end
    end
end

%求原始数据视电阻率
    
%数据参数
m1=[100;20;300];
rho1=m1（1:2:）;
h1=m1（3:）;
n1=length（rho1）;
    
%求视电阻率
for ifreq=1:nfreq
    w=2*pi*f（ifreq）;
    k=sqrt（i*w*u/rho1（n1））;
    Z=w*u/k;
    for j=n1:-1:2
        k=sqrt（i*w*u/rho1（j-1））;
        Z0=w*u/k;
        e1=1-exp（2*i*k*h1（j-1））;
        e2=1+exp（2*i*k*h1（j-1））;
        Z=Z0*（Z0*e1+Z*e2）/（Z0*e2+Z*e1）;
    end
    rhoas（ifreq）=（abs（Z））^2/（w*u）;
end

for iter=1:20000
    
    %m0模型视电阻率
    rho0=m0（1:3:）;
    h0=m0（4:5:）;
    n0=length（rho0）;
    
    %求视电阻率
    for ifreq=1:nfreq
        w=2*pi*f（ifreq）;
        k=sqrt（i*w*u/rho0（n0））;
        Z=w*u/k;
        for j=n0:-1:2
            k=sqrt（i*w*u/rho0（j-1））;
            Z0=w*u/k;
            e1=1-exp（2*i*k*h0（j-1））;
            e2=1+exp（2*i*k*h0（j-1））;
            Z=Z0*（Z0*e1+Z*e2）/（Z0*e