Skip to article frontmatterSkip to article content

Calcul de Trajectoire de Forage par la Méthode des Points Milieux

🎯 Objectif

Ce script a pour but de calculer et de visualiser la trajectoire 3D d’un forage en se basant sur une série de levés de déviation. La méthode employée est celle des points milieux, une technique standard pour modéliser le parcours d’un forage qui n’est jamais parfaitement droit.

⚙️ Fonctionnement

L’outil fonctionne en trois étapes simples :

  1. Saisie des données : L’utilisateur fournit les levés de déviation pour chaque station mesurée le long du forage :

    • Azimut : L’orientation de la boussole (direction horizontale).
    • Plongée : L’angle d’inclinaison par rapport à l’horizontale.
    • Profondeur mesurée : La distance le long du trou de forage jusqu’à la station.
  2. Calcul de la trajectoire : Le programme calcule les coordonnées cartésiennes (X,Y,ZX, Y, Z) de chaque station en appliquant l’algorithme des points milieux.

  3. Création de composites : Enfin, il est possible de calculer la position des “composites” : le script calcule les coordonnées 3D de points situés à un intervalle de longueur constant (LL) le long de la trajectoire du forage. Les coordonnées de chaque composite sont retournées en sortie.

Source
import numpy as np
import plotly.graph_objs as go
import ipywidgets as widgets
from IPython.display import display, clear_output, Markdown

def d2r(deg):
    return np.deg2rad(deg)

def compute_segments_midpoint(measurements, start_pos=(0,0,0), show_debug=False):
    n = len(measurements)
    depths = [m[0] for m in measurements]

    seg_lengths = []
    for i in range(n):
        if i == 0:
            seg_len = (depths[1] - depths[0]) / 2
        elif i == n - 1:
            seg_len = (depths[-1] - depths[-2]) / 2
        else:
            seg_len = (depths[i] - depths[i-1]) / 2 + (depths[i+1] - depths[i]) / 2
        seg_lengths.append(seg_len)

    positions = [start_pos]
    mids = []

    for i in range(n):
        md = seg_lengths[i]
        inc, azi = measurements[i][1], measurements[i][2]  # Correct order
        azi_rad, inc_rad = d2r(azi), d2r(inc)

        dx = md * np.cos(inc_rad) * np.sin(azi_rad)
        dy = md * np.cos(inc_rad) * np.cos(azi_rad)
        dz = -md * np.sin(inc_rad)

        last = positions[-1]
        new_pos = (last[0] + dx, last[1] + dy, last[2] + dz)
        positions.append(new_pos)

        md_mid = depths[i]
        mid_x = last[0] + dx / 2
        mid_y = last[1] + dy / 2
        mid_z = last[2] + dz / 2
        mids.append((md_mid, mid_x, mid_y, mid_z))

        if show_debug:
            print(f"Segment {i+1} : md segment={md:.2f} m, azimut={azi:.1f}°, inclinaison={inc:.1f}°")
            print(f"  dx={dx:.2f}, dy={dy:.2f}, dz={dz:.2f}")
            print(f"  Position finale segment: x={new_pos[0]:.2f}, y={new_pos[1]:.2f}, z={new_pos[2]:.2f}\n")

    return positions, mids

def interpolate_composites(measurements, positions, targets):
    depths = [m[0] for m in measurements]
    interp_XYZ = []

    for t in targets:
        if t <= depths[0]:
            inc, azi = measurements[0][1], measurements[0][2]
            azi_rad, inc_rad = d2r(azi), d2r(inc)
            md_extra = t - depths[0]
            x0, y0, z0 = positions[0]
            xt = x0 + md_extra * np.cos(inc_rad) * np.sin(azi_rad)
            yt = y0 + md_extra * np.cos(inc_rad) * np.cos(azi_rad)
            zt = z0 - md_extra * np.sin(inc_rad)
            interp_XYZ.append((t, xt, yt, zt))

        elif t >= depths[-1]:
            inc, azi = measurements[-1][1], measurements[-1][2]
            azi_rad, inc_rad = d2r(azi), d2r(inc)
            md_extra = t - depths[-1]
            x0, y0, z0 = positions[-1]
            xt = x0 + md_extra * np.cos(inc_rad) * np.sin(azi_rad)
            yt = y0 + md_extra * np.cos(inc_rad) * np.cos(azi_rad)
            zt = z0 - md_extra * np.sin(inc_rad)
            interp_XYZ.append((t, xt, yt, zt))

        else:
            for i in range(1, len(depths)):
                if depths[i-1] <= t <= depths[i]:
                    f = (t - depths[i-1]) / (depths[i] - depths[i-1])
                    x1, y1, z1 = positions[i]
                    x2, y2, z2 = positions[i+1]
                    xt = x1 + f * (x2 - x1)
                    yt = y1 + f * (y2 - y1)
                    zt = z1 + f * (z2 - z1)
                    interp_XYZ.append((t, xt, yt, zt))
                    break

    return interp_XYZ

