褶积方法制作合成地震记录c++

地震褶积方法制作合成地震记录

包括，（1）读取相模型，设置每种相的密度和速度，（2）计算反射系数，添加噪音，（3）设置子波，（4）进行褶积计算。具体的代码如下

void syntheticSeis(const string& faciesFileName, const string&synseisFileName, 
    vector<tuple<int, double, double>>faciesDenVelocity,
    double dt, double fm, int convLength)
{
    const double PI = 4 * (4 * atan(1.0f / 5) - atan(1.0f / 239));
    IModel3D<int> facies;
    facies.loadNumpy(faciesFileName);
    int ni = facies.getNi();
    int nj = facies.getNj();
    int nk = facies.getNk();
    auto faciesList = facies.getValueVec();
   
    for (auto item : faciesList) {
        bool find = false;
        for (auto fv : faciesDenVelocity) {
            if (get<0>(fv) == item) {
                find = true;
                break;
            }
        }
        if (find == false) {
            std::cout << "fail to find facies  " << item << 
                " in facies_velocity_set" << std::endl;
            return;
        }
    }
    
    // define the velocity of each facies
    auto impedance  = IModel3D<float>(nj,ni,nk,"veclocity");
    const int* faciesPtr = facies.grid().data();
    float* impPtr = impedance.grid().data();
    random_device dev;
    default_random_engine rnd(dev());
    uniform_real_distribution<double> u(0.95, 1);
    for (int i = 0; i < ni * nj * nk; i++) {
        int faciesV = faciesPtr[i];
        for (const auto& item : faciesDenVelocity) {
            if (get<0>(item) == faciesV) {
                //impPtr[i] = get<1>(item)* get<2>(item);
                impPtr[i] = get<1>(item) * get<2>(item) * u(rnd);  //with noise
                break;
            }
        }
    }

    // define the reflection coefficient
    auto refCoef = IModel3D<float>(nj, ni, nk, "reflectionCoefficient");
    for (int j = 0; j < nj; j++) {
        for (int i = 0; i < ni; i++) {
            refCoef.setValue(j, i, nk - 1, 0);
            for (int k = nk - 2; k >= 0; k--) {
                float imp1 = impedance.getValue(j, i, k + 1);
                float imp2 = impedance.getValue(j, i, k);
                float coef = (imp2 - imp1) / (imp2 + imp1);
                refCoef.setValue(j, i, k, coef);
            }
        }
    }
    float vMin = FLT_MAX, vMax = -FLT_MAX;
    refCoef.getMinMax(vMin, vMax);
    std::cout << "refcoef vmin= " << vMin << " vmax= " << vMax << std::endl;

    // calculate syntheic seismic
    int n = convLength;
    auto rickerWave = [=](int kk)->float {
        return (1 - 2 * pow(PI * fm * dt*kk, 2)) * exp(-pow(PI * fm * dt*kk, 2));
    };

    // define the synthetic seismic trace
    auto synSeis = IModel3D<float>(nj, ni, nk, "syntheticSeis");
    for (int j = 0; j < nj; j++) {
        for (int i = 0; i < ni; i++) {
            for (int k = nk - 1; k >= 0; k--) {
                // seismic convolution
                double trace = 0.0;
                for (int t = -int(n / 2); t<int(n / 2); t++) {
                    if (t + k < 0 || t + k >= nk)
                        continue;
                    else
                    {
                        double riker = rickerWave(t);
                        double coef = refCoef.getValue(j, i, k + t);
                        double value = riker * coef;
                        trace += value;
                    }
                }
                synSeis.setValue(j, i, k, trace);
            }
        }
    }
    synSeis.saveNumpy(synseisFileName);
    vMin = FLT_MAX;
    vMax = -FLT_MAX;
    synSeis.getMinMax(vMin, vMax);
    std::cout <<"create syntheitc seismic"<<synseisFileName<<  " synseis vmin= " << vMin << " vmax= " << vMax << std::endl;

}

效果如下