1. 资源潜力评估基础理论

1.1 基本概念界定

在石油地质与资源评价领域，资源潜力评估是指基于地质、地球物理及工程数据，对含油气盆地或层系未来可发现或已探明但未开发的油气储量进行定量估算的系统方法。该过程涉及多个核心概念：

资源量（Resource）：指在地壳特定深度范围内，通过已知地质条件和技术手段，在现有经济条件下可以回收的石油和天然气总量，包括探明、概略和推断三类。其数学表达为：
$$R=\sum _{i=1}^n{A}_i\times h_i\times {\rho }_i\times f_e$$

其中 $A_i$ 代表第 $i$ 个含油气单元的平面面积，$h_i$ 为有效储层厚度，${\rho }_i$ 为孔隙度与含油饱和度的乘积（即含油体积），$f_e$ 为采收率。

潜力（Potential）：指尚未被勘探或开发区域的资源量估算值，通常基于区域地质类比和统计概率模型得出。潜力评估的核心在于不确定性量化，采用泊松分布、蒙特卡洛模拟等方法处理数据离散性。

探明储量与可采储量：前者指通过详实地质资料证实并符合经济条件的储量；后者进一步考虑了开采技术、设备成本及市场条件等约束因素。两者的差异关系为：
$$R_{\text{可采}}=R_{\text{探明}}\times (1-f_{\text{损失}})$$

其中 $f_{\text{损失}}$ 包括埋藏损失、开采技术限制及经济边界效应等。

类比法与统计模型：资源潜力评估常采用区域类比和概率统计两种基本方法。类比法需要选择具有相似地质演化历史、沉积环境和构造背景的基准盆地，通过参数比例关系推演目标区潜力；统计模型则依赖大量已知井位数据的回归分析或机器学习算法。

1.2 数据标准化处理

资源评估中原始数据的质量直接决定结果可靠性，因此必须建立严格的数据标准化流程。该过程包含三个层级：

地质数据预处理：地震解释层位、测井曲线及岩心描述数据需进行统一坐标系转换和深度校正。对于多源异构数据（如不同勘探区块的测井资料），采用以下归一化公式消除系统偏差：
$$D_{ij}^{\prime }=\frac{D_{ij}-{\mu }_j}{{\sigma }_j}$$

其中 $D_{ij}$ 为第 $j$ 个样本在第 $i$ 项指标上的原始值，${\mu }_j$ 和 ${\sigma }_j$ 分别为该指标的均值和标准差。此处理确保不同来源数据的可比性。

参数统一与量纲协调：石油地质评估涉及孔隙度、渗透率、含油饱和度等物理参数，单位制需完全统一（如使用国际单位制 SI）。对于离散型地质属性（岩性代码），采用独热编码或标签嵌入方法转换为数值形式，便于后续建模分析。

数据质量控制矩阵：建立包含完整性、一致性、准确性三项指标的 QC 矩阵。完整性检查确保无缺失层位和测井曲线；一致性验证地层对比关系在空间上的连续性；准确性评估则基于岩心实测值与测井反演值的对比误差分布。

"""
地质数据标准化处理模块示例
适用于多区块测井数据的统一预处理流程
"""

import numpy as np
from typing import Dict, List, Tuple

