Forward and inverse modelling

In this tutorial, we will see how to use ResIPy API to do forward and inverse modelling of synthetic problems.

In [1]:
from resipy import R2
k = R2(typ='R2')
C:\Program Files (x86)\WPy64-3770\python-3.7.7.amd64\lib\site-packages\resipy\meshTools.py:50: UserWarning: pyvista not installed, 3D meshing viewing options will be limited
  warnings.warn('pyvista not installed, 3D meshing viewing options will be limited')
API path =  C:\Program Files (x86)\WPy64-3770\python-3.7.7.amd64\lib\site-packages\resipy
ResIPy version =  2.2.2
cR2.exe found and up to date.
R3t.exe found and up to date.
cR3t.exe found and up to date.
Working directory is: C:\Program Files (x86)\WPy64-3770\python-3.7.7.amd64\lib\site-packages\resipy
clearing dirname
In [2]:
import numpy as np # numpy for electrode generation

First we need to design some surface electrodes using numpy

In [3]:
elec = np.zeros((25,3))
elec[:,0] = np.arange(0, 25*5.0, 5.0) # with 25 electrodes at 5 m spacing 
k.setElec(elec)
print(k.elec)
        x    y    z  remote  buried label
0     0.0  0.0  0.0   False   False     1
1     5.0  0.0  0.0   False   False     2
2    10.0  0.0  0.0   False   False     3
3    15.0  0.0  0.0   False   False     4
4    20.0  0.0  0.0   False   False     5
5    25.0  0.0  0.0   False   False     6
6    30.0  0.0  0.0   False   False     7
7    35.0  0.0  0.0   False   False     8
8    40.0  0.0  0.0   False   False     9
9    45.0  0.0  0.0   False   False    10
10   50.0  0.0  0.0   False   False    11
11   55.0  0.0  0.0   False   False    12
12   60.0  0.0  0.0   False   False    13
13   65.0  0.0  0.0   False   False    14
14   70.0  0.0  0.0   False   False    15
15   75.0  0.0  0.0   False   False    16
16   80.0  0.0  0.0   False   False    17
17   85.0  0.0  0.0   False   False    18
18   90.0  0.0  0.0   False   False    19
19   95.0  0.0  0.0   False   False    20
20  100.0  0.0  0.0   False   False    21
21  105.0  0.0  0.0   False   False    22
22  110.0  0.0  0.0   False   False    23
23  115.0  0.0  0.0   False   False    24
24  120.0  0.0  0.0   False   False    25

Now we need to assign objects we wish to add and then create a mesh

In [19]:
geom_input = {'polygon1':[[30,45,45,30], [-5, -5, -30, -30]],
             'polygon2':[[75,90,90,75],[-5,-5,-30,-30]]}
k.createMesh(typ='trian', res0=200, geom_input=geom_input) # let's create the mesh based on these electrodes position
k.showMesh(xlim=[0,120], zlim=[-40,0])
Creating triangular mesh...Reading mesh.msh
Gmsh version == 3.x
reading node coordinates...
Determining element type...Triangle
Reading connection matrix...
ignoring 0 elements in the mesh file, as they are not required for R2/R3t
Finished reading .msh file
ResIPy Estimated RAM usage = 0.003983 Gb
done

Based on this mesh, we can define resistities for the regions

In [5]:
k.setStartingRes(regionValues={1:50, 2:500, 3:200}) # this assign a resistivity value (ohm.m) for each region defined.
k.showMesh(attr='res0')
regionValues: {1: 50, 2: 500, 3: 200}

We then need to define the sequence that we will use. We can easily create a dipole-dipole sequence using R2.createSequence() or import one using R2.importSequence().

