直流电法正演模拟

function T = Transfer_Fun(m)

global rho;

global h;

N = size(rho,2);

T = rho(N);

for i = N-1:-1:1

A = 1-exp(-2*m*h(i));

B = 1+exp(-2*m*h(i));

T = rho(i)*(rho(i)*A+T*B)/(rho(i)*B+T*A); end

T = T*m;

function z1 = hankel1_47(fun,B)

a = -3.05078187595e+00;

s = 1.10599010095e-01;

wt1 = [3.17926147465e-06

-9.73811660718e-06

1.64866227408e-05

-1.81501261160e-05

1.87556556369e-05

-1.46550406038e-05

1.53799733803e-05

-6.95628273934e-06

1.41881555665e-05

3.41445665537e-06

2.139********e-05

2.34962369042e-05

4.84340283290e-05

7.33732978590e-05

1.27703784430e-04

2.0812*******e-04

3.49803898913e-04

5.79107814687e-04

9.65887918451e-04

1.60401273703e-03

2.66903777685e-03

4.43111590040e-03

7.35631696247e-03

1.21782796293e-02

2.010********e-02

3.30096953061e-02

5.37143591532e-02

8.60516613299e-02

1.34267607144e-01

2.00125033067e-01

2.74027505792e-01

3.181********e-01

2.41655667461e-01

-5.40549161658e-02

-4.46912952135e-01

-1.92231885629e-01

5.52376753950e-01

-3.57429049025e-01

1.41510519002e-01

-4.61421935309e-02

1.48273761923e-02

-5.07479209193e-03

1.83829713749e-03

-6.67742804324e-04

2.21277518118e-04

-5.66248732755e-05

7.88229202853e-06];

z1 = 0;

for i = 1:size(wt1,1)

lamda = (1/B)*10^(a+(i-1)*s);

z1 = z1+feval(fun,lamda)*wt1(i);

end

z1 = z1/B;

function [rho_s] = DC1D_Forward(rho1,h1) global rho;

global h;

rho = rho1;

h = h1;

r = logspace(1,4,40);

for i = 1:size(r,2)

z(i) = hankel1_47('Transfer_Fun',r(i));

rho_s(i) = r(i)^2*z(i);

end