周报(7.29-8.4)

XLyu233

已于 2024-08-06 14:54:38 修改

阅读量508

点赞数 15

文章标签：学习笔记

于 2024-08-06 14:54:23 首次发布

本文链接：https://blog.csdn.net/XLyu233/article/details/140955985

版权

周报(7.29-8.4)

本周工作

学习了地震正演技术，但仍需加强学习

一维声波方程:
$\frac{∂^2u}{∂t^2}= a^2\frac{∂^2u}{∂x^2}+ f(x,t)$
其中： $u$ 是一维弦在纵轴上的波动位移， $t$ 是时间， $a$ 是加速度， $x$ 是弦的位置， $a^2=\frac{T}{\rho}$ ， $T$ 是每微分弦端点所被施加的力，经计算T是恒量， $\rho$ 是弦密度， $f(x,t)=\frac{1}{\rho}F(x,t)$ ， $F (x, t)$ 是某时施加在某处的外力。将一维方程扩展到二维，需要扩展对 $x$ 的偏导到求 $x, y$ 的偏导，再结合前述条件，则二维声波自由震动方程为
$\frac{1}{v^2}\frac{∂^2U}{∂t^2}=\frac{∂^2U}{∂x^2}+\frac{∂^2U}{∂z^2}$
在现实宏观情况下，声波的震动可无限细分，但在计算机模拟中，我们必须分时间区域来离散的计算。为此引入网格差分，能得到用于计算机模拟的式子：
$\frac{1}{v^2}\frac{U^{k+1}_{i,j}-2U^k_{i,j}+U^{k-1}_{i,j}}{\Delta t^2}=\frac{U^k_{i+1,j}-2U^k_{i,j}+U^k_{i-1,j}}{\Delta x^2}+\frac{U^k_{i,j+1}-2U^k_{i,j}+U^k_{i,j-1}}{\Delta z^2}$
由式子，编写代码如下:
```
import scipy.io
import numpy as np
import matplotlib.pyplot as plt
import math
import tqdm
# 模拟时长
EndTime = 0.5
# 间隔时间
delta_t = 0.0005
# 空间步长(m)
delta_x = 6
# 横向网格数
CNX = 301
# 纵向网格数
CNZ = 301
# 震源位置（首先设置到中心）
Sx=(CNX+1)//2
Sz=(CNZ+1)//2
# 震源频率
f0 = 30
# 速度模型
# v = seg_real_v_data_mat.copy()
v = np.ones((CNZ, CNX),dtype=np.int32)*1500
Unow = np.zeros((CNZ, CNX))
Uprev = np.zeros((CNZ, CNX))
Unext = np.zeros((CNZ, CNX))
data = []
times = []
time = 0
while time<EndTime:
    times.append(time)
    time+=delta_t
# while time<=EndTime:
# for time in range(0,EndTime,delta_t):
for time in tqdm.tqdm(times):
    for i in range(1, CNZ-1):
        for j in range(1, CNX-1):
            A = (-2*Unow[i,j]+Unow[i+1,j]+Unow[i-1,j])/(delta_x**2)
            B = (-2*Unow[i,j]+Unow[i,j+1]+Unow[i,j-1])/(delta_x**2)
            Unext[i,j] = 2*Unow[i,j]-Uprev[i,j]+(v[i,j]**2)*(A+B)*(delta_t**2)
    Unext[Sz,Sx]=5.76*(f0**2)*(1-16*((0.6*f0*time-1)**2))*math.exp(-8*(0.6*f0*time-1)**2)
    
    data.append(Unext.copy())
    
    Uprev = Unow.copy()
    Unow = Unext.copy()
    time+=delta_t
plt.imshow(data[399], cmap="gray")
```
结果如下:

