Frattalita' in misure geofisiche di pozzo (2)

Nella conclusione del post precedente era che i vari tentativi di ricondurre i dati sperimentali all'analisi frattale davano risultati contradittori...tutto cio' fino sono arrivato a leggere Fractal Behaviour of the Earth System V.P.Dimitri  Springer 2005 pag 63-81 Regularity Analysis Applied to Well Log Data (M.Fedi, D.Fiore, Mauro La Mann) dove parlando proprio dei dati di KST viene introdotto il concetto di multifrattalita' (per essere precisi di multifrattalita' aveva gia' parlato Mandelbrot ed in D. Turcotte  Fractals and Chaos in Geology and Geophysics si trova applicazione dei multifrattali applicata alle geoscienze a pag)

Prima di tutto si puo' effettuare un test se i dati del pozzo possono essere inquadrati all'interno della multifrattalita' 

Questo puo' essere fatto l'analisi di regolarita' locale tramite analisi di fluttuazione a momenti multipli (Multifractal Detrended Fluctuation Analysis, MF-DFA) 

Vengono testati valori del parametro q da -5 a 5 (escludendo lo 0) verificando la pendenza h(q) 

 

 

Se il segnale fosse monofrattale la curva sarebbe una linea orizzontale ma non e' il caso. Quindi il test mostra la presenza di trend multifrattale con q>0 (destra) influenzato dalle grandi fluttuazioni (picchi del grafico dei valori di Gamma o contrasti gelogici) o q<0 (sinistra) influenzati da regioni piu' omogenee o da rumore di fondo

La linea rossa (q=2) indica  la persistenza media dell'intera serie

nel grafico sottostante viene descritta l'intensita' della multifrattalita' (se non fossero dati multi frattali si ridurrebbe ad un punto. Piu' ampia e' l'apertura della curva maggiore e' l'intensita' della frattalita'

I dati quindi non possono essere descritti con un singolo esponente di Hurst 

 

import numpy as np

import pandas as pd

import matplotlib.pyplot as plt

from scipy.stats import linregress

# ======================================================

# 0 — FUNZIONE MF-DFA MANUALE

# ======================================================

def mfdfa(series, qs, scales, order=1):

"""

Calcola l'esponente di Hurst generalizzato h(q) utilizzando l'algoritmo MF-DFA.

"""

N = len(series)

# 1. Integrazione del segnale (necessaria per il detrending)

Y = np.cumsum(series - np.mean(series))

h_q = []

for q in qs:

F_q_2 = [] # Momenti di fluttuazione di q-esimo ordine

for s in scales:

s = int(s)

# Suddivisione in segmenti di lunghezza s

segments = N // s

# Calcolo della varianza per tutti i segmenti

F_q_s_list = []

for i in range(segments):

idx = slice(i * s, (i + 1) * s)

Y_segment = Y[idx]

# 2. Detrending locale (regressione polinomiale di ordine 'order')

x_local = np.arange(s)

p = np.polyfit(x_local, Y_segment, order)

Y_fit = np.polyval(p, x_local)

# 3. Calcolo della varianza locale (residui al quadrato)

variance = np.mean((Y_segment - Y_fit)**2)

F_q_s_list.append(variance)

# 4. Calcolo della funzione di fluttuazione di q-esimo ordine

F_q_s_list = np.array(F_q_s_list)

if q != 2.0:

# Esecuzione del momento q-esimo (q/2 sulle varianze)

F_q = np.mean(F_q_s_list**(q / 2.0))**(1.0 / q)

else:

# Caso speciale q=2 (momento quadratico)

F_q = np.sqrt(np.mean(F_q_s_list))

F_q_2.append(F_q)

# 5. Regressione log-log: F_q(s) ~ s^h(q)

log_s = np.log10(scales)

log_F_q = np.log10(F_q_2)

# Filtra i valori infiniti o NaN che possono emergere da log(0)

valid_indices = np.isfinite(log_s) & np.isfinite(log_F_q)

if np.sum(valid_indices) > 1:

slope, _, _, _, _ = linregress(log_s[valid_indices], log_F_q[valid_indices])

h_q.append(slope)

else:

h_q.append(np.nan)

return np.array(h_q)

# ======================================================

# 1 — LETTURA FILE CSV E PREPARAZIONE DATI

# ======================================================

filename = "well_log.csv"

RIGHE_LIMITE = 12700

COLONNA_INDICE = 1

try:

df = pd.read_csv(filename, nrows=RIGHE_LIMITE)

