Premessa: dopo questo post mi sono lanciato in questa avventura nei frattali nella geologia non avendo nessuna nei guai in cui mi sarei infilato, di quando fosse malevola la natura che non solo e' frattale ma multifrattale e come sarei stato inadeguato dal punto di vista matematico. Detto questo l'argomento e' cosi complesso che viene diviso in due post,il presente ed il successivo
-------------------------------------
Sempre leggendo Fractals and Chaos in Geology and Geophysics 2nd edition -- Donald
Lawson Turcotte -- 2nd ed_, Cambridge, U_K, New York, England, 1997 si ha che le misure geofisiche seguono una legge frattale
Per fare una prova mi sono scaricato un log di pozzo da KGS
I dati sono in formato LAS che non indicano il formato per laserscanner ma un formato testuale per log di pozzo per possono essere letti con la libreria lasio di Python
import lasio
# Load LAS file
las = lasio.read("1055385870.las")
# Convert to pandas DataFrame
df = las.df()
# Save to CSV
df.to_csv("well_log.csv", index=True)
io ho scaricato del tutto a caso il file 1055385870.las queste le prime righe
~Version Information
VERS. 2.0: CWLS LOG ASCII STANDARD -VERSION 2.0
WRAP. NO: ONE LINE PER DEPTH STEP
~WELL INFORMATION BLOCK
#MNEM.UNIT DATA DESCRIPTION OF MNEMONIC
#--------- ------------- ------------------------
STRT.FT 0.300 :START DEPTH
STOP.FT 1279.300 :STOP DEPTH
STEP.FT 0.100 :STEP UP_HOLE
NULL. -999.25 :NULL VALUE
COMP. LEIS, VICTOR J. :COMPANY
WELL. MENTZER #N-1 :WELL
FLD . NEOSHO FALLS - LEROY :FIELD\LOCATION
LOC . 26 23S 17E :LOCATION
CNTY. ALLEN :COUNTY
STAT. KANSAS :STATE
CTRY. USA :COUNTRY
SRVC. :SERVICE COMPANY
DATE. 05/06/2011 :LOG DATE
API . 15- 001-30189-00-00 :API NUMBER
LIC . :LICENSE NUMBER
~Curve Information Block
#MNEM.UNIT API CODE Curve Description
#--------- ------------- ------------------------
DEPT .FT 00 001 00 00 : 1 DEPTH
COND(DI) .MMHO/M 00 000 00 00 : 2 DEEP INDUCTION COND
RES(GRD) .OHM-M 00 220 00 00 : 3 GUARD RESISTIVITY
SP .MV 00 010 00 00 : 4 SPON POTENTIAL
RES(DI) .OHM-M 00 120 00 00 : 5 DEEP INDUCTION RES
RES(MI) .OHM-M 00 120 00 00 : 6 MED INDUCTION RES
GAMMA .API-GR 00 310 00 00 : 7 GAMMA RAY
COMP .G/CC 00 356 00 00 : 8 DEN COMPENSATION
CALIPER .INCH 00 280 00 00 : 9 SHORT ARM CALIPER
DEN(CDL) .G/CC 00 350 00 00 : 10 COMPENSATED DENSITY
TEMP .DEG F 00 660 00 00 : 11 TEMPERATURE
POR(DEN) .PERCENT 00 890 00 00 : 12 DENSITY POROSITY
DEL TEMP .DEG F 00 000 00 00 : 13 SAMPLE DELTA TEMP
~Parameter Information Block
#MNEM.UNIT Information Description
#---------- --------------- ------------------------
FILE. PROCESSED :File Type
FIID. 9541 :File Type Identifier
VERS. 3.60B :System Version
SER . 1 :System Serial Number
TRUK. .3879 :Truck Calibration Number
Ho scelto di lavorare sul parametro Gamma.
Ci sono metodi per evidenziare la frattalita'
1) e' basato sull'analisi della frequenza delle classi di una misura
2) e' basata S(f)∝f^-β
Le ipotesi di base per la verifica della frattalita' sono
a) non e' necessariamente stazionario, puo' avere un trend e la varianza puo' crescere nel tempo ma lo sono le differenze tra X(t) e X(t+1). Cio' implica l'autosomiglianza ovvero lo stesso comportamento a diverse scale
b) deve sottostare ad un legge di potenza. Questa puo ' essere intepretata come
- frequenza come potenza N(>x)∝x^(-alpha). la somma di elementi superiori ad una soglia e' direttamente proporzionale al valore della soglia alla potenza di -alpha
- power spectral power : diretta proporzionalita' tra la potenza di spettro ad una determinata frequenza e la frequenza elevata alla potenze di -Beta PSD(f)∝f^-β
Proviamo con i dati reali partendo prima dalla frequenza
Nel grafico soprastante ci sono le classi con il numero di misure inferiori ad una determinata soglia
Se si plotta in grafico log-log si ha un chiaro trend lineare nella parte finale del grafico
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
# ======================================================
# 1 — LETTURA FILE CSV E COLONNA DI INTERESSE
# ======================================================
filename = "well_log.csv"
df = pd.read_csv(filename)
# Usa la colonna RES(GRD) → indice 2
col_index = 6
x = df.iloc[:, col_index].values
print(x)
# RIMUOVI i NaN (prime 100 righe vuote)
x = x[~np.isnan(x)]
if len(x) < 10:
raise ValueError("Troppi valori NaN in RES(GRD). Servono più dati validi.")
# ======================================================
# 2 — RISCALAMENTO SE CI SONO VALORI NEGATIVI
# ======================================================
xmin_original = np.min(x)
if xmin_original <= 0:
shift = abs(xmin_original) + 1e-6
x = x + shift
print(f"Riscalamento effettuato: aggiunto shift = {shift:.6f}")
# ======================================================
# 3 — CREAZIONE DI 25 CLASSI TRA MINIMO E MASSIMO
# ======================================================
xmin = np.min(x)
xmax = np.max(x)
thresholds = np.linspace(xmin, xmax, 25)
# ======================================================
# 4 — PER OGNI SOGLIA, CONTA GLI ELEMENTI >= SOGLIA
# ======================================================
counts = np.array([np.sum(x >= thr) for thr in thresholds])
# FILTRO valori validi
valid = np.logical_and(counts > 0, thresholds > 0)
thresholds = thresholds[valid]
counts = counts[valid]
# controlla se sono rimasti valori
if len(thresholds) < 3:
raise ValueError("Troppi valori invalidi dopo il filtraggio. Controlla la colonna RES(GRD).")
# ======================================================
# 5 — LOGARITMI
# ======================================================
logT = np.log10(thresholds)
logN = np.log10(counts)
# ======================================================
# 6 — REGRESSIONE SOLO PER log(threshold) ∈ [2.25, 2.5]
# ======================================================
mask = (logT >= 1.9) & (logT <= 4)
if np.sum(mask) >= 2:
logT_fit = logT[mask].reshape(-1, 1)
logN_fit = logN[mask]
model = LinearRegression()
model.fit(logT_fit, logN_fit)
slope = model.coef_[0]
intercept = model.intercept_
print("\n=== RISULTATI REGRESSIONE ===")
print("Intervallo usato: log10(soglia) ∈ [1.9, 4.0]")
print(f"Coefficiente angolare (slope) = {slope}")
print(f"Intercetta = {intercept}")
logN_pred = model.predict(logT_fit)
else:
print("\n⚠ Nessun dato nel range log10(t) ∈ [1.9, 4.0].")
print("La regressione NON è stata eseguita.")
logT_fit, logN_fit, logN_pred = None, None, None
# ======================================================
# 7 — ISTOGRAMMA DELLE CLASSI
# ======================================================
plt.figure(figsize=(10,5))
plt.bar(range(len(counts)), counts, width=0.8)
plt.xlabel("Classe (0–24) Parametro Gamma")
plt.ylabel("Numero elementi ≥ soglia della classe")
plt.title("Istogramma delle 25 classi Parametro Gamma")
plt.grid(axis='y')
plt.show()
# ======================================================
# 8 — GRAFICO LOG–LOG CON REGRESSIONE
# ======================================================
plt.figure(figsize=(10,6))
plt.scatter(logT, logN, s=20, label="Dati (log-log)")
if logT_fit is not None:
plt.plot(logT_fit.flatten(), logN_pred, linewidth=2,
label="Regressione (2.25–2.5)")
plt.xlabel("log10(soglia)")
if logT_fit is not None:
# Plot regressione
plt.plot(logT_fit.flatten(), logN_pred, linewidth=2,
label="Regressione (2.25–2.5)")
# Testo sul grafico con slope (3 cifre significative)
txt = f"slope = {slope:.3g}"
# posizione: leggermente a destra del primo punto valido
x_pos = logT_fit.min()
y_pos = logN_pred.max()
plt.text(x_pos, y_pos, txt, fontsize=12, color="red")
plt.ylabel("log10(N ≥ soglia)")
plt.title("Analisi log–log con regressione su intervallo selezionato")
plt.legend()
plt.grid(True)
plt.show()
usando invece l'analisi di densita' spettrale
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from numpy.fft import rfft, rfftfreq
from scipy.signal import welch
# ======================================================
# CONFIGURAZIONE
# ======================================================
CSV_FILE = "well_log.csv" # <-- cambia questo
WINDOW = 400 # sliding window
STEP = 1
COLUMN_DEPTH = 0
COLUMN_RXORT = 6
# ======================================================
# 1. LETTURA CSV
# ======================================================
df = pd.read_csv(CSV_FILE)
depth = df.iloc[:, COLUMN_DEPTH].values
rxort = df.iloc[:, COLUMN_RXORT].values
# ======================================================
# 2. FILTRO VALORI INVALIDI
# ======================================================
bad_values = [-9999, -999.25]
valid_mask = (
~np.isnan(depth) &
~np.isnan(rxort) &
np.isfinite(depth) &
np.isfinite(rxort) &
~np.isin(rxort, bad_values)
)
depth = depth[valid_mask]
rxort = rxort[valid_mask]
print("Campioni dopo il filtro:", len(rxort))
# ======================================================
# 3. FREQUENZA DI CAMPIONAMENTO
# ======================================================
dz = np.mean(np.diff(depth))
fs = 1.0 / dz # samples per meter
# ======================================================
# 4. SPETTROGRAMMA FFT A FINESTRA MOBILE
# ======================================================
num_samples = len(rxort)
freqs = rfftfreq(WINDOW, d=1/fs)
spectra = []
centers = []
for i in range(0, num_samples - WINDOW, STEP):
segment = rxort[i:i+WINDOW]
segment = segment - np.mean(segment)
F = rfft(segment)
P = np.abs(F)**2
spectra.append(P)
centers.append(depth[i + WINDOW//2])
spectra = np.array(spectra)
centers = np.array(centers)
# ======================================================
# 5. PLOT SPETTROGRAMMA + LOG ORIGINALE
# ======================================================
fig = plt.figure(figsize=(12, 8))
gs = fig.add_gridspec(2, 1, height_ratios=[2, 1], hspace=0.25)
# ---- (A) Spectrogram ----
ax1 = fig.add_subplot(gs[0, 0])
im = ax1.imshow(
spectra.T,
aspect='auto',
origin='lower',
extent=[centers.min(), centers.max(), freqs.min(), freqs.max()]
)
ax1.set_ylabel("Frequency (cycles/m)")
ax1.set_title("Spectrogram (window = 400 samples)")
# plt.colorbar(im, ax=ax1).set_label("Power")
# ---- (B) Log ----
ax2 = fig.add_subplot(gs[1, 0], sharex=ax1)
ax2.plot(depth, rxort, linewidth=0.8)
ax2.set_xlabel("Depth (m)")
ax2.set_ylabel("")
ax2.grid(True)
plt.setp(ax1.get_xticklabels(), visible=False)
plt.show()
# ======================================================
# 6. PSD WELCH + FRACTAL FIT
# ======================================================
nper = min(1024, len(rxort))
freqs_psd, psd = welch(rxort, fs=fs, nperseg=nper)
print("PSD samples:", len(freqs_psd))
print("Freq range:", freqs_psd.min(), " - ", freqs_psd.max())
# ---- Intervallo utile per fit ----
mask = (freqs_psd > 0.001) & (freqs_psd < 0.2)
print("Punti nel fit:", np.sum(mask))
if np.sum(mask) < 2:
raise ValueError("Intervallo di frequenza troppo stretto o dati insufficienti!")
# ---- Fit lineare in scala log-log ----
p = np.polyfit(np.log(freqs_psd[mask]), np.log(psd[mask]), 1)
beta = -p[0]
H = (beta - 1) / 2
D = 2 - H
print("\n========== RISULTATI FRACTAL ==========")
print("Spectral β =", beta)
print("Hurst H =", H)
print("Fractal D =", D)
# ======================================================
# 7. PLOT PSD + FIT
# ======================================================
plt.loglog(freqs_psd, psd, label="PSD")
plt.loglog(freqs_psd[mask], np.exp(np.polyval(p, np.log(freqs_psd[mask]))),
label="Fit", linewidth=2)
plt.xlabel("Frequency (cycles/m)")
plt.ylabel("PSD")
plt.title("Spectral fractal analysis")
plt.legend()
plt.show()
Spectral β = -0.594 (pendenza della retta di approssimazione)
Per questo valore negativo indica che le basse frequenze hanno un maggior potenza rispetto alle alte. Dal punto di vista assoluto 0.6 e' molto vicino a zero ad indicare che i dati sono domininati dal rumore bianco.
L'esponente di Hurst viene calcolato come H=(β +1)/2 quindi 0.797 e la dimensione frattale D=2-β pari a 1.2
La varazione dei parametri frattali e' compresa
0<H<1
1<D<2
-1<β<3
Un ulteriore approccio e' quello di utilizzare l'analisi wavelet
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pywt
from scipy.stats import linregress
# ======================================================
# 1 — LETTURA FILE CSV E COLONNA DI INTERESSE
# ======================================================
filename = "well_log.csv"
df = pd.read_csv(filename,nrows=2500)
# Utilizzo la colonna GAMMA (indice 6) come nei codici precedenti
x = df.iloc[:, 6].values # Dati di log (GAMMA)
# ======================================================
# 2 — PULIZIA DATI E RISCALAMENTO
# ======================================================
# Converti in numerico e gestisci NaN
x = pd.to_numeric(x, errors='coerce')
# Rimuovi i valori NaN.
x = x[~np.isnan(x)]
# Riscalamento se ci sono negativi
xmin = np.min(x)
if xmin <= 0:
shift = abs(xmin) + 1e-6
x = x + shift
print(f"Riscalamento effettuato: shift = {shift:.6f}")
# ======================================================
# 3 — PARAMETRI WAVELET
# ======================================================
wavelet = "cmor1.5-1.0"
scales = np.arange(1, 100)
# ======================================================
# 4 — CALCOLO CWT e POTENZA
# ======================================================
coeffs, frequencies = pywt.cwt(x, scales, wavelet, sampling_period=1)
power = np.abs(coeffs)**2
# ======================================================
# 5 — CALCOLO DELL'ESPONENTE DI HURST (H)
# ======================================================
# Calcola la Potenza Media su ciascuna scala (Varianza Wavelet)
scale_mean_power = np.mean(power, axis=1)
# Trasformazione log-log
log_scales = np.log10(scales)
log_mean_power = np.log10(scale_mean_power)
# Selezione della regione di scaling (dall'10% all'80% delle scale)
min_scale_index = int(len(scales) * 0.1)
max_scale_index = int(len(scales) * 0.8)
H_log_scales = log_scales[min_scale_index:max_scale_index]
H_log_mean_power = log_mean_power[min_scale_index:max_scale_index]
# Regressione Lineare
slope, intercept, r_value, p_value, std_err = linregress(H_log_scales, H_log_mean_power)
# Calcolo di H
hurst_exponent = (slope + 1) / 2
fractal_dimension = 2 - hurst_exponent
print("-" * 50)
print("Risultati dell'Analisi di Frattalità Wavelet:")
print(f"Esponente di Hurst (H): {hurst_exponent:.4f}")
print(f"Dimensione Frattale (D): {fractal_dimension:.4f}")
print(f"Pendenza della regressione (2H - 1): {slope:.4f}")
print("-" * 50)
# ======================================================
# 6 — PLOT DEI RISULTATI (ALLINEAMENTO E RIMOZIONE COLORBAR)
# ======================================================
# NOTA: Usiamo sharex=True per allineare A e B. C (la regressione) deve essere separato.
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8),
sharex=True, gridspec_kw={'height_ratios': [3, 1]})
# --- A) SCALOGRAMMA WAVELET (SENZA SCALA DI POTENZA)
extent = [0, len(x), 0, len(scales)]
im = ax1.imshow(power, extent=extent, cmap="jet", aspect="auto", origin="lower")
ax1.set_title(f"Scalogramma Wavelet (Potenza per Scala)")
ax1.set_ylabel("Scala Wavelet")
# Rimozione della colorbar (scala di potenza)
# plt.colorbar(im, ax=ax1, orientation='vertical', label='Potenza') <-- RIMOSSA
# --- B) PLOT XY SOTTO (SEGNALE ORIGINALE)
ax2.plot(x, color="black")
ax2.set_xlabel("Indice campione (profondità)")
ax2.set_ylabel("Valore (GAMMA)")
ax2.grid(True)
plt.tight_layout()
plt.show()
# --- C) REGRESSIONE PER H (PLOT SEPARATO E NON ALLINEATO)
# Questo grafico non ha l'asse X della profondità, quindi deve rimanere separato
plt.figure(figsize=(8, 6))
plt.loglog(scales, scale_mean_power, 'o', markersize=3, label='Potenza Media')
plt.loglog(scales[min_scale_index:max_scale_index],
intercept + slope * H_log_scales,
'r--', label=f'Regressione: pendenza={slope:.2f}')
plt.xlabel('Scala Wavelet (log)')
plt.ylabel('Potenza Media (log)')
plt.title(f'Calcolo dell\'Esponente di Hurst (H={hurst_exponent:.4f})')
plt.legend()
plt.grid(True, which="both", ls="--")
plt.show()
Nessun commento:
Posta un commento