gimli-org / gimli

Geophysical Inversion and Modeling Library :earth_africa:
https://www.pygimli.org
Other
344 stars 131 forks source link

Separate IP inversion #652

Open jcmefra opened 3 months ago

jcmefra commented 3 months ago

Hello, hope you are doing good. I want to know if there is a way to do a separate IP inversion instead of a joint ip-ert inversion with the ERTIPmanager(). Thank you so much in advance.

Also, if it's possible, how do I access to the individual error metrics of each inversion in the joint inversion? I'm particularly interested in the ABS error % to compare with Res2Dinv. I know the inversion method is different but that's why I want to use both to see how much my results differ.

image image

halbmy commented 2 months ago

In contrast to frequency domain, where one could invert for complex resistivity, there is no "joint inversion" in time-domain IP. Instead, first a resistivity inversion is done (calling the normal ERTManager inversion) and then a IP inversion that can also be called by ERTIPManager.invertTDIP(). You should see things like RMS in the output but I just see that the inversion object is not kept for further analysis. I will change this.

jcmefra commented 2 months ago

So it's okay to save both inversion objects by doing inv_ert, inv_ip = mgr2.invert(args) ? I have access to the inversion verbose, but I don't know how to access to certain metrics specifically. Thank you.

Edit: does the inversion work for msec as well? I just changed the label, but I see that the ip domain is "mV/V".

Also, when I do this to invert only IP:

inv_ip = mgr.invertTDIP(secNodes=1, paraMaxCellSize=10.0,
                        zWeight=0.2, verbose=1)

I get this error:

RuntimeError: RegionManager knows no mesh.

But when I do the regular mgr.invert(secNodes=1, paraMaxCellSize=10.0, zWeight=0.2, verbose=1, it just does both inversions without defining any mesh.

halbmy commented 2 months ago
  1. No, the inversion object is not returned. We have to change this in the source code (to keep it as a member of the manager).
  2. Formally it is working and the results will be ok as a first-order approximation but the inversion is based on a instantaneous and not integrated chargeability.
  3. Before calling mgr.invertTDIP() you need to do a resistivity inversion through mgr.invert() (which will already do an IP inversion).

Thank you for pointing out the shortcomings in the code.

jcmefra commented 2 months ago

Alright, thank you so much. I've got another question, is there a way to remove those data points where data/response relation is very high or very low and therefore lead to a higher error%? Probably doing a scatterplot for data vs response and delete those exceeding certain thereshold. so I can reduce the chi2 directly with the inversion data.

I'm not pretty sure how to actually handle errors post-inversion. Thanks in advance.

halbmy commented 1 week ago

Yes, of course. The data vector is in data['rhoa'] and the response vector in mgr.inv.response. You can compute the rms or error-weighted misfit and look at the histogram or the misfit distribution and then remove points

misfit = (data['rhoa'] - mgr.inv.response) / data['err']
data.remove(np.abs(misfit) > 10)