except FileNotFoundError:

print(f"Errore: File '{filename}' non trovato.")

exit()

x = df.iloc[:, COLONNA_INDICE].values

x = pd.to_numeric(x, errors='coerce')

x = x[~np.isnan(x)]

# Riscalamento (opzionale)

xmin = np.min(x)

if xmin <= 0:

shift = abs(xmin) + 1e-6

x = x + shift

print(f"Riscalamento effettuato: shift = {shift:.6f}")

N = len(x)

print(f"Analisi eseguita su {N} campioni (Prime {RIGHE_LIMITE} righe del file).")

# ======================================================

# 2 — PARAMETRI MF-DFA

# ======================================================

# Scelta dell'intervallo q

Q_MIN = -5

Q_MAX = 5

Q_STEP = 0.5

qs = np.arange(Q_MIN, Q_MAX + Q_STEP, Q_STEP)

qs = qs[qs != 0] # Esclude q=0

# Scegliere un intervallo di scale logaritmico per una migliore regressione

min_scale = int(N / 50) # Scale minime per un detrending efficace (p.es. 50 campioni)

max_scale = int(N / 4) # Scala massima (1/4 della lunghezza totale)

scales = np.unique(np.logspace(np.log10(min_scale), np.log10(max_scale), num=50).astype(int))

ORDER = 1 # Ordine del polinomio di detrending locale

# ======================================================

# 3 — ESECUZIONE MF-DFA E CALCOLO h(q)

# ======================================================

print("Calcolo degli esponenti h(q)...")

hq = mfdfa(x, qs, scales, order=ORDER)

# Filtra i valori NaN (dovuti a una regressione fallita)

valid_hq = hq[np.isfinite(hq)]

valid_qs = qs[np.isfinite(hq)]

# ======================================================

# 4 — ANALISI E INTERPRETAZIONE DEI RISULTATI

# ======================================================

# 4a. Analisi monofattale (q=2)

if 2.0 in valid_qs:

q_mono_index = np.where(valid_qs == 2.0)[0][0]

hurst_monofractal = valid_hq[q_mono_index]

else:

# Interpolazione per q=2 se non è stato calcolato

hurst_monofractal = np.interp(2.0, valid_qs, valid_hq)

# 4b. Analisi della Multifrattalità

h_max = np.max(valid_hq)

h_min = np.min(valid_hq)

Delta_h = h_max - h_min

# ======================================================

# 5 — PLOT DEI RISULTATI

# ======================================================

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8), sharex=False)

# --- A) Grafico h(q) vs q (Test di Multifattalità)

ax1.plot(valid_qs, valid_hq, 'o-', color='darkblue', markersize=4)

ax1.axhline(hurst_monofractal, color='red', linestyle='--',

label=f'H monofattale (q=2): {hurst_monofractal:.3f}')

ax1.set_title('Esponente di Hurst Generalizzato $h(q)$')

ax1.set_xlabel('Momento di Fluttuazione $q$')

ax1.set_ylabel('Esponente di Scaling $h(q)$')

ax1.grid(True)

ax1.legend()

# --- B) Lo Spettro Multifrattale f(alpha) vs alpha (Singolarità)

# Calcolo dello Spettro tramite Trasformata di Legendre

tau_q = valid_qs * valid_hq - 1

alpha = np.gradient(tau_q, valid_qs)

f_alpha = valid_qs * alpha - tau_q

ax2.plot(alpha, f_alpha, 's-', color='green', markersize=4)

ax2.set_title('Spettro Multifrattale $f(\\alpha)$')

ax2.set_xlabel('Esponente di Singolarità $\\alpha$')

ax2.set_ylabel('Dimensione Frattale $f(\\alpha)$')

ax2.grid(True)

plt.tight_layout()

plt.show()

print("\n" + "=" * 50)

print(f"ANALISI MULTIFRATTALE (MF-DFA) - RISULTATI")

print(f"Δh = h_max - h_min = {Delta_h:.4f}")

if Delta_h > 0.1:

print("CONCLUSIONE: Il segnale è **fortemente multifrattale** (Δh > 0.1).")

print("La sua complessità non è uniforme: sono presenti strutture a diverse scale e diversi tipi di singolarità (anomalie e omogeneità).")

else:

print("CONCLUSIONE: Il segnale è prevalentemente **monofattale** (Δh <= 0.1).")

print(f"Esponente di Hurst Monofattale (q=2): H ≈ {hurst_monofractal:.3f}")

print("=" * 50)