def plotly_3d_path(positions, mids, interp, measurements):
    pos_arr = np.array(positions)
    fig = go.Figure()

    fig.add_trace(go.Scatter3d(
        x=pos_arr[:,0], y=pos_arr[:,1], z=pos_arr[:,2],
        mode='lines+markers',
        line=dict(color='blue', width=5),
        marker=dict(size=5, color='blue'),
        name='Trajectoire'
    ))

    md_vals = [m[0] for m in measurements]
    stations_coords = pos_arr[1:]
    fig.add_trace(go.Scatter3d(
        x=stations_coords[:,0], y=stations_coords[:,1], z=stations_coords[:,2],
        mode='markers+text',
        marker=dict(size=7, color='orange'),
        text=[f"MD {md:.1f} m" for md in md_vals],
        textposition='top center',
        name='Stations (mesures)'
    ))

    if mids:
        mids_np = np.array(mids)
        fig.add_trace(go.Scatter3d(
            x=mids_np[:,1], y=mids_np[:,2], z=mids_np[:,3],
            mode='markers+text',
            marker=dict(size=6, color='green'),
            text=[f"MD milieu {md:.1f} m" for md in mids_np[:,0]],
            textposition='bottom center',
            name='Points milieux'
        ))

    if interp:
        interp_np = np.array(interp)
        fig.add_trace(go.Scatter3d(
            x=interp_np[:,1], y=interp_np[:,2], z=interp_np[:,3],
            mode='markers+text',
            marker=dict(size=6, color='red'),
            text=[f"MD {md:.1f} m" for md in interp_np[:,0]],
            textposition='top center',
            name='Composites interpolés'
        ))

    fig.update_layout(
        scene=dict(
            xaxis_title='X (m)',
            yaxis_title='Y (m)',
            zaxis_title='Z (m, ↑)',
            aspectmode='data',
        ),
        title='Trajectoire 3D par méthode des points milieux',
        margin=dict(l=0, r=0, b=0, t=40),
        legend=dict(x=0.05, y=0.95)
    )
    fig.show()

# Widgets
num_points_w = widgets.BoundedIntText(value=3, min=2, max=50, description="Nombre stations:")
md_list_w = widgets.Text(value="0, 80, 180", description="MD (m):")
azi_list_w = widgets.Text(value="60, 55, 45", description="Azimut (\u00b0):")
inc_list_w = widgets.Text(value="90, 80, 70", description="Inclinaison (\u00b0):")
targets_w = widgets.Text(value="180, 183, 186", description="Composites (MD):")
collet_w = widgets.Text(value="0, 0, 0", description="Collet (X,Y,Z):")
show_debug_w = widgets.Checkbox(value=False, description="Afficher calculs intermédiaires")
button = widgets.Button(description="Calculer et Tracer", button_style='success', icon='play')
out = widgets.Output()

def on_button_clicked(b):
    with out:
        clear_output()
        try:
            md_list = list(map(float, md_list_w.value.strip().split(',')))
            azi_list = list(map(float, azi_list_w.value.strip().split(',')))
            inc_list = list(map(float, inc_list_w.value.strip().split(',')))
            collet = tuple(map(float, collet_w.value.strip().split(',')))

            targets_str = targets_w.value.strip()
            targets = list(map(float, targets_str.split(','))) if targets_str else []

            n = num_points_w.value
            if not (len(md_list) == len(azi_list) == len(inc_list) == n):
                print("⚠️ Erreur : Le nombre de valeurs pour MD, azimut et inclinaison doit correspondre au nombre de stations.")
                return

            measurements = list(zip(md_list, azi_list, inc_list))
            positions, mids = compute_segments_midpoint(measurements, start_pos=collet, show_debug=show_debug_w.value)
            interp = interpolate_composites(measurements, positions, targets)

            if interp:
                interp.sort(key=lambda x: x[0])
                md_output = "### Coordonnées des Composites\n\n"
                md_output += "| Profondeur (MD) | X | Y | Z |\n"
                md_output += "|:---:|:---:|:---:|:---:|\n"
                for t, x, y, z in interp:
                    md_output += f"| {t:.2f} m | {x:.2f} | {y:.2f} | {z:.2f} |\n"
                display(Markdown(md_output))
            elif targets:
                print("ℹ️ Aucun composite n'a pu être calculé.")

            plotly_3d_path(positions, mids, interp, measurements)

        except Exception as e:
            print(f"❌ Erreur lors de l'exécution : {str(e)}")

button.on_click(on_button_clicked)

info_text = widgets.Output()

with info_text:
    display(Markdown("""
    ### ℹ️ Conventions et paramètres
    - **Azimut (°)** : 0-360°, mesuré depuis le Nord dans le sens horaire.
    - **Inclinaison (°)** : angle entre le plan XY et l'axe Z descendant. 90° = forage vertical descendant.
    - **MD (m)** : profondeur mesurée le long du forage.
    - **Composites (MD)** : profondeurs où calculer les coordonnées interpolées.
    - **Collet (X,Y,Z)** : coordonnées du départ du forage.
    ---
    """))

# Regrouper les widgets d'entrée dans une VBox (ou HBox selon le besoin)
inputs = widgets.VBox([num_points_w, md_list_w, azi_list_w, inc_list_w, targets_w, collet_w, show_debug_w, button])

# Regrouper tout (infos + widgets + sortie)
main_ui = widgets.VBox([info_text, inputs, out])

display(main_ui)
Loading...

Exemple réaliste avec 100 stations

📝 Note : En construction