Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

3D gravity inversion of prism layers #322

Open
mdtanker opened this issue Jun 7, 2022 · 12 comments
Open

3D gravity inversion of prism layers #322

mdtanker opened this issue Jun 7, 2022 · 12 comments
Labels
enhancement Idea or request for a new feature

Comments

@mdtanker
Copy link
Member

mdtanker commented Jun 7, 2022

Description of the desired feature:

I'm working on a 3D gravity inversion for my Ph.D. to model the bathymetry beneath a float ice shelf. I'm using a vertical-prisms approach, and I'm mostly interested in a geometry inversion, where each prism's density remains unchanged, but the prisms' tops or bottoms are updated to minimize the observed gravity - forward gravity misfit. While this inversion is in the early stages, I'm already using several harmonica tools so I think it would be great to include a fully supported inversion module in harmonica.

Below is a description of my specific use case, a list of what I already have coded, and ideas for additional features to add.

Current state of inversion

Model setup

  • input an arbitrary number of layer topographies (i.e. ice surface, ice base, bathymetry, basement, Moho)
  • create a layer of vertical prisms between each set of grids
    • harmonica.prism_layer()
      • example: for water layer: surface = ice base, reference = bathymetry
    • assign a density to each prism (currently using constant density within each layer)
      • prism_layer(... properties={'density':water_density)
    • allow different cell sizes between layers
      • if lower grid cell size doesn't match upper grid's, use pygmt.grdtrack() to sample the lower grids values at the grid's prism locations

  • calculate gravity misfit
    • add forward gravities for all layers, and subtract from observed gravity

run geometry inversion

  • designate an active_layer ex: active_layer='bathymetry_prisms'
  • each iteration of the inversion will yield a surface_correction
    • set a max surface change per iteration
    • applied to the active_layer tops, as well as the bottoms of the layer above active_layer
    • recalculate the forward gravity of these 2 updated layers.
      • recalculte the misfit with these updated forward gravities before starting the 2nd iteration
  • currently this inversion works by finding the least-squares solution to the matrix equation Ax=b, where A is the Jacobian matrix of the vertical derivative of gravity acceleration, and b is the initial misfit between observed and forward gravity.
    • the least squares solution is found with scipy.sparse.linalg.lsqr()
    • the vertical derivative is approximated with the Hammer approximation of prisms, using annuluses as opposed to prisms.

This is just the method I already had implemented, suggest from my advisor, but I'm happy to alter this and would appreciate suggestions.

Features to add:

  • enforce constraint point (points of known bathymetry) by multiplying surface_correction by a constraints grid
    * 0's at constraints, linearly increasing to 1's at a specified distance from nearest constraints
    * scipy.spatial.cKDTree.query() to get grid of minimum distances to constraints
    * gmt grdmath LE and GE to set constraints to 0 and other points to 1
    * gmt grdblend to merge the clipped grids
    * and rioxarray.interpolate_na() to linearly interpolate between 0 and 1

  • allow irregular grid spacing within each layer (i.e. grid spacing increase for a buffer zone, outside the main area of interest)
  • allow depth-dependent density within each prism. Only for the forward modelling
  • density inversion
    • similar to the above geometry inversion, but updated each prism's density to minimize the misfit.
    • this is useful for the initial setup of the geometry inversion, allow the regional field to be accounted for in a layers density, leaving only the residual field left to invert on.

Discussion

  • Do we want to refer to this type of inversion as a geometry inversion, structural inversion, boundary inversion or something else? There doesn't seem to be much of a consensus on what you call this.
  • I'm pretty new to inversion theory, so if anyone has recommendations on where to learn this material please let me know!
  • I think some of the packages I have linked above could be replaced with Fatianado packages to keep dependencies minimal.
  • All the inputs in my inversion are stored in a nested dictionary layers, with a dictionary for each layer, which includes a .nc filename, a resolution, constant density values, and a dataframe and dataset for each layer. This makes it easy to apply functions to an arbitrary number of input layers (example below). If there's a better method for this let me know.
layers_list =[ 'ice', 'water', 'bathymetry', 'basement',]
spacing_list = [10e3, 10e3, 10e3, 10e3,]
rho_list = [ 920, 1030, 2600, 2800,]
fname_list=[
            '../inversion_layers/BedMachine_surface_filled_inv.nc',
            '../inversion_layers/BedMachine_icebase_filled_inv.nc',
            ../inversion_layers/BedMachine_bed_inv.nc',
            '../inversion_layers/ROSETTA_basement_BedMachine_bed_inv.nc',
            ]

# create nested dictionary of layer properties 
layers = {j:{'spacing':spacing_list[i], 
            'fname':fname_list[i], 
            'rho':rho_list[i]} for i, j in enumerate(layers_list)}

for i, j in enumerate(layers_list):
    # load each .nc as a xarray dataset
    layers[j]['grid']=xr.load_dataset(layers[j]['fname'])
    # create prism layer between each set of input grids
    layers[j]['prisms']=hm.prism_layer(
        coordinates=(list(layers[j]['grid'].x), list(layers[j]['grid'].y)),   
        surface=layers[j]['grid'].z, 
        reference=layers[reversed_layers_list[i-1]]['grid'].z,
        properties={'density':layers[j]['grid'].density})

# calculate forward gravity for each prism layer
for k, v in layers.items():
    df_grav[f'{k}_forward_grav'] = v['prisms'].prism_layer.gravity(
        coordinates = (df_grav.x, df_grav.y, df_grav.z),
        field = 'g_z', progressbar=True)
  • to remove edge effects, I'm using a buffer, currently 200km in each direction, as shown in the below figures by the black square. You can see the edge effects in the second to last panel, which is the lowest layer of prisms. I am doing all calculations on the full region, but only showing results within the buffer zone. Does this seem like a reasonable approach? It adds a lot of prisms, but this could be mitigated by using discretize.TreeMesh(). for the buffer region. I'm unsure if I should be calculating misfit for the whole region or just the inside. Also, I'm not sure if I should invert the prisms just within the inside, or for the whole region.

Are you willing to help implement and maintain this feature?

Yes I'm happy and excited to implement this, but I'm in the last year of my Ph.D. so the outcome of this inversion is more pressing than the implementation of the inversion into Harmonica. Hopefully they can happen simultaneously. @santisoler has suggested looking at both the old Fatiando module fatiando.inversion and some code related to Uieda and Barbosa 2016 for inspiration on how to implement this, which I will do.

@mdtanker mdtanker added the enhancement Idea or request for a new feature label Jun 7, 2022
@LL-Geo
Copy link
Member

LL-Geo commented Jun 22, 2022

Just find this paper, that might give an idea of how to include different density contrast (maybe different geology) into the inversion. Haas et al 2020. https://academic.oup.com/gji/article/221/3/1896/5808816?login=true

@craigmillernz
Copy link

This would be very useful feature. There is some overlap with SimPEG here. Perhaps SimPEG inversion engine could be used and linked to Harmonica for forward model with prisms. I think @santisoler is looking at this in SimPEG at the moment to rework SimPEG gravity forward modeling. It would save reinventing a lot of inverse code that is already well used and supported.

VPMG commercial code that does this style calls this a geometry inversion, VPMG also has the option to alternate 1 geometry iteration and one property (i.e density) iteration then back to geometry iteration etc to solve for both geometry and density.

@santisoler
Copy link
Member

Hi @craigmillernz. Thanks for joining the conversation. I agree that we should make efforts to avoid overlaps with SimPEG. And you're right, I'm working on building some common ground for prism forward modelling between SimPEG and Harmonica. I'm looking forward to getting more awesome results.

Nevertheless, I like the way that a geometry inversion could be built on top of Harmonica. SimPEG is mainly developed for inverting the physical properties of an already discretized space. I'm guessing that trying to implement something like the geometry inversion carried out by Uieda et al (2017) on top of the SimPEG framework could be challenging. But again, I don't mean to write a complete separate inversion framework in Harmonica just for that. Maybe we could start with specific inversion functions, like invert_moho for inverting the Moho or invert_basin for a sedimentary basin, more in sync with a use-case oriented development. What do you think?

@mdtanker
Copy link
Member Author

Thanks for pointing that out @LL-Geo, I'll look through it.

@mdtanker
Copy link
Member Author

Santi I like your idea of having a few specific inversions for common tasks, but it might end up with be a lot of duplicated code. What do you think are the benefits over having a generic invert_prism_geometry function, with specific gallery examples for applying it to common use cases?

@craigmillernz
Copy link

I tend to agree with Matt that a more generic function with use cases might be better. I was really just wondering how easy it would be to couple the prisms discretization onto the simpeg inversion framework, You'll know best about that! If they aren't really compatible, then a separate inversion within Harmonica would be appropriate.

@craigmillernz
Copy link

Matt you might find this useful for general inversion theory - from UBC. This used to be a standalone website, but it looks like you now have to download. Follow the instructions on the page.
https://gif.eos.ubc.ca/IAG

@mdtanker
Copy link
Member Author

mdtanker commented Dec 9, 2022

Hi everyone, sorry for the lack of updates here! I've been struggling to implement regularization for a while and was wondering if anyone had ideas to help. I've added some smoothness regularization within the least-squares solver with a damping term (with verde.base.least_squares), which helps stabilize the inversion. But, I'd like to add some hard constraints into the regularization as well, which seems to be referred to as "smallness" regularization (as discussed here).

There are a series of seismic constraints on the bathymetry depths which I would like to be unchanged (within error margins) after the inversion. Ideally, at each iteration, the correction value at these points would be 0, and the smoothness term will help ensure these constraint points don't turn into pedestals, as the unconstrained bathymetry surrounding them changes.

Does anyone have ideas for how to implement this? Would it be within the least squares solver?

I've tried to follow along with similar implementations in PyGIMLi, SimPEG, and the legacy Fatiando, but I've either only seen smoothness regularization, or I'm unable to understand how it's coded.

Any help or suggestions would be appreciated!

@santisoler
Copy link
Member

Hi @mdtanker. I'm glad you are getting into the rabbit hole of inversion, it's a fun one haha!

The good news here is that you are already using some simple kind of smallness regularization. The verde.base.least_squares function relies on the sklearn.linear_model.Ridge class that minimizes the objective function:

$$ \phi(\mathbf{m}) = \lVert \mathbf{d}_\text{obs} - \mathbf{G} \mathbf{m} \rVert^2 + \alpha \lVert \mathbf{m} \rVert^2 $$

The second term that involves the alpha (damping parameter) and the l2 norm of the parameters is a Tikhonov zeroth order regularization. SimPEG implements the smallness by defining this regularization term as the l2 norm of the difference between the model parameters and the reference model. In this reference model you can include you constrains:

$$ \phi(\mathbf{m}) = \lVert \mathbf{d}\text{obs} - \mathbf{G} \mathbf{m} \rVert^2 + \alpha \lVert \mathbf{m} - \mathbf{m}\text{ref} \rVert^2 $$

Here's a nice chapter to get into these topics: https://pubs.geoscienceworld.org/books/book/2111/chapter-abstract/114879875/Inversion-for-Applied-Geophysics-A-Tutorial?redirectedFrom=fulltext

The problem is that you won't be able to use the Ridge class for defining a reference model and you'll probably have to write you own implementation of the regularization.

There are also some lecture notes from @leouieda and @birocoles. They have really nice stuff there, particularly if you need to implement this stuff: https://www.semanticscholar.org/paper/T%C3%B3picos-de-invers%C3%A3o-em-geof%C3%ADsica-OliveiraVanderlei-Leonardo/fd651157443f37ad603057a10d683322ff123263

@WYJLCYWHZ
Copy link

Hi @mdtanker ,I came across your project on GitHub and I'm truly impressed by the remarkable feature you've been working on. Your research on 3D gravity inversion for modeling the bathymetry beneath a float ice shelf sounds fascinating. I'm particularly intrigued by your implementation of the vertical-prisms approach and the focus on geometry inversion to minimize the observed gravity misfit.

I'm genuinely interested in learning from your implementation and gaining insights into the techniques you've employed. It would greatly benefit my own learning and allow me to explore potential areas for improvement.I would be grateful if you could kindly share the code you've already developed. Thank you (^▽^)

@mdtanker
Copy link
Member Author

mdtanker commented Jul 4, 2023

Hi @WYJLCYWHZ. apologies for the slow reply. Thank you for your kind words on my work. The inversion code is included in my repository RIS_gravity_inversion, but it is written specifically for my use-case and is not very well documented or tested yet. In the coming months, I will work on making this inversion code more accessible and understandable for other users. I'll try and post updates in this thread or in a Fatiando discussion. Let me know if you have specific questions.

@WYJLCYWHZ
Copy link

Hi @WYJLCYWHZ. apologies for the slow reply. Thank you for your kind words on my work. The inversion code is included in my repository RIS_gravity_inversion, but it is written specifically for my use-case and is not very well documented or tested yet. In the coming months, I will work on making this inversion code more accessible and understandable for other users. I'll try and post updates in this thread or in a Fatiando discussion. Let me know if you have specific questions.

Thank you very much (^▽^)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
enhancement Idea or request for a new feature
Projects
None yet
Development

No branches or pull requests

5 participants