gimli-org / gimli

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

Questions - Geostatistical regularisation #600

Closed Bhergam closed 7 months ago

Bhergam commented 8 months ago

Questions regarding Geostat regularisation

Hello,I would like to know if someone could help me with Geostatistical regularization for 3D ERT, I’m currently trying to implement a simple simulation with a single prior data in a rectangular shape. I’m using a geostatistical constraint matrix to see its effects. I have encountered some things, I would like to disscus it :

Firstly while calling the eigenavlues ew and eigenvectors EV from a geostatistical matrix using the methods Cm05.ew and Cm05.EV respectively, the matrices derived are D and Q respectively according to the method numpy.linalg.eigh such as :

Cm = QDQ_t

Where Q_t is the inverse matrix of Q because this writting implies that Q belongs to the othogonal group.For a number of deined number of cells of 144 , the absolute value of the det( Q ) is 1 but for a higher number of cell ( 3500 ), the value of the determinant is -0.5 I believe that its due to some numerical approximation because I can see as well that some values in the matrix product of Q bu Q_transposed that should be 0 are closed to 10 **-16, do you have any ideas regarding this observation ? I tried to change the lower threeshold in the source code of Cm05 which I believe is the lower bound of the eigenvalues tolerable, I tried with –np.inf but the result seems the same so I don’t think it comes from this method but rather from some numerical approximation.

Depth Weighting

I have another question regarding the use of inv.setConstraintWeights(cWeight).. Halbmy metionned in another thread that « the weighting vector (containing the weight of the individual lines of the matrix)… defaults to one » if not specified ( link to the message ) but I would like to know how it works exactly.

Starting with its length, it is thus inuitive that it should be equal to the number of rows of the Cmatrix but one can compute a different weight vector with a different length and obtain a non absurd result. For the example, on the example given in 4 example gravimetry and magnetics , we can see that a depth weighting is perfom with the use of the altitude of each facet of cells that isn’t a bondary, the length of the array containing the weighting factors is related to the number of cells, but it is not neccesary the number of rows of the Cmatrix.

I also thought that there was a use of pg.matrix.MultLeftMatrix of Cmatrix by w_z but after some tests, I don’t think so. My skills in C++ aren’t sharp enough, I don’t understand all the code source, I found that in mesh.py in the function def_createParameterConstraintsLines that cWeights = np.ones(cMat.rows) if cWeights is None but I believe that inv.run doesn’t use such function directly so I don’t understand much about how the use of the weight matrix.

I hope I made myself clear enough, Kind regards.

halbmy commented 7 months ago

To the eigenvalue decomposition: Of course also eigenvalue decomposition is prone to numerical accuracy, which is in the order of 1e-16 for double values. So I don't see anything wrong in that and it does not practicably affect the inversion result, just like it doesn't when disregarding low eigenvalues.

While for a first-order derivative operator the weighting vector has the length of inner boundaries, it needs to be the number of cells for second-order derivative, damping and also geostatistic constraints. However, I am not sure whether it is wise to destroy this symmetric operator by some weighting. Nevertheless you can try setting it like in the magnetics example.