In [6]:
k.createSequence([('dpdp1', 1, 10)]) # create a dipole-dipole of diple spacing of 1 (=skip 0) with 10 levels
print(k.sequence) # the sequence is stored inside the R2 object
175 quadrupoles generated.
[[ 1  2  3  4]
 [ 2  3  4  5]
 [ 3  4  5  6]
 [ 4  5  6  7]
 [ 5  6  7  8]
 [ 6  7  8  9]
 [ 7  8  9 10]
 [ 8  9 10 11]
 [ 9 10 11 12]
 [10 11 12 13]
 [11 12 13 14]
 [12 13 14 15]
 [13 14 15 16]
 [14 15 16 17]
 [15 16 17 18]
 [16 17 18 19]
 [17 18 19 20]
 [18 19 20 21]
 [19 20 21 22]
 [20 21 22 23]
 [21 22 23 24]
 [22 23 24 25]
 [ 1  2  4  5]
 [ 2  3  5  6]
 [ 3  4  6  7]
 [ 4  5  7  8]
 [ 5  6  8  9]
 [ 6  7  9 10]
 [ 7  8 10 11]
 [ 8  9 11 12]
 [ 9 10 12 13]
 [10 11 13 14]
 [11 12 14 15]
 [12 13 15 16]
 [13 14 16 17]
 [14 15 17 18]
 [15 16 18 19]
 [16 17 19 20]
 [17 18 20 21]
 [18 19 21 22]
 [19 20 22 23]
 [20 21 23 24]
 [21 22 24 25]
 [ 1  2  5  6]
 [ 2  3  6  7]
 [ 3  4  7  8]
 [ 4  5  8  9]
 [ 5  6  9 10]
 [ 6  7 10 11]
 [ 7  8 11 12]
 [ 8  9 12 13]
 [ 9 10 13 14]
 [10 11 14 15]
 [11 12 15 16]
 [12 13 16 17]
 [13 14 17 18]
 [14 15 18 19]
 [15 16 19 20]
 [16 17 20 21]
 [17 18 21 22]
 [18 19 22 23]
 [19 20 23 24]
 [20 21 24 25]
 [ 1  2  6  7]
 [ 2  3  7  8]
 [ 3  4  8  9]
 [ 4  5  9 10]
 [ 5  6 10 11]
 [ 6  7 11 12]
 [ 7  8 12 13]
 [ 8  9 13 14]
 [ 9 10 14 15]
 [10 11 15 16]
 [11 12 16 17]
 [12 13 17 18]
 [13 14 18 19]
 [14 15 19 20]
 [15 16 20 21]
 [16 17 21 22]
 [17 18 22 23]
 [18 19 23 24]
 [19 20 24 25]
 [ 1  2  7  8]
 [ 2  3  8  9]
 [ 3  4  9 10]
 [ 4  5 10 11]
 [ 5  6 11 12]
 [ 6  7 12 13]
 [ 7  8 13 14]
 [ 8  9 14 15]
 [ 9 10 15 16]
 [10 11 16 17]
 [11 12 17 18]
 [12 13 18 19]
 [13 14 19 20]
 [14 15 20 21]
 [15 16 21 22]
 [16 17 22 23]
 [17 18 23 24]
 [18 19 24 25]
 [ 1  2  8  9]
 [ 2  3  9 10]
 [ 3  4 10 11]
 [ 4  5 11 12]
 [ 5  6 12 13]
 [ 6  7 13 14]
 [ 7  8 14 15]
 [ 8  9 15 16]
 [ 9 10 16 17]
 [10 11 17 18]
 [11 12 18 19]
 [12 13 19 20]
 [13 14 20 21]
 [14 15 21 22]
 [15 16 22 23]
 [16 17 23 24]
 [17 18 24 25]
 [ 1  2  9 10]
 [ 2  3 10 11]
 [ 3  4 11 12]
 [ 4  5 12 13]
 [ 5  6 13 14]
 [ 6  7 14 15]
 [ 7  8 15 16]
 [ 8  9 16 17]
 [ 9 10 17 18]
 [10 11 18 19]
 [11 12 19 20]
 [12 13 20 21]
 [13 14 21 22]
 [14 15 22 23]
 [15 16 23 24]
 [16 17 24 25]
 [ 1  2 10 11]
 [ 2  3 11 12]
 [ 3  4 12 13]
 [ 4  5 13 14]
 [ 5  6 14 15]
 [ 6  7 15 16]
 [ 7  8 16 17]
 [ 8  9 17 18]
 [ 9 10 18 19]
 [10 11 19 20]
 [11 12 20 21]
 [12 13 21 22]
 [13 14 22 23]
 [14 15 23 24]
 [15 16 24 25]
 [ 1  2 11 12]
 [ 2  3 12 13]
 [ 3  4 13 14]
 [ 4  5 14 15]
 [ 5  6 15 16]
 [ 6  7 16 17]
 [ 7  8 17 18]
 [ 8  9 18 19]
 [ 9 10 19 20]
 [10 11 20 21]
 [11 12 21 22]
 [12 13 22 23]
 [13 14 23 24]
 [14 15 24 25]
 [ 1  2 12 13]
 [ 2  3 13 14]
 [ 3  4 14 15]
 [ 4  5 15 16]
 [ 5  6 16 17]
 [ 6  7 17 18]
 [ 7  8 18 19]
 [ 8  9 19 20]
 [ 9 10 20 21]
 [10 11 21 22]
 [11 12 22 23]
 [12 13 23 24]
 [13 14 24 25]]

