In questo post vengono analizzati i valori di torbidita' della sorgente del Frigido (periodo 13/02/2020 01/12/2025 dati orari ) e Carbonera (26/07/2006 21/06/2013 dati ogni 3 minuti ) legate all'acquiferi carsico delle Alpi Apuane
 |
| Torbidita' del tempo Frigido |
Il parametro Esponente di Hurst e' stato introdotto Hurst, H. E. (1951). Long-term storage capacity of reservoirs. Transactions of the American Society of Civil Engineers, 116(1), 770-808.
E' una misura dell'inerzia del sistema idrogeologico (ovvero capacita' di stoccaggio) e della complessita' della rete carsica. Per i valore di 0.82 e 0.85 trovati si ha che il sistema non e' casuale (H=0.5), non e' deterministicamente semplice ed il comportamento e' correlato su differenti scale di tempo ed e' un sistema che ha "memoria" del proprio passato
Non e' possibile prevedere in maniera completamente deterministica il valore esatto del parametro indagato
 |
| Frigido |
Con modelli deterministici si avra' un discreto successo su scale di tempo brevi ma non sulla previsione a lunga scadenza. Il passaggio da regime deterministico a frattale e' indicato come tempo di crossover e si individua dai punti di flesso sull'andamento dei punti. Per Carbonera (grafico sottostante) il punto di cross over e' tra 1000 e 10000 unita' che tradotto nel passo di campionamento si arriva ad intervallo da le 5 e le 50 ore
 |
| Carbonera |
Valori di H vicini ad 1 indicano persistenza nei valori (sia che siano massimi che minimi) e questo puo' incidere anche sul valore massimo e non solo sulle code della distribuzione
La memoria non si esaurisce mai ma ha una legge di decadimento come una potenze (in questo caso 2*0.84 -2 = -0.32
Il fatto che il calcolo di H risulti sostanzialmente identico su serie con diverso passo di campionamento e' un indice della frattalita' del fenomeno osservato (e' indipendente dalla scala del tempo di osservazione)
il valore di Hurst prossimo ad 1 indica che e' probabile che la torbidita' aumenti nel prossimo lag (passo di campionamento) e' fortemente dipendente da cio' che accade nel lag precedente . Cio' sarebbe in contrasto con il concetto di massimo relativo della curva di scarica ma la statistica deve essere considerata valida al mantenersi delle condizioni al contorno e delle relative forzanti. Con un valore di 0.82 la coda di esaurimento sara' lunga
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# --- 1. FUNZIONE DEFINIZIONE per l'ANALISI R/S MANUALE ---
def compute_rs_data(series):
"""
Calcola la statistica Rango Riscalato (R/S) per una serie temporale.
Restituisce gli intervalli temporali (n) e i valori medi di R/S.
"""
N = len(series)
# Definisce gli intervalli temporali (n) in progressione geometrica fino a N/4
n_values = np.unique(np.logspace(1, np.log10(N // 4), num=50, dtype=int))
n_values = n_values[n_values >= 1]
rs_values = []
for n in n_values:
if n == 0: continue
num_blocks = N // n
rs_per_block = []
for i in range(num_blocks):
block = series[i * n : (i + 1) * n]
mean = np.mean(block)
std_dev = np.std(block)
if std_dev == 0:
continue
X = np.cumsum(block - mean)
R = np.max(X) - np.min(X)
rs_per_block.append(R / std_dev)
if rs_per_block:
rs_values.append(np.mean(rs_per_block))
else:
break
return n_values[:len(rs_values)], np.array(rs_values)
# --- 2. CARICAMENTO E PREPARAZIONE DEI DATI ADATTATI ---
FILE_NAME = "frigido.csv"
COLUMN_NAME = "Torbidita NTU"
FILTER_N_THRESHOLD = 100
try:
# ATTENZIONE: Caricamento con separatore PUNTO E VIRGOLA (;) e decimale VIRGOLA (,)
df = pd.read_csv(FILE_NAME, sep=';', decimal=',')
# Seleziona la colonna 'Torbidita NTU'
discharge_series = df[COLUMN_NAME]
# Conversione in numerico e gestione dei NaN (valori mancanti o non validi)
discharge_series = pd.to_numeric(discharge_series, errors='coerce').dropna()
# Converti in array NumPy
discharge_data = discharge_series.values
except FileNotFoundError:
print(f"Errore: Il file '{FILE_NAME}' non è stato trovato. Assicurati che sia nella directory corretta.")
exit()
except KeyError:
print(f"Errore: La colonna '{COLUMN_NAME}' non è stata trovata nel file. Controlla il nome esatto.")
print(f"Colonne disponibili: {df.columns.tolist()}")
exit()
if len(discharge_data) < 100:
print("Errore: Non ci sono abbastanza punti dati validi per eseguire l'analisi R/S.")
exit()
# --- 3. CALCOLO R/S ---
rs_n, rs_rs = compute_rs_data(discharge_data)
# Calcolo di H (Esponente di Hurst) per l'INTERVALLO COMPLETO
log_n_full = np.log10(rs_n)
log_rs_full = np.log10(rs_rs)
coefficients_full = np.polyfit(log_n_full, log_rs_full, 1)
H_full = coefficients_full[0]
c_full = 10**coefficients_full[1]
# --- 4. PLOT: ANALISI R/S INTERVALLO COMPLETO ---
plt.figure(figsize=(8, 5))
plt.plot(rs_n, c_full * rs_n**H_full, color="deepskyblue", label=f'Pendenza Frattale (H={H_full:.2f})')
plt.scatter(rs_n, rs_rs, color="navy", alpha=0.5, label='Dati R/S empirici')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('Intervallo Temporale (n)')
plt.ylabel('Rango Riscalato R/S')
plt.title(f'Scaling Frattale della Torbidità (Full Range)')
plt.legend()
plt.savefig('Frigido_Full_R_S_Plot.png')
# --- 5. PLOT: INTERVALLO FILTRATO (n > 100) ---
filter_mask = rs_n > FILTER_N_THRESHOLD
X_filtered = rs_n[filter_mask]
Y_filtered = rs_rs[filter_mask]
# Calcolo di H per il solo intervallo filtrato
if len(X_filtered) > 1:
log_n_filtered = np.log10(X_filtered)
log_rs_filtered = np.log10(Y_filtered)
coefficients_filtered = np.polyfit(log_n_filtered, log_rs_filtered, 1)
H_filtered = coefficients_filtered[0]
c_filtered = 10**coefficients_filtered[1]
plt.figure(figsize=(8, 5))
plt.scatter(X_filtered, Y_filtered, color="darkgreen", alpha=0.7, label=f'Dati R/S empirici (n > {FILTER_N_THRESHOLD})')
plt.plot(X_filtered, c_filtered * X_filtered**H_filtered, color="red", linestyle='--',
label=f'Pendenza Frattale Filtrata (H={H_filtered:.2f})')
plt.xscale('log')
plt.yscale('log')
plt.xlabel('Intervallo Temporale (n)')
plt.ylabel('Rango Riscalato R/S')
plt.title(f'Scaling Frattale della Torbidità (Intervallo Temporale > {FILTER_N_THRESHOLD})')
plt.legend()
plt.savefig('Frigido_Filtered_R_S_Plot_gt_100.png')
print("\n--- Risultati del Calcolo ---")
print(f"Hurst Esponent (Torbidità - Intervallo Completo): {H_full:.4f}")
print(f"Hurst Esponent (Torbidità - Filtrato n>{FILTER_N_THRESHOLD}): {H_filtered:.4f}")
else:
print(f"Nota: Non ci sono abbastanza punti dati rimasti dopo il filtraggio (n > {FILTER_N_THRESHOLD}) per generare il grafico filtrato.")
Nessun commento:
Posta un commento