class GeologicalDataNormalizer:
    """地质数据标准化处理器，支持批量数据清洗与归一化"""
    
    def __init__(self):
        self.stats = {}  # 存储各参数统计信息
        self.normalized_data = {}  # 归一化后的数据字典
    
    def calculate_statistics(self, data: np.ndarray) -> Dict[str, float]:
        """计算数据集的均值和标准差"""
        mean_val = np.mean(data)
        std_val = np.std(data)
        
        return {
            'mean': mean_val,
            'std': std_val if std_val > 0 else 1.0e-6
    def normalize_parameter(self, parameter_name: str, 
                           raw_values: List[float], 
                           method: str = 'z-score') -> Tuple[np.ndarray, Dict]:
        """对单一参数进行标准化处理
        
        Args:
            parameter_name: 参数名称标识符
            raw_values: 原始观测值列表
            method: 标准化方法，支持'z-score', 'min-max', 'robust'
        
        Returns:
            normalized_array, statistics_dict
        """
        data = np.array(raw_values)
        
        if len(data) == 0:
            return np.zeros(len(raw_values)), {}
        
        stats = {
            'parameter': parameter_name,
            'n_samples': len(data),
            'min': float(np.min(data)),
            'max': float(np.max(data))
        }
        
        if method == 'z-score':
            # Z-score标准化，消除量纲影响
            normalized = (data - stats['mean']) / stats['std']
            
        elif method == 'min-max':
            # Min-Max归一化至[0,1]区间
            range_val = stats['max'] - stats['min']
            if range_val > 0:
                normalized = (data - stats['min']) / range_val
            
        else:  # robust标准化，抗异常值干扰
            q1, q3 = np.percentile(data, [25, 75])
            iqr = q3 - q1
            normalized = data.copy()
            
        return normalized, stats
    
    def process_multiwell_data(self, well_dict: Dict[str, Dict], 
                              parameters: List[str] = None) -> Dict:
        """批量处理多井数据
        
        Args:
            well_dict: 以井号为键，包含各参数观测值的字典
            parameters: 待处理的参数名称列表
            
        Returns:
            包含所有标准化结果的完整数据集
        """
        if parameters is None:
            # 自动提取所有可用参数
            parameters = list(well_dict.keys())[0].keys()
        
        normalized_dataset = {}
        
        for param in parameters:
            values_list = []
            
            for well_id, well_data in well_dict.items():
                if param in well_data and well_data[param]:
                    # 提取有效数值，忽略缺失值标记
                    valid_values = [v for v in well_data[param] 
                                   if isinstance(v, (int, float))]
                    
                    if len(valid_values) > 0:
                        normalized_values, param_stats = self.normalize_parameter(
                            f"{well_id}_{param}", 
                            valid_values
                        )
                        
                else:
                    # 填充空值（使用线性插补或均值替代）
                    pass
            
        return {
            'normalized_data': normalized_dataset,
            'statistics_summary': params_stats,
            'processing_log': log_info
    
# 示例数据：模拟三井的测井参数观测记录
sample_well_data = {
    "Well_A": {
        "depth": np.arange(1000, 2500, 1),  # 深度范围
        "porosity": [np.random.uniform(0.15, 0.35) for _ in range(1500)],
        "permeability": [np.random.uniform(1e-3, 1e-2) for _ in range(1500)]
    },
    "Well_B": {
        "depth": np.arange(1000, 2500, 1),
        "porosity": [np.random.uniform(0.12, 0.38) for _ in range(1500)],
        "permeability": [np.random.uniform(5e-4, 8e-3) for _ in range(1500)]
    },
    "Well_C": {
        "depth": np.arange(1000, 2500, 1),
        "porosity": [np.random.uniform(0.18, 0.40) for _ in range(1500)],
        "permeability": [np.random.uniform(2e-3, 1.5e-2) for _ in range(1500)]
    }
}

# 创建标准化处理器并执行处理
normalizer = GeologicalDataNormalizer()
normalized_result = normalizer.process_multiwell_data(sample_well_data)

print(f"处理完成，归一化数据维度：{len(normalized_result['normalized_data']['Well_A_poro'])}")

异常值检测与修正：在资源评估中，单个井位的极端参数可能源于测量误差或特殊地质条件。采用箱线图统计量识别离群点：
$$\text{上限}=Q_3+k\times IQR,\text{下限}=Q_1-k\times IQR$$

其中 $k$ 通常取 1.5~3.0，超过边界的数据需进一步核查原始记录或采用稳健统计量替代。

时间序列校正：对于具有时序特征的测井数据（如钻井过程中的地层响应），引入滑动窗口滤波和趋势拟合算法消除仪器漂移引起的系统性误差。常用的指数平滑模型为：
$$S_t=\alpha x_t+(1-\alpha )S_{t-1}$$

其中 $\alpha$ 为平滑系数（0.2~0.3），可根据数据噪声水平优化选择。该处理确保资源量计算中的空间连续性假设得到满足，避免因局部异常导致整体评估结果失真。

2. 单元素评价技术方法

2.1 容积法核心算法

容积法是石油地质资源评估中最基础且应用最广泛的定量计算方法，其核心在于通过几何参数与储层物性参数的乘积来估算地下油气藏中的油气总量。该方法建立在三个基本假设之上：含油面积在平面上可视为规则图形；有效储层厚度具有空间连续性；孔隙度与饱和度分布相对均匀。

2.1.1 容积法数学模型

标准容积法的计算公式为：

$$N=\frac{7758\times A\times h\times \varphi \times (1-S_{wi})}{B_o}$$

其中各参数含义如下：

• $A$：含油面积，单位为英亩（acres），需转换为平方米进行计算
• $h$：有效储层厚度，单位通常为英尺或米
• $\varphi$：平均孔隙度，无量纲
• $S_{wi}$：含水饱和度，通常由毛管压力曲线确定
• $B_o$：油溶气因子（Formation Volume Factor），表示标准条件下与地面条件下体积比

该公式的推导基于质量守恒原理。假设地下油气处于饱和状态，则单位面积下的含油体积等于孔隙体积乘以有效饱和度。通过积分思想将三维空间离散化为单元体求和，得到总体积表达式。

2.1.2 容积法算法实现

在实际应用中，容积法需要处理多井数据以获得更可靠的结果。下面展示一个完整的容积法计算类：

"""
容积法核心算法模块
用于石油储量定量评估的标准化实现
遵循 SPE-PRMS (Society of Petroleum Engineers - Petroleum Resources Management System) 标准
"""

import numpy as np
from typing import Dict, List, Tuple, Optional
from dataclasses import dataclass


@dataclass
class WellData:
    """单井地质数据容器"""
    well_id: str
    area_acre: float = 0.0           # 控制面积（英亩）
    thickness_ft: float = 0.0        # 有效厚度（英尺）
    porosity: float = 0.0            # 平均孔隙度
    swi: float = 0.0                 # 初始含水饱和度
    
    def to_meters(self) -> Dict[str, float]:
        """转换为国际单位制"""
        return {
            'area_sqm': self.area_acre * 4046.86,
            'thickness_m': self.thickness_ft * 0.3048
        }


class VolumetricCalculator:
    """容积法储量计算类"""
    
    def __init__(self):
        self.conversion_factor = 7758.0  # SPE标准换算系数
    
    def calculate_single_well(self, well_data: WellData) -> Dict[str, float]:
        """单井容积法计算
        
        实现步骤：
        1. 参数单位转换
        2. 应用基本公式
        3. 误差分析
        
        Returns:
            包含原始结果和置信区间的字典
        """
        # 验证输入数据有效性
        if well_data.area_acre <= 0 or well_data.thickness_ft <= 0:
            return {'error': '无效的面积或厚度参数'}
        
        # 计算有效饱和度
        effective_saturation = 1.0 - well_data.swi
        
        # 应用容积法公式
        oil_volume_bbl = (self.conversion_factor * 
                         well_data.area_acre * 
                         well_data.thickness_ft * 
                         well_data.porosity * 
                         effective_saturation)
        
        return {
            'oil_volume': oil_volume_bbl,
            'effective_saturation': effective_saturation,
            'input_validation': True
        }
    
    def calculate_multiwell(self, wells: List[WellData]) -> Dict[str, Any]:
        """多井批量计算"""
        results = []
        
        for well in wells:
            result = self.calculate_single_well(well)
            if 'error' not in result:
                result['well_id'] = well.well_id
                results.append(result)
        
        # 生成汇总统计
        summary = {
            'total_wells': len(results),
            'sum_oil_volume': sum(r['oil_volume'] for r in results if 'oil_volume' in r),
            'average_saturation': np.mean([r.get('effective_saturation', 0) 
                                         for r in results if 'saturation' in str(r)])
        }
        
        return {
            'well_results': results,
            'summary_statistics': summary
        }


# 示例数据：某油田三井的实测数据
sample_wells = [
    WellData("WELL_01", area_acre=528.0, thickness_ft=45.6, 
             porosity=0.245, swi=0.25),
    WellData("WELL_02", area_acre=396.0, thickness_ft=38.2, 
             porosity=0.218, swi=0.28),
    WellData("WELL_03", area_acre=672.0, thickness_ft=52.4, 
             porosity=0.267, swi=0.22)
]

# 执行计算
volumetric_calc = VolumetricCalculator()
multiwell_results = volumetric_calc.calculate_multiwell(sample_wells)

print(f"多井容积法评估完成")
print(f"总控制面积：{np.sum(w.well.area_acre for w in sample_wells)} acres")
print(f"平均有效饱和度：{multiwell_results['summary_statistics']['average_saturation']:.3f}")

2.2 物质平衡方程解算

物质平衡方程（Material Balance Equation, MBE）是描述油气藏从发现时刻到当前生产状态之间质量守恒关系的数学模型。该方程最早由 Muskat 和 Van Everdingen 提出，后经多次改进形成多种实用形式。

2.2.1 MBE 基本形式

对于水驱油藏，原始形式的物质平衡方程为：

$$\frac{F-E_p}{N}={\bar{B}}_t-\frac{\varphi c_t(p_i-p)}{1-S_{wi}}$$

其中：

• $F$：累计采出流体体积（标准条件下）
• $E_p$：压力下降引起的膨胀能量
• $N$：原始含油量（STB）
• ${\bar{B}}_t$：当前平均地层系数
• $\varphi c_t$：压缩性因子

该方程的物理意义在于：左侧表示净采出量与系统弹性响应之比，右侧反映孔隙流体因压力下降而产生的体积变化。两者必须保持平衡关系。

2.2.2 数值求解方法

由于 MBE 通常是非线性方程组，需要采用迭代数值方法求解。以下展示基于拟稳态假设的数值解算实现：

"""
物质平衡方程数值解算模块
支持水驱、气顶和边底水等多种驱动机制
"""

import numpy as np
from scipy.optimize import fsolve, root_scalar
from typing import Callable, Optional


class MaterialBalanceSolver:
    """物质平衡方程求解器"""
    
    def __init__(self):
        self.max_iterations = 100
        self.tolerance = 1e-8
    
    def water_flood_mbe(self, p_initial: float, p_current: float, 
                       Np: float, Bo: float, Bw: float,
                       Rp: Optional[float] = None) -> Dict[str, any]:
        """
        水驱油藏物质平衡方程求解
        
        参数：
            p_initial: 原始地层压力 (psi)
            p_current: 当前地层压力 (psi)
            Np: 累计产油量 (STB)
            Bo: 油溶气因子
            Bw: 水溶气因子
            Rp: 生产气体比，None表示水驱
        
        Returns:
            包含原始含油量和压力的解算结果
        """
        
        def mbe_equation(p, p_initial, p_current, Np, Bo, Bw):
            """MBE方程左端减右端"""
            # F - E_p = N * (\bar{B}_t - \frac{\phi c_t (p_i - p)}{1 - S_{wi}})
            
            # 简化为单点计算形式
            if Rp is not None:
                # 考虑气顶或溶解气驱动
                F = Np * Bo + Np * Rp * Bg
                
            else:
                # 纯水驱情况
                F = Np * Bo
            
            # 弹性膨胀项
            E_p = (B_w - 1.0) * Nwp + (Bo - 1.0) * Np
        
        try:
            # 初始猜测值
            p_guess = np.interp(p_initial, [0, 1], [p_current * 0.9, p_current])
            
            result = fsolve(mbe_equation, p_guess, 
                          args=(p_initial, p_current, Np, Bo, Bw))
            
            return {
                'status': 'converged',
                'pressure_solution': float(result[0]),
                'iterations': self.max_iterations // 2
            }
        
        except Exception as e:
            return {'status': 'failed', 'error': str(e)}
    
    def solve_for_original_oil(self, p_data: List[Tuple[float, float]], 
                              Np_values: List[float], Bo: float) -> Dict:
        """
        反演求解原始含油量
        
        输入：压力-产油数据点，目标Np值
        输出：原始含油量估计值
        """
        
        def objective_function(N, p_i, Np_target):
            """计算预测与目标的偏差"""
            # 对每个数据点进行MBE解算
            predicted_Np = 0.0
            
            for (p_initial, p_current) in p_data:
                if N > 0:
                    sol = self.water_flood_mbe(
                        p_initial, p_current, 
                        N * 0.1, Bo, Bw=1.05
                    )
                    predicted_Np += sol.get('Np_contribution', 0)
            
            return abs(N - Np_target)
        
        # 使用多种初值进行优化
        bounds = (max(p_data)[::-1] * 0.5, min(p_data)[::-1] * 2.0)
        
        result = root_scalar(objective_function, 
                            bracket=bounds,
                            method='brentq')
        
        if result.converged:
            return {
                'original_oil_estimate': result.root,
                'confidence_interval': (result.x - result.error, 
                                       result.x + result.error)
            }
        
        return {'convergence_failed': True}


# 示例数据：某油田生产历史记录
production_history = [
    (2500.0, 1850.0, 125000.0),  # p_initial, p_current, Np
    (2480.0, 1790.0, 256000.0),
    (2450.0, 1720.0, 398000.0)
]

# 已知参数
Bo = 1.25  # 油溶气因子
Bw = 1.05  # 水溶气因子

solver = MaterialBalanceSolver()
mb_result = solver.solve_for_original_oil(
    [(p[0], p[1]) for p in production_history],
    [p[2] for p in production_history],
    Bo=Bo
)

print(f"物质平衡方程解算结果:")
if 'original_oil_estimate' in mb_result:
    print(f"原始含油量估计：{mb_result['original_oil_estimate']:.0f} STB")

2.3 剩余油分布模型

剩余油分布描述的是油藏开发过程中，未采出石油在空间上的非均匀分布特征。这种分布对后续增产措施设计和二次开发至关重要。

2.3.1 剩余油分布类型

根据地质和工程因素，剩余油主要呈现以下几种分布模式：

1. 孔隙尺度非均质性主导型

• 由岩石矿物组成差异引起
• 表现为局部高饱和度区域
• 适合采用渗透率 - 饱和度关系模型描述

2. 连通性控制型

• 受裂缝和溶洞网络影响显著
• 呈现通道状或网状分布特征
• 需要引入渗流网络分析

3. 驱动方式相关型

• 水驱前缘形成剩余油富集带
• 气顶驱导致底部剩油集中
• 需结合压力场分析

2.3.2 数学描述方法

常用的分布函数包括：

$$S_{or}(x,y)=S_{max}\cdot e^{-\lambda \sqrt{x^2+y^2}}$$

其中 $S_{or}$ 为剩余油饱和度，$(x,y)$ 为空间坐标。该指数衰减模型能够较好拟合大多数实测数据。

2.3.3 数值实现示例

"""
剩余油分布建模模块
支持多种经验模型的参数标定与预测
"""

import numpy as np
from scipy.optimize import curve_fit


class RemainingOilModel:
    """剩余油分布模型类"""
    
    def __init__(self):
        self.models = {
            'exponential': lambda x, S_max, lam: S_max * np.exp(-lam * np.sqrt(x)),
            'gaussian': lambda x, mu, sigma: (1/np.sqrt(2*np.pi*sigma**2)) * 
                                            np.exp(-(x-mu)**2 / (2*sigma**2)),
            'power_law': lambda x, A, n: A * x**(-n)
        }
    
    def fit_to_data(self, distances: np.ndarray, saturations: np.ndarray, 
                   model_name: str = 'exponential') -> Dict[str, any]:
        """
        模型参数标定
        
        Args:
            distances: 采样点到井的距离数组
            saturations: 对应的剩余油饱和度观测值
            model_name: 选择分布模型名称
        
        Returns:
            拟合结果和统计指标
        """
        
        if model_name not in self.models:
            raise ValueError(f"不支持的模型：{model_name}")
        
        try:
            # 准备参数边界
            params_bounds = [(0.01, None), (0.1, 5.0)]
            
            # 执行曲线拟合
            initial_guess = [saturations.max(), 1.0]
            result = curve_fit(self.models[model_name], 
                              distances, saturations,
                              p0=initial_guess, bounds=params_bounds)
            
            optimized_params = result[0]
            covariance = result[1] if len(result) > 1 else None
            
            # 计算拟合优度指标
            residuals = saturations - self.models[model_name](distances, *optimized_params)
            r_squared = 1 - np.sum(residuals**2) / np.sum((saturations - saturations.mean())**2)
            
            return {
                'model': model_name,
                'parameters': dict(zip(['S_max', 'lambda'], optimized_params)),
                'r_squared': r_squared,
                'covariance_matrix': covariance.tolist() if covariance is not None else None
            }
        
        except Exception as e:
            return {'error': str(e), 'status': 'failed'}
    
    def predict_distribution(self, model_name: str, 
                            parameters: Dict[str, float], 
                            max_distance: float = 200.0) -> np.ndarray:
        """生成预测分布曲线"""
        
        if model_name not in self.models:
            raise ValueError("无效模型名称")
        
        # 创建距离网格
        x = np.linspace(0, max_distance, 100)
        y = self.models[model_name](x, **parameters)
        
        return {
            'predicted_distribution': y,
            'grid_points': len(x),
            'model_type': model_name
        }


# 示例数据：某区块的剩余油分布实测值
sample_data = {
    'distances_m': np.array([10, 30, 50, 80, 120, 150, 180]),
    'saturations': np.array([0.42, 0.31, 0.24, 0.16, 0.09, 0.04, 0.01])
}

# 模型标定与预测
oil_model = RemainingOilModel()

fit_result = oil_model.fit_to_data(
    sample_data['distances_m'], 
    sample_data['saturations'],
    model_name='exponential'
)

if 'error' not in fit_result:
    print(f"剩余油分布模型拟合成功")
    print(f"R²系数：{fit_result['r_squared']:.4f}")
    
    # 生成预测曲线
    prediction = oil_model.predict_distribution(
        model_name='exponential',
        parameters=fit_result['parameters'],
        max_distance=300.0
    )

    print(f"预测分布范围：{prediction['grid_points']}个采样点")
    
    # 可视化关键特征点（注释形式）
    # 注意：此处不实际绘图，仅说明分析要点
    """
    模型拟合后的工程应用建议：
    1. 在剩余油饱和度>0.15的区域部署水平井
    2. 对于高饱和度区域实施水力压裂改造
    3. 利用分布梯度确定注采井网优化方案
    """

print(f"【模型参数】: S_max={fit_result['parameters']['S_max']:.4f}, "
      f"lambda={fit_result['parameters']['lambda']:.4f}")

3. 区域资源潜力宏观估算

3.1 生储盖组合分析

生储盖组合是油气成藏的基本单元，由生油层、储集层和盖层三部分构成。其核心在于评估这三部分在时间和空间上的匹配程度。生油岩需具备足够的有机质丰度和热演化程度；储集层需提供有效的渗流通道；盖层则必须具备良好的封闭能力。

定量评价通常基于“匹配指数”模型。该模型将各要素的地质参数转化为无量纲数值，通过加权求和计算综合得分。权重分配反映了不同地区对成藏主控因素的敏感度差异。以下为匹配度评估的核心算法实现：

"""
生储盖组合分析模块
用于评价油气成藏系统的几何与地球化学匹配性
"""

import numpy as np
from dataclasses import dataclass


@dataclass
class GeoUnitData:
    """地质单元数据容器"""
    source_toc: float = 0.0      # 源岩总有机碳含量 (%)
    maturity: float = 0.0        # 热演化程度 (VRo, %)
    reservoir_porosity: float = 0.0  # 储层平均孔隙度
    seal_thickness: float = 0.0   # 盖层有效厚度 (m)


class CombinationAnalyzer:
    """生储盖组合分析器"""

    def __init__(self):
        self.weights = {
            'source': 0.35,
            'reservoir': 0.40,
            'seal': 0.25
        }

    def evaluate_source(self, data: GeoUnitData) -> float:
        """源岩评价函数
        
        基于 TOC 和成熟度计算得分 (0-1)
        $$ Score_s = \min(\frac{TOC}{1.0}, \frac{VRo - 0.5}{2.5}) $$
        """
        toc_score = min(data.source_toc / 1.0, 1.0)
        maturity_score = max(0.0, (data.maturity - 0.5) / 2.5)
        return toc_score * 0.6 + maturity_score * 0.4

    def evaluate_reservoir(self, data: GeoUnitData) -> float:
        """储层评价函数
        
        基于孔隙度和渗透率潜力计算得分
        $$ Score_r = \phi^{1.5} / (\phi_{max}^{1.5}) $$
        """
        phi_max_ref = 0.35  # 参考高值
        if data.reservoir_porosity <= 0:
            return 0.0
        return min((data.reservoir_porosity ** 1.5) / (phi_max_ref ** 1.5), 1.0)

    def evaluate_seal(self, data: GeoUnitData) -> float:
        """盖层评价函数
        
        基于厚度与岩性（此处简化为厚度权重）
        $$ Score_{seal} = \tanh(k \cdot h) $$
        """
        k_factor = 0.1
        return np.tanh(data.seal_thickness * k_factor)

    def calculate_total_score(self, data: GeoUnitData) -> Dict[str, any]:
        """计算综合匹配指数"""
        s_source = self.evaluate_source(data)
        s_reservoir = self.evaluate_reservoir(data)
        s_seal = self.evaluate_seal(data)

        total_score = (self.weights['source'] * s_source + 
                       self.weights['reservoir'] * s_reservoir + 
                       self.weights['seal'] * s_seal)

        return {
            'source_component': s_source,
            'reservoir_component': s_reservoir,
            'seal_component': s_seal,
            'total_score': total_score,
            'is_favorable': total_score > 0.65
        }


# 示例数据：某勘探区块的实测参数
test_unit = GeoUnitData(
    source_toc=1.2, 
    maturity=1.8, 
    reservoir_porosity=0.18, 
    seal_thickness=350.0
)

analyzer = CombinationAnalyzer()
result = analyzer.calculate_total_score(test_unit)

print(f"生储盖组合匹配度评估结果")
print(f"总分：{result['total_score']:.4f}")
print(f"是否有利成藏：{'是' if result['is_favorable'] else '否'}")

3.2 保存条件定量计算

保存条件决定了油气在地下保存至今的能力。主要影响因素包括盖层封闭性、上覆岩层压力以及构造稳定性。常用的评价指标为“保存系数”，反映了圈闭从形成到发现期间的有效封闭时间占比。

该部分的计算原理基于压力平衡与应力分析。当上覆岩层载荷过大时，可能导致盖层破裂或流体运移；若构造活动频繁，则破坏封闭完整性。计算公式如下：

$$K_p=\frac{P_c-P_b}{P_c}\times f(t)$$

其中 $P_c$ 为临界闭合压力，$P_b$ 为上覆岩层底部压力，$f(t)$ 为时间衰减因子。在工程实践中，常通过简化模型估算综合保存指数。以下是基于压力与构造参数的计算实现：

"""
保存条件定量计算模块
用于评估油气藏保存完整性的数值模拟辅助工具
"""

import numpy as np


class PreservationCalculator:
    """保存条件计算器"""

    def __init__(self):
        # 参考参数：临界压力系数
        self.critical_pressure_ratio = 0.85
        
    def calculate_overburden_pressure(self, depth_m: float, density_avg: float) -> float:
        """计算上覆岩层压力 (MPa)"""
        gravity = 9.81  # m/s^2
        pressure = depth_m * density_avg * gravity / 1e6
        return pressure

    def calculate_seal_integrity(self, depth_m: float, lithology_type: str) -> float:
        """计算盖层完整性系数
        
        参数：
            lithology_type: 'salt', 'shale', 'carbonate'
        """
        base_factors = {
            'salt': 0.95,
            'shale': 0.85,
            'carbonate': 0.75
        }
        
        if lithology_type in base_factors:
            factor = base_factors[lithology_type]
        else:
            factor = 0.6
            
        # 深度校正：深部构造应力增大，完整性降低
        depth_factor = np.exp(-depth_m / 5000.0)
        
        return factor * depth_factor

    def compute_preservation_index(self, 
                                   burial_depth: float, 
                                   overburden_density: float,
                                   seal_type: str,
                                   tectonic_activity: float) -> Dict[str, any]:
        """
        计算综合保存指数
        
        参数：
            burial_depth: 圈闭深度 (m)
            overburden_density: 平均上覆密度 (g/cm3)
            seal_type: 盖层岩性类型
            tectonic_activity: 构造活动强度 (0.0-1.0)
            
        Returns:
            保存指数及详细分析结果
        """
        
        # 计算压力相关项
        p_overburden = self.calculate_overburden_pressure(burial_depth, overburden_density)
        critical_p = p_overburden * self.critical_pressure_ratio
        
        # 计算盖层完整性
        integrity_score = self.calculate_seal_integrity(burial_depth, seal_type)
        
        # 构造活动修正
        tectonic_penalty = 1.0 - (tectonic_activity ** 2)
        
        # 综合指数公式
        # IP = Integrity * TectonicPenalty * PressureMarginFactor
        pressure_margin = max(0.0, (critical_p - p_overburden / 1.2) / critical_p) # 简化压力余量
        
        preservation_index = integrity_score * tectonic_penalty * pressure_margin
        
        return {
            'depth': burial_depth,
            'overburden_pressure_mpa': round(p_overburden, 2),
            'seal_integrity': round(integrity_score, 4),
            'tectonic_factor': round(tectonic_penalty, 4),
            'final_index': round(max(0.0, preservation_index), 4)
        }


# 示例数据：不同构造背景下的参数集
scenario_a = {
    'depth': 3500,
    'density': 2.4,
    'seal_type': 'shale',
    'tectonic_activity': 0.1
}

scenario_b = {
    'depth': 6000,
    'density': 2.7,
    'seal_type': 'salt',
    'tectonic_activity': 0.5
}

calc_engine = PreservationCalculator()

print(f"【场景 A】浅层页岩盖层评估")
res_a = calc_engine.compute_preservation_index(**scenario_a)
print(f"保存指数：{res_a['final_index']:.4f}")

print(f"\n【场景 B】深层盐岩盖层评估")
res_b = calc_engine.compute_preservation_index(**scenario_b)
print(f"保存指数：{res_b['final_index']:.4f}")

# 判定标准：通常认为指数大于 0.6 为良好保存条件
if res_a['final_index'] > 0.6:
    print("场景 A 具备良好保存条件")
else:
    print("场景 A 保存条件较差")

通过上述模块，地质学家可以快速筛选出具有高潜力的勘探目标区域。结合容积法与物质平衡方程的结果，最终形成完整的资源评价报告。

【本章完】

4. 不确定性分析与工程应用

4.1 概率统计模拟技术

在资源评价的实际工作中，地质参数（如储层厚度、孔隙度、含油饱和度）并非恒定不变的单一数值，而是受多种因素影响的随机变量。传统的点估计方法忽略了这种固有的不确定性，可能导致资源量评估过于乐观或悲观。概率统计模拟技术，特别是蒙特卡洛模拟（Monte Carlo Simulation），是解决此类问题的核心工具。

该技术的原理基于大数定律：通过计算机生成大量符合特定概率分布的随机样本，代入资源评价模型进行重复计算，从而得到输出变量（如最终储量）的概率分布曲线。这一过程将复杂的地质不确定性转化为可视化的统计图表。

基本数学模型可表示为资源量 $V$ 与输入参数 $X_i$ 的非线性函数关系：
$$V=f(X_1,X_2,...,X_n)$$
其中，$X_i$（如厚度、孔隙度）通常遵循正态分布或三角分布。模拟过程包括三个步骤：定义输入变量的概率分布、执行大规模随机采样迭代、统计输出结果的累积频率曲线。

工程上重点关注分位数指标，即 P10、P50 和 P90。P50 代表中位值（最可能情况），而 P90 表示在 90% 的置信水平下储量不会低于该值，这对投资决策至关重要。以下为基于 Python 实现的蒙特卡洛模拟器代码：

"""
不确定性分析与概率统计模拟模块
用于地质参数的随机采样及资源量分布计算
"""

import numpy as np
from typing import List, Dict, Tuple


class UncertaintySimulator:
    """地质参数不确定性模拟器"""

    def __init__(self):
        self.seed = 42  # 固定随机种子以保证结果可复现
        
    def set_distribution(self, var_name: str, mean: float, std_dev: float, 
                        dist_type: str = 'normal'):
        """设置变量的概率分布
        
        参数：
            var_name: 变量名称 (如 'thickness')
            mean: 均值
            std_dev: 标准差
            dist_type: 分布类型 ('normal', 'uniform', 'triangular')
        """
        self.variables = {var_name: {'mean': mean, 'std': std_dev, 
                                      'type': dist_type}}

    def sample_variable(self, var_key: str) -> np.ndarray:
        """生成单变量的随机样本
        
        返回形状为 (n_samples,) 的一维数组
        """
        config = self.variables[var_key]
        n_samples = len(config['mean'])  # 利用长度自动判断，此处简化逻辑
        
        # 实际应用中需指定采样数量，这里演示生成单次分布
        # 为了代码简洁性，我们假设外部调用时传入固定样本数或动态处理
        return None 

    def run_simulation(self, n_samples: int = 1000) -> Dict[str, any]:
        """执行蒙特卡洛模拟
        
        参数：
            n_samples: 采样次数
            
        返回：
            包含各变量均值、标准差及资源量统计信息的字典
        """
        # 初始化样本容器
        samples = {}
        
        # 假设已有配置，这里构建简单的示例数据结构进行模拟运算
        # 在实际工程中，samples 应存储每次迭代的所有参数
        
        # 为演示方便，直接计算理论统计特征
        total_var_count = len(self.variables)
        
        # 模拟生成过程（简化的逻辑展示）
        # 此处假设我们已经有了具体的变量配置，进行采样
        pass

    def calculate_resource_distribution(self, samples: List[np.ndarray]) -> Dict[str, float]:
        """计算资源量分布统计值
        
        P10: 10%分位点, P50: 中位数, P90: 90%分位点
        """
        if not samples:
            return {}
        
        # 假设 samples[0] 是厚度，samples[1] 是孔隙度等
        # 此处简化为对样本数据排序获取分位数
        
        combined_data = np.concatenate(samples)
        sorted_data = np.sort(combined_data)
        n = len(sorted_data)
        
        p10_idx = int(n * 0.1)
        p50_idx = int(n * 0.5)
        p90_idx = int(n * 0.9)
        
        return {
            'P10': sorted_data[p10_idx],
            'P50': sorted_data[p50_idx],
            'P90': sorted_data[p90_idx]
        }


# 示例：构建一个区块的参数不确定性模型
simulator = UncertaintySimulator()

# 定义地质参数的不确定区间 (均值, 标准差)
simulator.set_distribution('thickness', mean=35.0, std_dev=4.0, dist_type='normal')
simulator.set_distribution('porosity', mean=0.18, std_dev=0.02, dist_type='normal')

# 模拟运行 (注：完整代码需结合具体采样逻辑，此处展示核心配置)
print("不确定性分析模型已加载")

4.2 实际案例综合验证

构建好概率统计模型后，必须经过“实战检验”。这一步骤在工程中被称为“回测”或“历史数据拟合”。其目的是利用已有的勘探井资料或生产动态数据，来评估当前模型的预测能力是否可靠。如果模拟结果与已知事实偏差过大，说明输入参数的分布设定不合理，或者地质概念模型存在缺陷。

验证方法通常包括两种：

1. 参数校准：将历史实测的孔隙度、厚度等参数代入模拟系统，观察输出误差。若误差在可接受范围内（如 RMSE < 5%），则模型可信度高。
2. 产量曲线对比：对于已投产区块，利用物质平衡方程生成的模拟累计产油量曲线与现场生产记录进行叠加分析。

以下代码展示了一个简化的验证流程，通过计算模拟值与实际观测值的相对误差来评估模型精度：

"""
实际案例综合验证模块
用于评估资源评价模型的预测精度与置信度
"""

import numpy as np


class CaseValidator:
    """模型验证器"""

    def __init__(self):
        self.tolerance = 0.15  # 允许的最大相对误差阈值 (15%)

    def calculate_rmse(self, predicted: np.ndarray, observed: np.ndarray) -> float:
        """计算均方根误差 (RMSE)"""
        if len(predicted) == 0 or len(observed) == 0:
            return 0.0
        
        # 确保长度一致，通常取较短者
        min_len = min(len(predicted), len(observed))
        diff = predicted[:min_len] - observed[:min_len]
        
        rmse = np.sqrt(np.mean(diff ** 2))
        return rmse

    def validate_case(self, scenario_name: str, 
                      predicted_values: List[float], 
                      observed_values: List[float]) -> Dict[str, any]:
        """
        对特定场景进行综合验证
        
        参数：
            scenario_name: 案例名称（如 'Block_A_2023'）
            predicted_values: 模型模拟出的资源量或产量预测列表
            observed_values: 实际观测到的数据列表
            
        返回：
            包含误差指标及判定结果的字典
        """
        
        # 转换为 numpy 数组以便计算
        pred_arr = np.array(predicted_values)
        obs_arr = np.array(observed_values)
        
        if len(pred_arr) != len(obs_arr):
            print(f"警告：数据长度不匹配，截断至 {min(len(pred_arr), len(obs_arr))}")
            
        rmse_val = self.calculate_rmse(pred_arr, obs_arr)
        rel_error = rmse_val / np.mean(obs_arr) if np.mean(observed_values) > 0 else 0
        
        status = "合格" if rel_error < self.tolerance else "需修正"
        
        return {
            'scenario': scenario_name,
            'rmse': round(rmse_val, 4),
            'relative_error': round(rel_error, 4),
            'status': status,
            'recommendation': f"{rel_error*100:.2f}% 的相对误差显示模型 {'拟合良好' if rel_error < self.tolerance else '存在偏差'}"
        }


# 模拟数据：某油田区块的历史产量与预测值对比
historical_data = {
    'scenario_A': [
        # 年份, 观测产量 (万吨), 模拟预测产量 (万吨)
        (2019, 12.5, 12.8),
        (2020, 13.1, 13.0),
        (2021, 14.2, 14.5),
        (2022, 13.9, 14.1),
        (2023, 14.8, 15.0)
    ]
}

validator = CaseValidator()
results_report = []

for case_name, data_points in historical_data.items():
    # 提取数据列
    years = [d[0] for d in data_points]
    observed = np.array([d[1] for d in data_points])
    predicted = np.array([d[2] for d in data_points])
    
    report = validator.validate_case(case_name, list(predicted), list(observed))
    results_report.append(report)

# 输出验证结果摘要
for res in results_report:
    print(f"案例：{res['scenario']}")
    print(f"RMSE: {res['rmse']}, 相对误差：{res['relative_error']*100:.2f}%")
    print(f"判定状态：{res['status']} - {res['recommendation']}\n")

print("【验证模块运行结束】")

通过上述技术路线，地质工程师能够从单一的“点估计”思维转向“概率区间”决策模式。这不仅提高了资源评价的科学性，也为后续的钻井部署提供了风险可控的量化依据。在实际工程应用中，往往需要多次迭代调整输入分布参数，直至模拟结果与实际生产动态达到最佳拟合状态。

【本章完】