第九届“数维杯”大学生数学建模挑战赛（C题）深度剖析|建模完整过程+详细思路+代码全解析

问题1

问题1的建模过程如下：

设勘探区域为$D$，勘探井位数量为$n$，每个勘探井位的坐标为$(x_i,y_i),i=1,2,...,n$。根据勘探数据，假设该区域内天然气水合物资源的分布范围为一个矩形区域，记为$R=[x_{min},x_{max}]\times [y_{min},y_{max}]$。

根据天然气水合物资源的分布特点，可以采用高斯核函数来估计资源的分布情况，即
$$f(x,y)=\sum _{i=1}^n\frac{1}{2\pi {\sigma }^2}{e}^{-\frac{(x-x_i{)}^2+(y-y_i{)}^2}{2{\sigma }^2}},(x,y)\in D$$
其中，$f(x,y)$表示资源分布的概率密度函数，$\sigma$为高斯核函数的标准差，可以根据具体情况来确定。

根据资源的分布概率密度函数，可以求得资源的累积分布函数为
$$F(x,y)={\int }_{-\infty }^x{\int }_{-\infty }^{y}f(u,v)dudv$$
根据题目给出的勘探数据，可以计算出每个勘探井位对应的天然气水合物资源的概率分布，即每个点$(x_i,y_i)$对应的概率$f(x_i,y_i)$。因此，可以根据累积分布函数的性质来确定天然气水合物资源的分布范围，即求解下面的不等式组
{ F ( x , y ) ≤ 0.95 F ( x , y ) ≥ 0.05 ( x , y ) ∈ R

$$\begin{cases} F(x,y)\le 0.95\\ F(x,y)\ge 0.05\\ (x,y)\in R \end{cases}$$ ⎩ ⎨ ⎧​F(x,y)≤0.95F(x,y)≥0.05(x,y)∈R​
其中，$F(x,y)$为累积分布函数。

解出上述不等式组即可得到天然气水合物资源的分布范围$R$。

根据所给勘探数据，我们可以确定天然气水合物资源的分布范围为勘探区域内的14口井。每口井都有相应的深度信息和测量的孔隙度和天然气水合物饱和度信息。因此，我们可以根据勘探数据中的深度信息和勘探区域的地质信息，确定每口井的坐标，并绘制出勘探区域的地图。

python示例代码实现：

import pandas as pd
import matplotlib.pyplot as plt
import geopandas as gpd

# 读取勘探井位信息数据
well_data = pd.read_excel('well_data.xlsx')

# 绘制勘探井位地图
fig, ax = plt.subplots(figsize=(8, 8))
ax.scatter(well_data['经度'], well_data['纬度'], color='red', s=5)
ax.set_xlabel('经度')
ax.set_ylabel('纬度')
ax.set_title('勘探井位地图')
plt.show()

# 连接勘探井位线段
fig, ax = plt.subplots(figsize=(8, 8))
ax.plot(well_data['经度'], well_data['纬度'], color='red', linewidth=2)
ax.set_xlabel('经度')
ax.set_ylabel('纬度')
ax.set_title('勘探井位连线图')
plt.show()

# 利用地理信息系统（GIS）绘制资源分布范围
gdf = gpd.GeoDataFrame(well_data, geometry=gpd.points_from_xy(well_data['经度'], well_data['纬度']))
fig, ax = plt.subplots(figsize=(10, 10))
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
ax.set_aspect('equal')
world.plot(ax=ax, color='lightgray', edgecolor='black')
gdf.plot(ax=ax, color='red', markersize=10)
plt.show()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

问题2

问题2的建模过程如下：

假设研究区域内的有效厚度、地层孔隙度和饱和度分别为$T(x,y)$、$P(x,y)$和$S(x,y)$，其中$x$和$y$分别表示水平和垂直方向的空间坐标。则在勘探区域内，这些参数的概率分布可以表示为：

有效厚度的概率分布：

$$P(T)=\frac{N_T}{N}$$

其中，$N_T$表示研究区域内有效厚度大于零的点的数量，$N$表示研究区域内总的点的数量。

地层孔隙度的概率分布：

$$P(P)=\frac{N_P}{N}$$

其中，$N_P$表示研究区域内地层孔隙度大于零的点的数量，$N$表示研究区域内总的点的数量。

饱和度的概率分布：

$$P(S)=\frac{N_S}{N}$$

其中，$N_S$表示研究区域内饱和度大于零的点的数量，$N$表示研究区域内总的点的数量。

这些参数的变化规律可以通过统计分析得到，比如通过计算不同区域内的平均值、方差等。

3)给出天然气水合物的概率分布，以及估计天然气水合物资源量。

天然气水合物的概率分布可以表示为：

$$P(H)=P(T)\times P(P)\times P(S)$$

其中，$H$表示天然气水合物的存在概率，$P(T)$、$P(P)$和$P(S)$分别表示有效厚度、地层孔隙度和饱和度的概率分布。

天然气水合物的资源量可以通过计算总的天然气水合物储量得到，即：

$$Q=\sum _{i=1}^N{H}_i\times V_i\times F$$

其中，$N$表示勘探区域内的点的数量，$H_i$表示第$i$个点的天然气水合物的存在概率，$V_i$表示该点的储量，$F$表示产气量因子。

问题3

问题三的建模过程如下：

1）天然气水合物资源的概率分布：
天然气水合物资源的概率分布可以通过计算各个勘探点的水合物饱和度和有效厚度的乘积来得到。假设勘探点数为n，勘探点i的水合物饱和度为$S_i$，有效厚度为h_i，则该点的天然气水合物资源量为$V_i=S_i\ast h_i$。根据中心极限定理，当样本量较大时，各个勘探点的天然气水合物资源量近似为正态分布，因此可以对各个勘探点的资源量进行求和，得到整个研究区域内天然气水合物资源量的概率分布，即：
$$P(V)=\sum _{i=1}^n\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{(V-V_i{)}^2}{2{\sigma }^2}}$$
其中，$\sigma$为资源量的标准差，可以通过样本数据计算得到。

2）估计天然气水合物资源量：
天然气水合物资源量的估计可以通过计算平均资源量和总资源量来得到。平均资源量为各个勘探点的资源量的平均值，即：
$$\bar{V}=\frac{1}{n}\sum _{i=1}^n{V}_i$$
总资源量为各个勘探点的资源量的总和，即：
$$V_{total}=\sum _{i=1}^n{V}_i$$

3）具体计算：
假设勘探点数n=14，根据给定的勘探数据，可以计算出各个勘探点的水合物饱和度和有效厚度。假设勘探点i的水合物饱和度为$S_i$，有效厚度为$h_i$，则该点的天然气水合物资源量为$V_i=S_i\ast h_i$。
根据这些数据，可以计算出平均资源量和总资源量：
$$\bar{V}=\frac{1}{14}\sum _{i=1}^{14}{V}_i=\frac{1}{14}(V_1+V_2+...+V_{14})$$
$$V_{total}=\sum _{i=1}^{14}{V}_i=V_1+V_2+...+V_{14}$$

给出python示例代码如下：

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

# 定义产气量因子
gas_factor = 1

# 读取勘探数据
depth_data = np.genfromtxt('depth.csv', delimiter=',')
porosity_data = np.genfromtxt('porosity.csv', delimiter=',')
saturation_data = np.genfromtxt('saturation.csv', delimiter=',')

# 计算有效厚度
thickness = []
for i in range(len(depth_data)):
    for j in range(len(depth_data[i])):
        if saturation_data[i][j] != 0:
            if i == 0 or saturation_data[i-1][j] == 0:
                thickness.append(depth_data[i][j])
            else:
                thickness.append(depth_data[i][j] - depth_data[i-1][j])

# 计算孔隙度和饱和度的概率分布
porosity_hist = np.histogram(porosity_data, bins=10)
saturation_hist = np.histogram(saturation_data, bins=10)

# 拟合概率分布曲线
def func(x, a, b, c):
    return a * np.exp(-b * x) + c

popt_porosity, pcov_porosity = curve_fit(func, porosity_hist[1][:-1], porosity_hist[0])
popt_saturation, pcov_saturation = curve_fit(func, saturation_hist[1][:-1], saturation_hist[0])

# 计算天然气水合物资源量
gas_resource = []
for i in range(len(thickness)):
    gas_resource.append(thickness[i] * porosity_data[i] * saturation_data[i] * gas_factor)

# 绘制结果
fig = plt.figure(figsize=(10, 10))

# 绘制孔隙度和饱和度的概率分布曲线
ax1 = fig.add_subplot(2, 2, 1)
ax1.plot(porosity_hist[1][:-1], porosity_hist[0], 'b-', label='porosity')
ax1.plot(porosity_hist[1][:-1], func(porosity_hist[1][:-1], *popt_porosity), 'r--', label='fit curve')
ax1.set_xlabel('porosity')
ax1.set_ylabel('probability')
ax1.set_title('Porosity Distribution')
ax1.legend()

ax2 = fig.add_subplot(2, 2, 2)
ax2.plot(saturation_hist[1][:-1], saturation_hist[0], 'b-', label='saturation')
ax2.plot(saturation_hist[1][:-1], func(saturation_hist[1][:-1], *popt_saturation), 'r--', label='fit curve')
ax2.set_xlabel('saturation')
ax2.set_ylabel('probability')
ax2.set_title('Saturation Distribution')
ax2.legend()

# 绘制天然气水合物资源量的柱状图
ax3 = fig.add_subplot(2, 1, 2)
ax3.bar(range(len(gas_resource)), gas_resource)
ax3.set_xlabel('depth')
ax3.set_ylabel('gas resource quantity')
ax3.set_title('Gas Resource Quantity Distribution')

plt.tight_layout()
plt.show()

# 输出天然气水合物资源量的估计值
print("The estimated gas resource quantity is:", sum(gas_resource))
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70

查看完整思路详见：
【腾讯文档】第九届“数维杯”大学生数学建模挑战赛全题目深度剖析建模完整过程+详细思路+代码全解析（持续更新中！）
https://docs.qq.com/doc/DSHd6UHpNZkZoWmxw