基于物理特征融合与机器学习的多井协同钻井速率实时预测与优化(

开始│├─ 数据加载与预处理│   ├─ 加载多井钻井数据│   ├─ 选择关键特征列│   ├─ 合并多井数据集│   ├─ 处理缺失值与异常值│   └─ 创建时间戳特征│├─ 特征工程│   ├─ 数据变换（Box-Cox处理偏度）│   └─ 计算物理特征（平均比能MSE）│├─ 数据分析│   ├─ 分布可视化│   ├─ 相关性分析│   └─ 特征选择│├─ 模型构建│   ├─ 数据分割（训练集/测试集）│   ├─ 特征归一化│   ├─ XGBoost模型训练│   └─ 支持向量机模型训练（对比）│├─ 模型评估│   ├─ R^2分数计算│   ├─ 误差指标分析（MAE/MSE/RMSE）│   ├─ 预测结果可视化│   └─ 特征重要性分析│└─ 模型应用    ├─ 模型保存与加载    └─ 新数据预测结束

# 导入必要的库import numpy as np  # 数值计算库import pandas as pd  # 数据处理库import matplotlib.pyplot as plt  # 数据可视化库import seaborn as sns  # 高级数据可视化库import warnings  # 警告处理warnings.filterwarnings('ignore')  # 忽略所有警告
# 加载MIP-3HA井的钻井数据MIP3A = pd.read_csv(r'datasets\MIP-3H Top.csv')  # 读取CSV文件
# 数据预处理函数：选择关键特征列def slice_columns(data):    # 定义需要保留的关键特征列    columns = ['Hole Depth', 'Rate Of Penetration', 'Bit Depth', 'Hook Load',                'Standpipe Pressure', 'Rotary RPM', 'Rotary Torque', 'Weight on Bit',                'Block Height', 'On Bottom Hours', 'Circulating Hours', 'Tool Face',                'Inclination', 'Azimuth', 'Differential Pressure', 'YYYY/MM/DD', 'HH:MM:SS']    df = data[columns]  # 选择指定列    return df
# 多井数据合并函数def combine_and_transform():    # 定义要合并的井数据集    wells = [MIP3A, MIP3B, MIP5A, MIP5B, MIPSW]      well_names = ['MIP3A', 'MIP3B', 'MIP5A', 'MIP5B', 'MIPSW']  # 井标识    combined_df = pd.DataFrame()  # 创建空DataFrame
    # 遍历每口井的数据    for well, well_name in zip(wells, well_names):        df = slice_columns(well)  # 选择关键特征        df['Well ID'] = well_name  # 添加井标识列        combined_df = pd.concat([combined_df, df], ignore_index=True)  # 合并数据
    # 创建时间戳列：合并日期和时间    combined_df['TimeStamp'] = pd.to_datetime(combined_df['YYYY/MM/DD'] + ' ' + combined_df['HH:MM:SS'])    return combined_df
# 调用函数合并多井数据df = combine_and_transform()
# 数据清洗：替换无效值为NaNdf = df.replace(-999.25, np.NaN)
# 删除工具面列（数据无变化）和空值行df = df.drop(columns='Tool Face')df = df.dropna()
# 异常值处理函数def process_outliers(df):    # 转速过滤：保留0-100 RPM之间的合理值    df = df[(df['Rotary RPM'] > 0) & (df['Rotary RPM'] < 100)]
    # 钻压过滤：保留小于40 Klbs的值    df = df[df['Weight on Bit'] < 40]
    # 扭矩过滤：保留小于4900 ft-lbs的值    df = df[df['Rotary Torque'] < 4900]
    # 差压过滤：保留190-770 psi之间的值    df = df[(df['Differential Pressure'] > 190) & (df['Differential Pressure'] < 770)]
    # 钻速过滤：保留小于385 ft/hr的值    df = df[df['Rate Of Penetration'] < 385]    return df
# 应用异常值处理df_clean = process_outliers(df)
# Box-Cox变换函数（处理数据偏度）from scipy.stats import boxcoxdef boxcox_transform(df, column_name):    # 确保所有值为正（Box-Cox要求）    df[column_name] = df[column_name] - df[column_name].min() + 0.1      # 应用Box-Cox变换    transformed, _ = boxcox(df[column_name])    df[column_name] = transformed  # 替换原数据    return df
# 对钻压和扭矩进行Box-Cox变换df_clean = boxcox_transform(df_clean, 'Weight on Bit')df_clean = boxcox_transform(df_clean, 'Rotary Torque')
# 平均比能(MSE)特征工程import mathdef create_mse(df):    mse_values = []  # 存储MSE计算结果    hole_diameter = 8.5  # 井眼直径(英寸)
    # 遍历每行数据计算MSE    for i in range(len(df)):        rpm = df['Rotary RPM'].iloc[i]  # 转速        torque = df['Rotary Torque'].iloc[i]  # 扭矩        rop_avg = df['Rate Of Penetration'].iloc[i]  # 平均钻速        wob = df['Weight on Bit'].iloc[i]  # 钻压
        # MSE计算公式（避免除零错误）        if rop_avg != 0:            mse = ((480 * torque * rpm) / (hole_diameter**2 * rop_avg)) + \                  ((4 * wob) / (hole_diameter**2 * math.pi))        else:            mse = 0        mse_values.append(mse)
    df['MSE'] = mse_values  # 添加MSE列    return df
# 应用MSE特征工程df_clean = create_mse(df_clean)
# 特征选择selected = ['Rate Of Penetration', 'Rotary Torque', 'Rotary RPM', 'Weight on Bit',             'Differential Pressure', 'Hook Load', 'MSE', 'Hole Depth']df_final = df_clean[selected]
# 数据标准化和相关性分析from sklearn.preprocessing import StandardScalerscaler = StandardScaler()datanorm = scaler.fit_transform(df_final)
# 计算斯皮尔曼相关系数from scipy import statsrho, pval = stats.spearmanr(datanorm)
# 数据分割（训练集70%，测试集30%）from sklearn.model_selection import train_test_splity = df_final[['Rate Of Penetration']]  # 目标变量（钻速）X = df_final.drop(['Rate Of Penetration'], axis=1)  # 特征变量X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.30, random_state=1000)
# 特征归一化（MinMax缩放）from sklearn import preprocessingX_train_scaler = preprocessing.MinMaxScaler(feature_range=(0, 1))X_train = X_train_scaler.fit_transform(X_train)y_train_scaler = preprocessing.MinMaxScaler(feature_range=(0, 1))y_train = y_train_scaler.fit_transform(y_train)
# XGBoost模型构建from xgboost import XGBRegressorxgb_model = XGBRegressor(    objective='reg:squarederror',  # 回归任务    n_estimators=200,  # 树的数量    reg_lambda=1,  # L2正则化    max_depth=3,  # 树的最大深度    learning_rate=0.1,  # 学习率    reg_alpha=0.1  # L1正则化)
# 模型训练xgb_model.fit(X_train, y_train)
# 模型评估from sklearn.metrics import r2_score, mean_absolute_error, mean_squared_errory_pred = xgb_model.predict(X_test)
# 计算评估指标r2 = r2_score(y_test, y_pred)mae = mean_absolute_error(y_test, y_pred)mse = mean_squared_error(y_test, y_pred)rmse = np.sqrt(mse)
# 特征重要性可视化feature_imp = pd.Series(xgb_model.feature_importances_, index=X.columns)feature_imp.sort_values().plot(kind='barh')
# 模型保存与加载import joblibjoblib.dump(xgb_model, 'rop_model_xgb.pkl')  # 保存模型loaded_model = joblib.load('rop_model_xgb.pkl')  # 加载模型
# 新数据预测new_data = pd.read_csv('freshdrillingdata.csv')  # 加载新数据# 注意：新数据需要与训练数据相同的预处理predicted_rop = loaded_model.predict(new_data)  # 预测ROP

https://www.zhihu.com/consult/people/792359672131756032?isMe=1