.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/Ex_TL_inversion.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_Ex_TL_inversion.py: Ex. Time-Lapse ERT Inversion Techniques =================================== This example demonstrates different approaches for time-lapse electrical resistivity tomography (ERT) inversion using PyHydroGeophysX. The example includes: 1. Full time-lapse inversion with temporal regularization 2. Windowed time-lapse inversion for large datasets 3. L1-norm regularized inversion for sharp boundary recovery 4. Comparison of different inversion strategies 5. Visualization of resistivity evolution over time Time-lapse ERT inversion is crucial for monitoring subsurface water content changes and understanding hydrological processes in watersheds. The temporal regularization helps maintain consistency between consecutive time steps while allowing for realistic changes. .. GENERATED FROM PYTHON SOURCE LINES 20-45 .. code-block:: default import os import sys import numpy as np import matplotlib.pyplot as plt import pygimli as pg from pygimli.physics import ert from mpl_toolkits.axes_grid1 import make_axes_locatable # Setup package path for development try: # For regular Python scripts current_dir = os.path.dirname(os.path.abspath(__file__)) except NameError: # For Jupyter notebooks current_dir = os.getcwd() # Add the parent directory to Python path parent_dir = os.path.dirname(current_dir) if parent_dir not in sys.path: sys.path.append(parent_dir) # Import PyHydroGeophysX modules from PyHydroGeophysX.inversion.time_lapse import TimeLapseERTInversion from PyHydroGeophysX.inversion.windowed import WindowedTimeLapseERTInversion .. GENERATED FROM PYTHON SOURCE LINES 46-68 .. code-block:: default data_dir = os.path.join(current_dir, "results","TL_measurements","appres") #"C:/Users/HChen8/Documents/GitHub/PyHydroGeophysX/examples/results/TL_measurements/appres" # List of ERT data files testing monthly time-lapse inversion ert_files = [ "synthetic_data30.dat", "synthetic_data60.dat", "synthetic_data90.dat", "synthetic_data120.dat", "synthetic_data150.dat", "synthetic_data180.dat", "synthetic_data210.dat", "synthetic_data240.dat", "synthetic_data270.dat", "synthetic_data300.dat", "synthetic_data330.dat", "synthetic_data360.dat", ] .. GENERATED FROM PYTHON SOURCE LINES 69-70 ## 1.Full L2 Time-Lapse Inversion .. GENERATED FROM PYTHON SOURCE LINES 72-73 #################### FUll Time-Lapse Inversion ##################### .. GENERATED FROM PYTHON SOURCE LINES 73-113 .. code-block:: default # Full paths to data files data_files = [os.path.join(data_dir, f) for f in ert_files] # Measurement times (can be timestamps or any sequential numbers representing time) measurement_times = [1, 2, 3, 4, 5, 6, 7 ,8, 9, 10, 11, 12] # Adjust based on your actual acquisition times # Create a mesh for the inversion (or load an existing one) data = ert.load(data_files[0]) ert_manager = ert.ERTManager(data) mesh = ert_manager.createMesh(data=data, quality=34) # Set up inversion parameters inversion_params = { "lambda_val": 50.0, # Regularization parameter "alpha": 10.0, # Temporal regularization parameter "decay_rate": 0.0, # Temporal decay rate "method": "cgls", # Solver method ('cgls', 'lsqr', etc.) "model_constraints": (0.001, 1e4), # Min/max resistivity values (ohm-m) "max_iterations": 15, # Maximum iterations "absoluteUError": 0.0, # Absolute data error (V) "relativeError": 0.05, # Relative data error (5%) "lambda_rate": 1.0, # Lambda reduction rate "lambda_min": 1.0, # Minimum lambda value "inversion_type": "L2" # 'L1', 'L2', or 'L1L2' } # Create the time-lapse inversion object inversion = TimeLapseERTInversion( data_files=data_files, measurement_times=measurement_times, mesh=mesh, **inversion_params ) # Run the inversion print("Starting time-lapse inversion...") result = inversion.run() print("Inversion complete!") .. GENERATED FROM PYTHON SOURCE LINES 114-126 .. code-block:: default from palettable.lightbartlein.diverging import BlueDarkRed18_18 fixed_cmap = BlueDarkRed18_18.mpl_colormap fig = plt.figure(figsize=[16,6]) # True resistivity model for i in range(12): ax = fig.add_subplot(3,4,i+1) ax, cbar = pg.show(result.mesh,result.final_models[:,i],pad=0.3,orientation="vertical",cMap=fixed_cmap,cMin= 100,cMax = 3000 , ylabel="Elevation (m)",label=' Resistivity ($\Omega$ m)',ax=ax,logScale=False,coverage=result.all_coverage[i]>-1) cbar.remove() .. GENERATED FROM PYTHON SOURCE LINES 127-140 Full Time-Lapse ERT Inversion Results ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The time-lapse sequence reveals the temporal evolution of subsurface water content over 12 months. Notice how the resistivity distribution changes systematically, with high resistivity zones (dry areas) gradually becoming more conductive (wet) as water infiltrates through the soil layers. The temporal regularization ensures smooth transitions between consecutive time steps while capturing realistic subsurface changes. .. image:: /auto_examples/images/Ex_TL_inversion_fig_01.png :align: center :width: 900px .. GENERATED FROM PYTHON SOURCE LINES 143-144 ## 2. Window L2 time-lapse inversion .. GENERATED FROM PYTHON SOURCE LINES 146-147 Measurement times (can be timestamps or any sequential numbers representing time) .. GENERATED FROM PYTHON SOURCE LINES 147-188 .. code-block:: default measurement_times = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12] # Adjust based on your actual acquisition times # Create a mesh for the inversion (or load an existing one) data = ert.load(os.path.join(data_dir, ert_files[0])) ert_manager = ert.ERTManager(data) mesh = ert_manager.createMesh(data=data, quality=34) # Set up inversion parameters inversion_params = { "lambda_val": 10.0, # Regularization parameter "alpha": 10.0, # Temporal regularization parameter "decay_rate": 0.0, # Temporal decay rate "method": "cgls", # Solver method ('cgls', 'lsqr', etc.) "model_constraints": (0.001, 1e4), # Min/max resistivity values (ohm-m) "max_iterations": 15, # Maximum iterations "absoluteUError": 0.0, # Absolute data error (V) "relativeError": 0.05, # Relative data error (5%) "lambda_rate": 1.0, # Lambda reduction rate "lambda_min": 1.0, # Minimum lambda value "inversion_type": "L2" # 'L1', 'L2', or 'L1L2' } # Define the window size (number of timesteps to process together) window_size = 3 # A window size of 3 means each window includes 3 consecutive measurements # Create the windowed time-lapse inversion object inversion = WindowedTimeLapseERTInversion( data_dir=data_dir, # Directory containing ERT data files ert_files=ert_files, # List of ERT data filenames measurement_times=measurement_times, # List of measurement times window_size=window_size, # Size of sliding window mesh=mesh, # Mesh for inversion **inversion_params # Pass the same inversion parameters ) # Run the inversion, optionally in parallel print("Starting windowed time-lapse inversion...") result = inversion.run() print("Inversion complete!") .. GENERATED FROM PYTHON SOURCE LINES 189-192 .. code-block:: default result.final_models = np.array(result.final_models) result.final_models.shape .. GENERATED FROM PYTHON SOURCE LINES 193-254 .. code-block:: default from palettable.lightbartlein.diverging import BlueDarkRed18_18 import matplotlib.pyplot as plt import numpy as np import matplotlib.pylab as pylab params = {'legend.fontsize': 13, #'figure.figsize': (15, 5), 'axes.labelsize': 13, 'axes.titlesize':13, 'xtick.labelsize':13, 'ytick.labelsize':13} pylab.rcParams.update(params) plt.rcParams["font.family"] = "Arial" fixed_cmap = BlueDarkRed18_18.mpl_colormap fig = plt.figure(figsize=[16, 6]) # Use tight_layout with adjusted parameters to reduce space plt.subplots_adjust(wspace=0.05, hspace=0.05) # True resistivity model for i in range(12): row, col = i // 4, i % 4 ax = fig.add_subplot(3, 4, i+1) # Add common ylabel only to leftmost panels ylabel = "Elevation (m)" if col == 0 else None # Add resistivity label only to the middle-right panel (row 1, col 3) resistivity_label = ' Resistivity ($\Omega$ m)' if (i == 7) else None # Only show axis ticks on leftmost and bottom panels if col != 0: ax.set_yticks([]) if row != 2: # Not bottom row ax.set_xticks([]) else: # Add "distance (m)" label to bottom row panels ax.set_xlabel("Distance (m)") # Create the plot ax, cbar = pg.show(result.mesh, result.final_models[:,i], pad=0.3, orientation="vertical", cMap=fixed_cmap, cMin=100, cMax=3000, ylabel=ylabel, label=resistivity_label, ax=ax, logScale=False, coverage=result.all_coverage[i]>-1.2) # Only keep colorbar for the middle-right panel (row 1, col 3) # This corresponds to panel index 7 in a 0-based indexing system if i != 7: # Keep only the colorbar for panel 7 cbar.remove() plt.tight_layout() .. GENERATED FROM PYTHON SOURCE LINES 255-268 Windowed Time-Lapse Inversion Results ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The windowed approach processes consecutive time steps in overlapping groups, which is computationally more efficient for large datasets while maintaining temporal coherence. This method shows similar overall patterns to the full time-lapse inversion but with reduced computational cost. The windowing approach is particularly valuable when processing years of monitoring data or when computational resources are limited. .. image:: /auto_examples/images/Ex_TL_inversion_fig_02.png :align: center :width: 900px .. GENERATED FROM PYTHON SOURCE LINES 270-271 ## 3. Full L1 Time-lapse Inversion .. GENERATED FROM PYTHON SOURCE LINES 276-279 .. code-block:: default ax, cbar = pg.show(result.mesh,result.final_models[:,i],pad=0.3,orientation="vertical",cMap=fixed_cmap,cMin= 100,cMax = 3000 , ylabel="Elevation (m)",label=' Resistivity ($\Omega$ m)',logScale=False,coverage=result.all_coverage[i]>-1) .. GENERATED FROM PYTHON SOURCE LINES 283-284 Full paths to data files .. GENERATED FROM PYTHON SOURCE LINES 284-322 .. code-block:: default data_files = [os.path.join(data_dir, f) for f in ert_files] # Measurement times (can be timestamps or any sequential numbers representing time) measurement_times = [1, 2, 3, 4, 5, 6, 7 ,8, 9, 10, 11, 12] # Adjust based on your actual acquisition times # Create a mesh for the inversion (or load an existing one) data = ert.load(data_files[0]) ert_manager = ert.ERTManager(data) mesh = ert_manager.createMesh(data=data, quality=34) # Set up inversion parameters inversion_params = { "lambda_val": 50.0, # Regularization parameter "alpha": 10.0, # Temporal regularization parameter "decay_rate": 0.0, # Temporal decay rate "method": "cgls", # Solver method ('cgls', 'lsqr', etc.) "model_constraints": (0.001, 1e4), # Min/max resistivity values (ohm-m) "max_iterations": 15, # Maximum iterations "absoluteUError": 0.0, # Absolute data error (V) "relativeError": 0.05, # Relative data error (5%) "lambda_rate": 1.0, # Lambda reduction rate "lambda_min": 1.0, # Minimum lambda value "inversion_type": "L1" # 'L1', 'L2', or 'L1L2' } # Create the time-lapse inversion object inversion = TimeLapseERTInversion( data_files=data_files, measurement_times=measurement_times, mesh=mesh, **inversion_params ) # Run the inversion print("Starting time-lapse inversion...") result_L1 = inversion.run() print("Inversion complete!") .. GENERATED FROM PYTHON SOURCE LINES 323-335 .. code-block:: default from palettable.lightbartlein.diverging import BlueDarkRed18_18 fixed_cmap = BlueDarkRed18_18.mpl_colormap fig = plt.figure(figsize=[16,6]) # True resistivity model for i in range(12): ax = fig.add_subplot(3,4,i+1) ax, cbar = pg.show(result_L1.mesh,result_L1.final_models[:,i],pad=0.3,orientation="vertical",cMap=fixed_cmap,cMin= 100,cMax = 3000 , ylabel="Elevation (m)",label=' Resistivity ($\Omega$ m)',ax=ax,logScale=False,coverage=result.all_coverage[i]>-1) cbar.remove() .. GENERATED FROM PYTHON SOURCE LINES 336-349 L1-Norm Regularized Inversion Results ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The L1-norm regularization produces sharper boundaries compared to L2 regularization, which is beneficial for detecting distinct geological layers or sharp interfaces. Notice how the resistivity transitions are more abrupt in this inversion, making it easier to identify distinct zones of different water content. This approach is particularly useful when the subsurface is expected to have layered structures rather than gradual transitions. .. image:: /auto_examples/images/Ex_TL_inversion_fig_03.png :align: center :width: 900px .. GENERATED FROM PYTHON SOURCE LINES 352-376 Summary and Recommendations ~~~~~~~~~~~~~~~~~~~~~~~~~~ This example demonstrated three approaches to time-lapse ERT inversion: **Key Findings:** * **Temporal regularization** is essential for realistic time-lapse results * **Windowed inversion** provides computational efficiency with minimal quality loss * **L1 regularization** enhances boundary detection in layered media * **Parameter tuning** (λ, α) significantly affects result quality **Recommendations:** * Use **full L2** for high-quality results with moderate datasets * Apply **windowed L2** for large datasets or real-time processing * Choose **L1 regularization** when sharp interfaces are expected * Always validate results against known geological information **Next Steps:** * Combine with seismic constraints (see structure-constrained examples) * Apply uncertainty quantification (see Monte Carlo examples) * Integrate with hydrological models for enhanced interpretation .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.000 seconds) .. _sphx_glr_download_auto_examples_Ex_TL_inversion.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: Ex_TL_inversion.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: Ex_TL_inversion.ipynb `