这是模拟大小为(301, 301)，速度为1500m/s，震源在中心的地层中，以6m、0.0005s的精度模拟0.5s后的波形传播结果的图像。虽然numpy底层使用c++进行运算，但在python中频繁进行的各项调用也拖慢了代码运行的速度，使得模拟时间达到了5min，因此尝试用c语言重写代码：
```
#define PY_SSIZE_T_CLEAN
#include <Python.h>
static PyObject *BasicSpawn(PyObject *self, PyObject *args, PyObject *keywds)
{
    // 模拟时长
    float EndTime = 0.5F;
    // 间隔时间
    float Delta_t = 0.0005F;
    // 空间步长(m)
    float Delta_x = 6;
    // 横向网格数
    int CNX = 301;
    // 纵向网格数
    int CNZ = 301;
    // 震源位置（首先设置到中心）
    int Sx = (CNX + 1) / 2;
    int Sz = 0;
    // 震源频率
    int F0 = 30;
    // 速度模型
    // int v = 1500;
    int *V = NULL;
    int VFree = 0;
    PyObject * PyV = NULL;
    static char *kwlist[] = {"EndTime", "Delta_t", "Delta_x", "CNX", "CNZ", "Sx", "Sz", "F0", "V", NULL};
    int ok = PyArg_ParseTupleAndKeywords(args, keywds, "|fffiiiiiY", kwlist, &EndTime, &Delta_t, &Delta_x, &CNX, &CNZ, &Sx, &Sz, &F0, &PyV);
    if (!ok)
    {
        return NULL;
    }
    if (PyV == NULL){
        V = malloc(sizeof(int)*CNX*CNZ);
        VFree = 1;
        for (int i=0;i<CNZ;i++){
            for (int j=0;j<CNX;j++){
                V[i*CNZ+j] = 1500;
            }
        }
    } else if (PyByteArray_CheckExact(PyV)){
        // 长度不正确
        if (PyByteArray_Size(PyV)!=CNZ*CNX*sizeof(int32_t)){
            PyErr_SetString(PyExc_ValueError, "Input shape not equal to CNX*CNZ.");
            return NULL;
        }
        V = (int *)PyByteArray_AsString(PyV);
    }
    // 模拟次数
    int count = EndTime / Delta_t+1;

    float *Unow, *Uprev, *Unext, *Ures;
    Unow = malloc(sizeof(float) * CNZ * CNX);
    Uprev = malloc(sizeof(float) * CNZ * CNX);
    Unext = malloc(sizeof(float) * CNZ * CNX);
    
    memset(Unow, 0, sizeof(float) * CNZ * CNX);
    memset(Uprev, 0, sizeof(float) * CNZ * CNX);
    memset(Unext, 0, sizeof(float) * CNZ * CNX);

    Ures = malloc(sizeof(float) * count * CNX);
    float time = 0.0;
    for (int timeCount = 0; timeCount < count; time += Delta_t, timeCount++)
    {
        for (int i = 1; i < CNZ - 1; i++)
        {
            for (int j = 1; j < CNX - 1; j++)
            {
                float A = (-2.0 * Unow[i * CNX + j] + Unow[(i + 1) * CNX + j] + Unow[(i - 1) * CNX + j]) / pow(Delta_x, 2.0);
                float B = (-2.0 * Unow[i * CNX + j] + Unow[i * CNX + j + 1] + Unow[i * CNX + j - 1]) / pow(Delta_x, 2.0);
                Unext[i * CNX + j] = 2.0 * Unow[i * CNX + j] - Uprev[i * CNX + j] + pow(V[i*CNX+j], 2.0) * (A + B) * pow(Delta_t, 2.0);
            }
        }
        Unext[Sz * CNX + Sx] = 5.76 * pow(F0, 2) * (1 - 16 * pow(0.6 * F0 * time - 1, 2)) * exp(-8.0 * pow(0.6 * F0 * time - 1, 2));

        memcpy(Uprev, Unow, sizeof(float) * CNZ * CNX);
        memcpy(Unow, Unext, sizeof(float) * CNZ * CNX);

        for (int i = 0; i< CNX;i++){
            printf("count:%d, timeCount:%d, i:%d, time:%f\n",count, timeCount, i, time);
            Ures[timeCount*CNX+i] = Unow[0*CNX+i];
        }
    }
    printf("123\n");
    free(Uprev);
    free(Unext);
    free(Unow);
    printf("????\n");
    PyObject *res = PyByteArray_FromStringAndSize((char *)Ures, sizeof(float) * count * CNX);
    free(Ures);

    if (VFree)
        free(V);
    return res;
}

static PyMethodDef PySEGSpawnMethods[] = {
    {"BasicSpawn", BasicSpawn, METH_VARARGS | METH_KEYWORDS, "a BasicSpawn function"},
    {NULL, NULL, 0, NULL}};

static struct PyModuleDef PySEGSpawnModule = {
    PyModuleDef_HEAD_INIT,
    "PySEGSpawn",
    NULL,
    -1,
    PySEGSpawnMethods};
PyMODINIT_FUNC PyInit_PySEGSpawn(void)
{
    return PyModule_Create(&PySEGSpawnModule);
}
```
为了方便python进行调用，我对c语言代码进行特殊处理，该代码将被编译为python的一个库，可被python代码调用，并且速度大幅提升，仅需十余秒就可以完成模拟。可以通过传入速度图像字节串和其他参数控制生成，结果类型也为字节串，可以直接载入到numpy进行显示。调用代码如下，结果与python代码相同。
```
import PySEGSpawn
import numpy as np
import matplotlib.pyplot as plt
test = PySEGSpawn.BasicSpawn(EndTime=0.5, Delta_t=0.0005)
data = np.frombuffer(test,dtype=np.float32)
plt.imshow(data.reshape(301,301), cmap="gray")
```
学习DD-Net网络

下周任务

复现DD-Net网络代码。
选择具体方向，开展研究。
继续学习地震反演原理，深度学习原理。

XLyu233

关注

15
点赞
踩
12

收藏

觉得还不错? 一键收藏
0
评论
周报(7.29-8.4)

学习了地震正演技术，但仍需加强学习∂t2∂2ua2∂x2∂2ufxtu是一维弦在纵轴上的波动位移，t是时间，a是加速度，x是弦的位置，a2ρTT是每微分弦端点所被施加的力，经计算T是恒量，ρ是弦密度，fxtρ1FxtFxt是某时施加在某处的外力。将一维方程扩展到二维，需要扩展对x的偏导到求xy的偏导，再结合前述条件，则二维声波自由震动方程为v21∂t2∂2U∂x2∂2U。
复制链接

扫一扫