bostick反演

1.bostick反演的基本公式

$$H=\sqrt{\frac{{\rho }_a}{\omega \mu }}$$

$$\rho (H)={\rho }_a(\omega )(\frac{\pi }{2\phi (\omega )}-1)$$

2.代码实现

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
# 设置中文显示
font = {'family':'SimHei'}
plt.rcParams["axes.unicode_minus"] = False

# 定义正演函数
def mt1d_fwd(rho, h):
    mu = (4e-7)*np.pi
    T = np.logspace(-3, 4, 40)
    m = len(rho) # 层数
    
    k = np.zeros((len(rho), len(T)), dtype = np.complex)
    for i in range(len(rho)):
        k[i] = np.sqrt((-1j*2*np.pi*mu)/(T*rho[i]))
        
    Z = -(2j*mu*np.pi)/(T*k[m-1]) # python 的索引从 0 开始
    layer = np.arange(m-1, 0, -1) # np.arange() 不包含 stop 参数
    for n in layer:
        A = -(1j*mu*2*np.pi)/(T*k[n-1])
        B = np.exp(-2*k[n-1]*h[n-1])
        Z = A*(A*(1-B)+Z*(1+B))/(A*(1+B)+Z*(1-B))
    
    rho_a = (T/(mu*2*np.pi))*(np.abs(Z)**2)
    phase = -np.arctan(Z.imag/Z.real)*180/np.pi
    return [rho_a, phase]

# 定义 bostick 反演函数
def bostick(rho_a, ph):
    mu = np.pi*(4e-7)
    T = np.logspace(-3, 4, 40)
    D = np.sqrt((rho_a*T)/(2*np.pi*mu))
    rho_h = rho_a*(180/(2*ph)-1)
    return D, rho_h

# 设置理论模型
T = np.logspace(-3, 4, 40)
depth = np.array([1, 500, 2000, 100000])
rho = [10, 200, 10, 10]
# 使用正演得出的反演计算的原始数据
rho_a, ph = mt1d_fwd([10, 200, 10], [500, 2000])

# 反演计算
D, rho_h = bostick(rho_a, ph)

# 绘制反演结果图
plt.figure()
plt.plot(depth/1000, rho, drawstyle = 'steps-post',
        label = '原始模型')
plt.plot(D/1000, rho_h, drawstyle = 'steps-post',
        label = '反演模型')
plt.legend(prop = font)
plt.xlabel('深度(km)', font)
plt.ylabel(r'电阻率($\Omega\cdot m$)', font)
plt.loglog()
plt.savefig('反演结果.png', dpi = 300)

#绘制反演对比图
plt.figure()
h = np.diff(D)  #计算出bostick反演的层厚度
rho_a1, ph1 = mt1d_fwd(rho_h, h) # 将反演结果进行正演
# 绘制电阻率响应模型
plt.subplot(211)
plt.plot(T, rho_a, '*', label = '原始模型响应')
plt.plot(T, rho_a1, label = '反演模型响应')
plt.xlabel('T(s)')
plt.ylabel(r'$\rho_a(\Omega\cdot m)$')
plt.semilogx()
plt.legend(prop = font)
# 绘制相位响应模型
plt.subplot(212)
plt.plot(T, ph, '*', label = '原始模型响应')
plt.plot(T, ph1, label = '反演模型响应')
plt.xlabel('T(s)')
plt.ylabel('Phase(°)')
plt.semilogx()
plt.legend(prop = font, loc=4)
plt.savefig('反演响应数据对比')

plt.show()

3.反演结果

反演结果：

反演数据对比

参考资料：柳建新，童孝忠，MATLAB在地球物理中的应用。