In ScottHMcKean/duvernay_geomech_model: A Geostatistical Model of Geomechanical Properties in the Duvernay Play

library(duvernaygeomechmodel)
set.seed(1)

Load datasets, encode factors, and set brittleness with residual data

mldf_nores = read_csv("../data/merged_tests_no_residual.csv") %>%
  mutate(horiz_orient = ifelse(orientation == 'horizontal', 1, 0)) %>%
  mutate(shale_group = ifelse(group == 'shale', 1, 0)) %>%
  mutate(poisson_ratio = ifelse(poisson_ratio > 0.5, 0.5, poisson_ratio)) %>%
  as.data.frame()

mldf_wres = read_csv("../data/merged_tests_w_residual.csv") %>%
  mutate(residual_stress_dev = ifelse(residual_stress_dev < 0, 0, residual_stress_dev)) %>%
  mutate(brittleness = (peak_stress_dev - residual_stress_dev)/peak_stress_dev) %>%
  mutate(horiz_orient = ifelse(orientation == 'horizontal', 1, 0)) %>%
  mutate(shale_group = ifelse(group == 'shale', 1, 0)) %>%
  mutate(poisson_ratio = ifelse(poisson_ratio > 0.5, 0.5, poisson_ratio)) %>%
  mutate(brittleness = ifelse(brittleness > 1, 1, brittleness)) %>%
  as.data.frame()

feat_set = c('pred_vp', 'bulk_density', 'horiz_orient', 'shale_group')
feat_pref = 'vp_rho'

generate_models(mldf_nores, feats = feat_set, 
  targets = c("youngs_modulus", "poisson_ratio"), feat_prefix = feat_pref)

generate_models(mldf_wres, feats = feat_set, 
  targets = c("brittleness"), feat_prefix = feat_pref)

Get the Wave Equation results - this is a bit of a 'model' hack as we simply optimize the correction factor and bootstrap to get uncertainty. This provides the fairest comparison to statistical models.

n = 1000
mc_cv = crossv_mc(mldf_nores, n, test=0.9)

dfs <- vector(mode = "list", length = n)
for (iter in 1:n){
  train = mldf_nores[mc_cv$train[[iter]]$idx,]
  test = mldf_nores[mc_cv$test[[iter]]$idx,]
  opt_pr_correction = optimize(pr_correction_loss, c(-1,1), test)$minimum
  train_pr_pred = wave_pr_w_correction(train, opt_pr_correction)
  test_pr_pred = wave_pr_w_correction(test, opt_pr_correction)
  train_pr_mae = mean(abs(train_pr_pred - train$poisson_ratio))
  train_pr_rmse = sqrt(mean((train_pr_pred - train$poisson_ratio)^2))
  test_pr_mae = mean(abs(test_pr_pred - test$poisson_ratio))
  test_pr_rmse = sqrt(mean((test_pr_pred - test$poisson_ratio)^2))
  opt_ym_correction = optimize(ym_correction_loss, c(-100,100), test)$minimum
  train_ym_pred = wave_ym_w_correction(train, opt_ym_correction)
  test_ym_pred = wave_ym_w_correction(test, opt_ym_correction)
  train_ym_mae = mean(abs(train_ym_pred - train$youngs_modulus))
  train_ym_rmse = sqrt(mean((train_ym_pred - train$youngs_modulus)^2))
  test_ym_mae = mean(abs(test_ym_pred - test$youngs_modulus))
  test_ym_rmse = sqrt(mean((test_ym_pred - test$youngs_modulus)^2))
  dfs[[iter]] = data.frame(
    iter = c(iter,iter, iter, iter),
    correction = c(opt_pr_correction, opt_pr_correction, opt_ym_correction, opt_ym_correction),
    mae = c(train_pr_mae, test_pr_mae, train_ym_mae, test_ym_mae),
    rmse = c(train_pr_rmse, test_pr_rmse, train_ym_rmse, test_ym_rmse),
    set = c('train', 'test', 'train', 'test'),
    model_type = c('WAVE', 'WAVE','WAVE', 'WAVE'),
    target = c("pr", "pr", "ym", "ym"),
    feat_set = c("vp_vs_rho", "vp_vs_rho", "vp_vs_rho", "vp_vs_rho")
  )
}

wave_results = bind_rows(dfs)
write_rds(wave_results, '../output/statistical_models/wave/resample_results.rds')

feat_set = c('pred_vp', 'bulk_density', 'pred_vs', 'horiz_orient', 'shale_group')
feat_pref = 'vp_vs_rho'

generate_models(mldf_nores, feats = feat_set, 
  targets = c("youngs_modulus", "poisson_ratio"), feat_prefix = feat_pref)

generate_models(mldf_wres, feats = feat_set, 
  targets = c("brittleness"), feat_prefix = feat_pref)

feat_set = c('deviatoric_stress', 'pred_vp', 'bulk_density', 'horiz_orient', 'shale_group')
feat_pref = 'dev_vp_rho'

generate_models(mldf_nores, feats = feat_set, 
  targets = c("youngs_modulus", "poisson_ratio"), feat_prefix = feat_pref)

generate_models(mldf_wres, feats = feat_set, 
  targets = c("brittleness"), feat_prefix = feat_pref)

feat_set = c('deviatoric_stress', 'pred_vs', 'pred_vp', 'bulk_density', 'horiz_orient', 'shale_group')
feat_pref = 'dev_vp_vs_rho'

generate_models(mldf_nores, feats = feat_set, 
  targets = c("youngs_modulus", "poisson_ratio"), feat_prefix = feat_pref)

generate_models(mldf_wres, feats = feat_set, 
  targets = c("brittleness"), feat_prefix = feat_pref)

feat_set = c('deviatoric_stress', 'confining_pressure', 'pred_vp', 'bulk_density', 'horiz_orient', 'shale_group')
feat_pref = 'dev_conf_vp_rho'

generate_models(mldf_nores, feats = feat_set, 
  targets = c("youngs_modulus", "poisson_ratio"), feat_prefix = feat_pref)

generate_models(mldf_wres, feats = feat_set, 
  targets = c("brittleness"), feat_prefix = feat_pref)

feat_set = c('deviatoric_stress', 'confining_pressure', 'pred_vs', 'pred_vp', 'bulk_density', 'horiz_orient', 'shale_group')
feat_pref = 'dev_conf_vp_vs_rho'

generate_models(mldf_nores, feats = feat_set, 
  targets = c("youngs_modulus", "poisson_ratio"), feat_prefix = feat_pref)

generate_models(mldf_wres, feats = feat_set, 
  targets = c("brittleness"), feat_prefix = feat_pref)

ScottHMcKean/duvernay_geomech_model documentation built on Sept. 12, 2022, 10:08 a.m.