Questions - Geostatistical regularisation

[Bhergam]

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]

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.

[makeabhishek]

Thanks for the explanation @halbmy. I would like to know a bit more about smoothening and damping. As I followed the example on regularization, which states that cType = 1: # first order smoothness; cType = 2: # 2nd order smoothing; cType = 10: # mix of 1st order smoothing and damping.

As cType =0 defines the damping in (traveltime) inversion. What exactly defines the damping, what is the formulation for that?

What is the starting model mRef used by inverse model and how does inversion choose mRef, if we don't specify it explicitly?

[halbmy]

For damping (cType=0) the regularization matrix is the identity matrix, so $|m-mRef|$ is minimized.

mRef is by default zero but you can use the starting model through isReference=True.

The starting model for traveltime is

• a gradient model derived from vTop and vBottom (useGradient=True)
• a constant slowness derived from mean apparent velocity (useGradient=False)

[makeabhishek]

Thanks @halbmy.
(1) If its a unity matric for cType=0, what does damping means here? what damping parameter is assigned?

(2) What if we don't have the vTop and vBottom defined (in case of circular geometry) for example, (lam=12, useGradient=1, zWeight=15.0). How does the initial model defined in default case?

(3) Furthermore, how does cType=10 combine 1st order smoothing and damping? Is there any reference for this? What is the mathemtical formulation?

[halbmy]

1. Damping means minimizing the norm of the model (or its difference to the reference)
2. If vTop and vBottom are equal, the model is homogeneous. But then you can also use startModel. By default, the median apparent velocity is used.
3. cType=10 adds cType=0 and cType=1, i.e. adds the norm of the model and the model roughness (1st order). See regularization tutorial https://www.pygimli.org/_tutorials_auto/3_inversion/plot_5_Regularization.html