基于天然气水合物勘探的多功能CPTU数据解释与聚类算法改进【附代码】
✨ 长期致力于天然气水合物、多功能静力触探、聚类、分层、浅层气研究工作，擅长数据搜集与处理、建模仿真、程序编写、仿真设计。
✅ 专业定制毕设、代码
✅如需沟通交流，点击《获取方式》
（1）孔穴扩张理论与有限元耦合的锥尖阻力解释模型：
针对海洋粘土中CPTU锥尖阻力，提出孔穴扩张理论与大应变有限元耦合解释方法。建立轴对称有限元模型，土体采用修正剑桥模型，参数包括刚度指数Ir=200、应变软化系数δ=0.8。通过改变探头粗糙度系数（0.1~0.9），获得锥形因子Nc与刚度指数的关系式：Nc = 4.2 + 0.9*ln(Ir) - 0.3*α，其中α为粗糙度。利用该关系式，从实测锥尖阻力反演不排水抗剪强度，误差小于8%。同时，模拟了土体塑性区扩展，当锥尖阻力达到5MPa时，塑性区半径约为探头半径的5.2倍，与孔穴扩张理论预测一致。
（2）层次聚类与Robertson分类图结合的自动分层算法：
提出一种基于CPTU数据的层次聚类分层流程，命名为CPTU-HClust。首先对原始数据进行标准化（z-score），然后采用Ward距离作为聚类准则。定义最大聚类个数Ncmax=15，通过轮廓系数评估确定最优聚类个数。在聚类结果中，结合Robertson和Campanella分类图，将每一类映射到土类代码（如砂土、粘土、粉土）。算法用Matlab实现并封装为函数。在南海天然气水合物勘探区CPTU实测数据处理中，聚类算法成功识别出6个主层和3个过渡夹层，与钻孔岩心对比，分层深度误差小于0.3米。特别是识别出水合物储层上覆的致密层，其锥尖阻力值在8~12MPa之间，聚类单独成类。
（3）多功能参数联合识别水合物与浅层气储层：
利用CPTU的多功能参数（电阻率Rt、声波速度AC、归一化土类指数ISBT）进行聚类初始参数选择。针对水合物储层，选用电阻率、声波时差和自然伽马作为初始聚类参数，聚类结果清晰区分出水合物层（电阻率>50Ω·m，声波时差<200μs/m）和泥岩层。针对浅层气储层，则采用ISBT、Rt和AC三个参数，含气层表现为低电阻率（<5Ω·m）和低声波速度（<1500m/s）且ISBT指数大于3.0。在印尼某浅层气勘探区应用，聚类算法成功划分出含气层厚度2.5米，与后期钻探验证一致。该方法减少了人工解释的主观性，将分层效率提高5倍。
import numpy as np from scipy.cluster.hierarchy import fcluster, linkage, dendrogram from sklearn.preprocessing import StandardScaler def cptu_hclust(data, max_clusters=15): # 数据标准化 scaler = StandardScaler() data_norm = scaler.fit_transform(data) # Ward 链接 Z = linkage(data_norm, method='ward') # 轮廓系数选择最优聚类数 from sklearn.metrics import silhouette_score best_n = 2 best_score = -1 for n in range(2, max_clusters+1): labels = fcluster(Z, n, criterion='maxclust') if len(set(labels)) > 1: score = silhouette_score(data_norm, labels) if score > best_score: best_score = score best_n = n final_labels = fcluster(Z, best_n, criterion='maxclust') return final_labels, best_n def cone_factor(Ir, roughness=0.5): # 锥形因子 Nc 经验公式 Nc = 4.2 + 0.9 * np.log(Ir) - 0.3 * roughness return Nc def classify_soil(qc, fs, u2): # 简化的 Robertson 分类 (基于 qc, fs) from math import log10 # 归一化 pa = 101.3 # kPa Qtn = (qc - u2) / pa Fr = (fs / (qc - u2)) * 100 # 典型分类逻辑 if Qtn > 100 and Fr < 1: return '砂土' elif Qtn < 50 and Fr > 5: return '粘土' else: return '粉土' def identify_gas_hydrate_cluster(data_multi): # 假设 data_multi 包含 Rt, AC, ISBT 三列 # 特征加权: 电阻率权重高 weights = np.array([2.0, 1.5, 1.0]) weighted_data = data_multi * weights labels, _ = cptu_hclust(weighted_data) # 分析每类的中心 unique_labels = np.unique(labels) class_centers = [] for lab in unique_labels: center = np.mean(data_multi[labels==lab], axis=0) class_centers.append(center) # 识别含气层: 低 Rt & 低 AC hydrate_cluster = None for i, center in enumerate(class_centers): if center[0] > 50 and center[1] < 200: # 水合物特征 hydrate_cluster = unique_labels[i] break return labels, hydrate_cluster if __name__ == '__main__': # 模拟 CPTU 数据: 深度, qc, fs, u2 depth = np.linspace(0, 20, 200) qc = 5 + 2*np.sin(depth) + np.random.randn(200)*0.5 fs = 0.1 * qc + np.random.randn(200)*0.02 u2 = 100 + 5*depth + np.random.randn(200)*10 data = np.column_stack([qc, fs, u2]) labels, best_n = cptu_hclust(data) print('最优聚类数:', best_n) # 锥形因子示例 Nc_val = cone_factor(Ir=200, roughness=0.5) print('锥形因子 Nc = {:.2f}'.format(Nc_val)) # 分类示例 soil = classify_soil(qc[0], fs[0], u2[0]) print('第一层土类:', soil)