探地雷达均匀与随机介质波场模拟（基于有限差分法）

探地雷达均匀与随机介质波场模拟（基于有限差分法）

---

---

前言

最近开一个新系列，是关于探地雷达的波场正演模拟。毕业那段时间弄了很多地震波场正演模拟的东西，感觉用神经网络模拟确实比较玄学，彻底击碎了我想利用深度学习的方法做正演模拟的幻想，暑假期间痛定思痛去学了一下有限差分法和波场模拟的相关内容，结合探地雷达正演，做了一些有趣的正演模拟。
我把相关代码都放在博客里面了，觉得好用的可以拿去用，正演代码是基于冯德山的《探地雷达数值模拟及程序实现》一书中的matlab程序改进而来，我只是将这些代码改成Python语言并使用Torch实现加速功能；随机介质的构造算法来源于郭士礼的博士论文《基于随机介质的高速公路探地雷达检测理论研究》，代码全部由我个人编写。
本篇博客重在展示结果，其中有些过程现阶段我认为还是很难描述明白，以后可能分期详细描述

---

一、自相关函数的构造

强烈建议去重新阅读《基于随机介质的高速公路探地雷达检测理论研究》中有关随机介质构造的理论部分，这里不多做赘述。我们直接上代码。

import numpy as np
import matplotlib.pyplot as plt
import os
os.environ['KMP_DUPLICATE_LIB_OK']="True"
#建造一个计算自相关函数的类，需要可以在这里随意建造
#传入的参数有水平方向，竖直方向的分割数组，自相关角度，粗糙因子，水平和竖直方向的自相关角度
class genereate_autocor():
    def __init__(self,x,z,theta,r,a,b):
        self.x = x
        self.z = z
        self.a = a
        self.b = b
        self.theta = theta
        self.r = r
        self.X,self.Z = np.meshgrid(x,y)
        self.A = (np.cos(self.theta)*self.X + np.sin(self.theta)*self.Z)/self.a
        self.B = (np.sin(self.theta)*self.X - np.cos(self.theta)*self.Z)/self.b
#高斯模型       
    def gauss_model(self):
        phi = np.exp(-(self.A**2+self.B**2))
        return phi
#指数模型
    def exp_model(self):
        phi = np.exp(-np.sqrt(self.A**2+self.B**2))
        return phi
#混合模型
    def mix_model(self):
        phi = np.exp(-(self.A**2+self.B**2)**(1/(1+self.r)))
        return phi    

#建造一个生成随机扰动模型的类，输入的参数有空间的划分网格数，自相关函数
class genereate_SRD():
    def __init__(self,space_step,auto_func):
        self.space_step = space_step
        self.auto_func = auto_func

    def srd(self):
        #先将自相关函数进行反向平移
        fft_auto = np.fft.ifftshift(self.auto_func)
        #做傅里叶变换
        fft_auto = np.fft.fft2(fft_auto)
        #生成随机相位谱
        random_model = [np.random.uniform(0,2*np.pi) for i in range(int(self.space_step**2))]
        phase = []
        for i in range(0,len(random_model),1):
            phase.append(complex(np.cos(random_model[i]),np.sin(random_model[i])))
        phase = np.array(phase).reshape((self.space_step,self.space_step))
        #得到随机功率谱
        random_power_spec = abs(fft_auto)*phase
        #做反傅里叶变换
        ifft_random_power_spec = np.fft.ifft2(random_power_spec)
        #得到的随机扰动进行均值和方差的标准化
        std_random_dis = (abs(ifft_random_power_spec) - np.mean(abs(ifft_random_power_spec)))/np.std(abs(ifft_random_power_spec))

        return std_random_dis

#写一个局部区域腐蚀的函数
#骨料eps1,沥青eps2,空隙eps3，骨料率p1,孔隙率p3,局部区域大小r边长
def corrode(eps1,eps2,eps3,p1,p3,r,mat):
    mat = mat.flatten()
    index = np.linspace(1,len(mat),len(mat))
    eps = np.zeros(shape = mat.shape)
    
    mat = mat.reshape((len(mat),1))
    index = index.reshape((len(index),1))
    eps = eps.reshape((len(eps),1))
    
    combine = np.hstack((index,eps,mat))
    combine = combine[np.lexsort(-combine.T)]
    num_p1,num_p3 = 0,0
    for i in range(0,len(mat),1):
        if num_p1/len(mat) <= p1:
            combine[i,1] = eps1
            num_p1 = num_p1+1
    for j in range(-1,-len(mat)-1,-1):
        if num_p3/len(mat) <= p3:
            combine[j,1] = eps3
            num_p3 = num_p3+1
    for k in range(0,len(mat),1):
        if combine[k,1]!=eps1 and combine[k,1]!=eps3:
            combine[k,1] = eps2
    combine = combine[np.lexsort(combine[:,::-1].T)]
    eps,index = combine[:,1].reshape((r,r)),combine[:,0].reshape((r,r))
    return eps  

#最后写一个构造整个介质的程序
def generate_Dp_media(input_gaint_mat,eps1,eps2,eps3,p1,p3,r):
    input_size = len(input_gaint_mat[0,:])
    e = np.zeros(shape = input_gaint_mat.shape)
    for i in range(0,int(input_size/r),1):
        for j in range(0,int(input_size/r),1):
            e0 = corrode(eps1,eps2,eps3,p1,p3,r,input_gaint_mat[i*r:i*r+r,j*r:j*r+r])
            e[i*r:i*r+r,j*r:j*r+r] = e0
    return e
#从这里我们开始演示利用自相关函数构造某种介质的介电常数
#用一下论文里的数据
y = np.linspace(-1.75,1.75,350)
x = np.linspace(-1.75,1.75,350)
t = np.pi/4*0
X,Y = np.meshgrid(x,y)
r = 2.7
a,b = 4.42*0.01,4.62*0.01
#先生成三种自相关函数
obj_1 = genereate_autocor(x,y,t,r,a,b)
gauss = obj_1.gauss_model()
exp = obj_1.exp_model()
mix = obj_1.mix_model()
#再分别生成三种自相关函数下的随机扰动
std_gauss = genereate_SRD(350,gauss).srd()
std_exp = genereate_SRD(350,exp).srd()
std_mix = genereate_SRD(350,mix).srd()
#生成三种随机扰动下的介质介电常数
eps_gauss = generate_Dp_media(std_gauss,