Then comes the forward modelling part itself. The forward modelling will run R2, cR2, ... in forward mode inside a fwd directory inside the working directory. The resulting apparent resistivity are then embeded inside a Survey object and directly available for inversion for instance.

In [7]:
k.forward(noise=0.05, iplot=True) # forward modelling with 5 % noise added to the output
Writing .in file and mesh.dat... done!
Writing protocol.dat... done!
Running forward model... 

 >> R  2    R e s i s t i v i t y   I n v e r s i o n   v4.0 <<

 >> D a t e : 06 - 08 - 2020
 >> My beautiful survey                                                             
 >> F o r w a r d   S o l u t i o n   S e l e c t e d <<
 >> Determining storage needed for finite element conductance matrix
 >> Generating index array for finite element conductance matrix
 >> Reading start resistivity from resistivity.dat               

 Measurements read:   175     Measurements rejected:     0

 >> Total Memory required is:          0.035 Gb
0/175 reciprocal measurements found.
0/175 reciprocal measurements found.
Forward modelling done.
In [8]:
k.param['b_wgt'] = 0.05 # we are now informing R2 that the relative error is 0.05 for inversion
k.param['a_wgt'] = 0.0 # the offset error is zero, to be consistent with the error added to the forward model

We can already see that the pseudo-section already show clearly the import on the region we defined. We can now invert these apparent resistivities. Inverting the forward models allow the user to see if the parameters of the surveys (the sequence and electrode spacing) were optimium to resolve the target. If needed the user can change them and do the whole process again.

In [9]:
k.invert()
Writing .in file and protocol.dat... All non fixed parameters reset to 100 Ohm.m and 0 mrad, as the survey to be inverted is from a forward model.
done!
--------------------- MAIN INVERSION ------------------


 >> R  2    R e s i s t i v i t y   I n v e r s i o n   v4.0 <<

 >> D a t e : 06 - 08 - 2020
 >> My beautiful survey                                                             
 >> I n v e r s e   S o l u t i o n   S e l e c t e d <<
 >> Determining storage needed for finite element conductance matrix
 >> Generating index array for finite element conductance matrix
 >> Reading start resistivity from res0.dat                      
 >> R e g u l a r i s e d   T y p e <<
 >>   L i n e a r    F i l t e r    <<
 >> L o g - D a t a   I n v e r s i o n <<
 >> N o r m a l   R e g u l a r i s a t i o n <<
 >> D a t a   w e i g h t s   w i l l   b e  m o d i f i e d <<


 Processing dataset   1


 Measurements read:   175     Measurements rejected:     0
   Geometric mean of apparent resistivities:  0.18258E+03

 >> Total Memory required is:          0.039 Gb

   Iteration   1
     Initial RMS Misfit:        13.67       Number of data ignored:     0
     Alpha:         420.865   RMS Misfit:        2.04  Roughness:        1.898
     Alpha:         195.348   RMS Misfit:        1.46  Roughness:        3.221
     Alpha:          90.673   RMS Misfit:        1.15  Roughness:        4.778
     Alpha:          42.087   RMS Misfit:        1.12  Roughness:        6.424
     Alpha:          19.535   RMS Misfit:        1.23  Roughness:        8.074
     Step length set to   1.00000
     Final RMS Misfit:        1.12
     Attempted to update data weights and caused overshoot
     treating as converged

 Solution converged - Outputing results to file

 Calculating sensitivity map


 Processing dataset   2


 End of data:  Terminating
1/1 results parsed (1 ok; 0 failed)
In [10]:
k.showResults(index=0, color_map='jet') # show the initial model
k.showResults(index=1, sens=False, color_map='jet') # show the inverted model
ERROR: No sensitivity attribute found
In [ ]: