Note
Go to the end to download the full example code.
3D magnetics modelling and inversion#
Based on the synthetic model of Li & Oldenburg (1996), we demonstrate 3D inversion of magnetic data. The forward operator bases on the formula given by Holstein et al. (2007).
In the following, we will build the model, create synthetic data, and do inversion using a depth-weighting function as outlined in the paper.
import numpy as np
import matplotlib.pyplot as plt
import pygimli as pg
from pygimli.viewer import pv
from pygimli.physics.gravimetry import MagneticsModelling
Synthetic model and data generation#
The grid is 1x1x0.5 km with a spacing of 50 m.
Mesh: Nodes: 4851 Cells: 4000 Boundaries: 12800
We create a 3D matrix that is later filled as vector into the grid. The model consists of zeros and patches of 5x6 cells per depth slice that are shifted by one cell for subsequent cells.
We show the model making use of the pyvista package that can be called by
pygimli.show(). The mesh itself is shown as a wireframe, the anomaly
is plotted as a surface plot using a threshold filter. After the first call
with hold=True, the plotter is used to draw any subsequent plots that can
also be slices or clips. Moreover, the camera position is set so that the
vertical axis is going downwards (x is Northing and y is Easting as common in
magnetics).
pl, _ = pg.show(grid, style="wireframe", hold=True)
pv.drawMesh(pl, grid, label="synth", style="surface", cMap="Spectral_r",
filter={"threshold": dict(value=0.05, scalars="synth")})
pl.camera_position = "yz"
pl.camera.roll = 90
pl.camera.azimuth = 180 - 15
pl.camera.elevation = 10
pl.camera.zoom(1.2)
_ = pl.show()
For the computation of the total field, we define the global magnetic
field using the IGRF (total field, inclination and declination) settings
given in the paper. Any global field can also be retrieved by the
pyIGRF module.
F, I, D = 50000, 75, 25 # total field in nT
H = F * np.cos(np.deg2rad(I))
X = H * np.cos(np.deg2rad(D))
Y = H * np.sin(np.deg2rad(D))
Z = F * np.sin(np.deg2rad(I))
igrf = [D, I, H, X, Y, Z, F]
# Alternatively one could use pyIGRF at a specific position
# import pyIGRF
# igrf = pyIGRF.igrf_value(lat=50.59465, lon=12.64139)
py, px = np.meshgrid(x, y)
px = px.ravel()
py = py.ravel()
points = np.column_stack((px, py, np.ones_like(px)*20))
# The forward operator
cmp = ["TFA"] # ["Bx", "By", "Bz"]
fop = MagneticsModelling(mesh=grid, points=points, cmp=cmp, igrf=igrf)
model = pg.Vector(grid.cellCount(), 1.0)
data = fop.response(grid["synth"])
[ 0% ] 1 of 441 complete
[ 0% ] 2 of 441 complete
[ 1% ] 3 of 441 complete
[ 1% ] 4 of 441 complete
[ 1% ] 5 of 441 complete
[ 1% ] 6 of 441 complete
[+ 2% ] 7 of 441 complete
[+ 2% ] 8 of 441 complete
[+ 2% ] 9 of 441 complete
[+ 2% ] 10 of 441 complete
[+ 2% ] 11 of 441 complete
[+ 3% ] 12 of 441 complete
[+ 3% ] 13 of 441 complete
[+ 3% ] 14 of 441 complete
[+ 3% ] 15 of 441 complete
[++ 4% ] 16 of 441 complete
[++ 4% ] 17 of 441 complete
[++ 4% ] 18 of 441 complete
[++ 4% ] 19 of 441 complete
[++ 5% ] 20 of 441 complete
[++ 5% ] 21 of 441 complete
[++ 5% ] 22 of 441 complete
[++ 5% ] 23 of 441 complete
[++ 5% ] 24 of 441 complete
[++ 6% ] 25 of 441 complete
[++ 6% ] 26 of 441 complete
[++ 6% ] 27 of 441 complete
[++ 6% ] 28 of 441 complete
[+++ 7% ] 29 of 441 complete
[+++ 7% ] 30 of 441 complete
[+++ 7% ] 31 of 441 complete
[+++ 7% ] 32 of 441 complete
[+++ 7% ] 33 of 441 complete
[+++ 8% ] 34 of 441 complete
[+++ 8% ] 35 of 441 complete
[+++ 8% ] 36 of 441 complete
[+++ 8% ] 37 of 441 complete
[+++ 9% ] 38 of 441 complete
[+++ 9% ] 39 of 441 complete
[+++ 9% ] 40 of 441 complete
[+++ 9% ] 41 of 441 complete
[++++ 10% ] 42 of 441 complete
[++++ 10% ] 43 of 441 complete
[++++ 10% ] 44 of 441 complete
[++++ 10% ] 45 of 441 complete
[++++ 10% ] 46 of 441 complete
[++++ 11% ] 47 of 441 complete
[++++ 11% ] 48 of 441 complete
[++++ 11% ] 49 of 441 complete
[++++ 11% ] 50 of 441 complete
[+++++ 12% ] 51 of 441 complete
[+++++ 12% ] 52 of 441 complete
[+++++ 12% ] 53 of 441 complete
[+++++ 12% ] 54 of 441 complete
[+++++ 12% ] 55 of 441 complete
[+++++ 13% ] 56 of 441 complete
[+++++ 13% ] 57 of 441 complete
[+++++ 13% ] 58 of 441 complete
[+++++ 13% ] 59 of 441 complete
[+++++ 14% ] 60 of 441 complete
[+++++ 14% ] 61 of 441 complete
[+++++ 14% ] 62 of 441 complete
[+++++ 14% ] 63 of 441 complete
[++++++ 15% ] 64 of 441 complete
[++++++ 15% ] 65 of 441 complete
[++++++ 15% ] 66 of 441 complete
[++++++ 15% ] 67 of 441 complete
[++++++ 15% ] 68 of 441 complete
[++++++ 16% ] 69 of 441 complete
[++++++ 16% ] 70 of 441 complete
[++++++ 16% ] 71 of 441 complete
[++++++ 16% ] 72 of 441 complete
[++++++ 17% ] 73 of 441 complete
[++++++ 17% ] 74 of 441 complete
[++++++ 17% ] 75 of 441 complete
[++++++ 17% ] 76 of 441 complete
[++++++ 17% ] 77 of 441 complete
[+++++++ 18% ] 78 of 441 complete
[+++++++ 18% ] 79 of 441 complete
[+++++++ 18% ] 80 of 441 complete
[+++++++ 18% ] 81 of 441 complete
[+++++++ 19% ] 82 of 441 complete
[+++++++ 19% ] 83 of 441 complete
[+++++++ 19% ] 84 of 441 complete
[+++++++ 19% ] 85 of 441 complete
[++++++++ 20% ] 86 of 441 complete
[++++++++ 20% ] 87 of 441 complete
[++++++++ 20% ] 88 of 441 complete
[++++++++ 20% ] 89 of 441 complete
[++++++++ 20% ] 90 of 441 complete
[++++++++ 21% ] 91 of 441 complete
[++++++++ 21% ] 92 of 441 complete
[++++++++ 21% ] 93 of 441 complete
[++++++++ 21% ] 94 of 441 complete
[++++++++ 22% ] 95 of 441 complete
[++++++++ 22% ] 96 of 441 complete
[++++++++ 22% ] 97 of 441 complete
[++++++++ 22% ] 98 of 441 complete
[++++++++ 22% ] 99 of 441 complete
[+++++++++ 23% ] 100 of 441 complete
[+++++++++ 23% ] 101 of 441 complete
[+++++++++ 23% ] 102 of 441 complete
[+++++++++ 23% ] 103 of 441 complete
[+++++++++ 24% ] 104 of 441 complete
[+++++++++ 24% ] 105 of 441 complete
[+++++++++ 24% ] 106 of 441 complete
[+++++++++ 24% ] 107 of 441 complete
[+++++++++ 24% ] 108 of 441 complete
[++++++++++ 25% ] 109 of 441 complete
[++++++++++ 25% ] 110 of 441 complete
[++++++++++ 25% ] 111 of 441 complete
[++++++++++ 25% ] 112 of 441 complete
[++++++++++ 26% ] 113 of 441 complete
[++++++++++ 26% ] 114 of 441 complete
[++++++++++ 26% ] 115 of 441 complete
[++++++++++ 26% ] 116 of 441 complete
[++++++++++ 27% ] 117 of 441 complete
[++++++++++ 27% ] 118 of 441 complete
[++++++++++ 27% ] 119 of 441 complete
[++++++++++ 27% ] 120 of 441 complete
[++++++++++ 27% ] 121 of 441 complete
[+++++++++++ 28% ] 122 of 441 complete
[+++++++++++ 28% ] 123 of 441 complete
[+++++++++++ 28% ] 124 of 441 complete
[+++++++++++ 28% ] 125 of 441 complete
[+++++++++++ 29% ] 126 of 441 complete
[+++++++++++ 29% ] 127 of 441 complete
[+++++++++++ 29% ] 128 of 441 complete
[+++++++++++ 29% ] 129 of 441 complete
[+++++++++++ 29% ] 130 of 441 complete
[+++++++++++ 30% ] 131 of 441 complete
[+++++++++++ 30% ] 132 of 441 complete
[+++++++++++ 30% ] 133 of 441 complete
[+++++++++++ 30% ] 134 of 441 complete
[++++++++++++ 31% ] 135 of 441 complete
[++++++++++++ 31% ] 136 of 441 complete
[++++++++++++ 31% ] 137 of 441 complete
[++++++++++++ 31% ] 138 of 441 complete
[++++++++++++ 32% ] 139 of 441 complete
[++++++++++++ 32% ] 140 of 441 complete
[++++++++++++ 32% ] 141 of 441 complete
[++++++++++++ 32% ] 142 of 441 complete
[++++++++++++ 32% ] 143 of 441 complete
[+++++++++++++ 33% ] 144 of 441 complete
[+++++++++++++ 33% ] 145 of 441 complete
[+++++++++++++ 33% ] 146 of 441 complete
[+++++++++++++ 33% ] 147 of 441 complete
[+++++++++++++ 34% ] 148 of 441 complete
[+++++++++++++ 34% ] 149 of 441 complete
[+++++++++++++ 34% ] 150 of 441 complete
[+++++++++++++ 34% ] 151 of 441 complete
[+++++++++++++ 34% ] 152 of 441 complete
[+++++++++++++ 35% ] 153 of 441 complete
[+++++++++++++ 35% ] 154 of 441 complete
[+++++++++++++ 35% ] 155 of 441 complete
[+++++++++++++ 35% ] 156 of 441 complete
[++++++++++++++ 36% ] 157 of 441 complete
[++++++++++++++ 36% ] 158 of 441 complete
[++++++++++++++ 36% ] 159 of 441 complete
[++++++++++++++ 36% ] 160 of 441 complete
[++++++++++++++ 37% ] 161 of 441 complete
[++++++++++++++ 37% ] 162 of 441 complete
[++++++++++++++ 37% ] 163 of 441 complete
[++++++++++++++ 37% ] 164 of 441 complete
[++++++++++++++ 37% ] 165 of 441 complete
[++++++++++++++ 38% ] 166 of 441 complete
[++++++++++++++ 38% ] 167 of 441 complete
[++++++++++++++ 38% ] 168 of 441 complete
[++++++++++++++ 38% ] 169 of 441 complete
[+++++++++++++++ 39% ] 170 of 441 complete
[+++++++++++++++ 39% ] 171 of 441 complete
[+++++++++++++++ 39% ] 172 of 441 complete
[+++++++++++++++ 39% ] 173 of 441 complete
[+++++++++++++++ 39% ] 174 of 441 complete
[+++++++++++++++ 40% ] 175 of 441 complete
[+++++++++++++++ 40% ] 176 of 441 complete
[+++++++++++++++ 40% ] 177 of 441 complete
[+++++++++++++++ 40% ] 178 of 441 complete
[++++++++++++++++ 41% ] 179 of 441 complete
[++++++++++++++++ 41% ] 180 of 441 complete
[++++++++++++++++ 41% ] 181 of 441 complete
[++++++++++++++++ 41% ] 182 of 441 complete
[++++++++++++++++ 41% ] 183 of 441 complete
[++++++++++++++++ 42% ] 184 of 441 complete
[++++++++++++++++ 42% ] 185 of 441 complete
[++++++++++++++++ 42% ] 186 of 441 complete
[++++++++++++++++ 42% ] 187 of 441 complete
[++++++++++++++++ 43% ] 188 of 441 complete
[++++++++++++++++ 43% ] 189 of 441 complete
[++++++++++++++++ 43% ] 190 of 441 complete
[++++++++++++++++ 43% ] 191 of 441 complete
[+++++++++++++++++ 44% ] 192 of 441 complete
[+++++++++++++++++ 44% ] 193 of 441 complete
[+++++++++++++++++ 44% ] 194 of 441 complete
[+++++++++++++++++ 44% ] 195 of 441 complete
[+++++++++++++++++ 44% ] 196 of 441 complete
[+++++++++++++++++ 45% ] 197 of 441 complete
[+++++++++++++++++ 45% ] 198 of 441 complete
[+++++++++++++++++ 45% ] 199 of 441 complete
[+++++++++++++++++ 45% ] 200 of 441 complete
[+++++++++++++++++ 46% ] 201 of 441 complete
[+++++++++++++++++ 46% ] 202 of 441 complete
[+++++++++++++++++ 46% ] 203 of 441 complete
[+++++++++++++++++ 46% ] 204 of 441 complete
[+++++++++++++++++ 46% ] 205 of 441 complete
[+++++++++++++++++ 47% ] 206 of 441 complete
[+++++++++++++++++ 47% ] 207 of 441 complete
[+++++++++++++++++ 47% ] 208 of 441 complete
[+++++++++++++++++ 47% ] 209 of 441 complete
[+++++++++++++++++ 48% ] 210 of 441 complete
[+++++++++++++++++ 48% ] 211 of 441 complete
[+++++++++++++++++ 48% ] 212 of 441 complete
[+++++++++++++++++ 48% ] 213 of 441 complete
[+++++++++++++++++ 49% ] 214 of 441 complete
[+++++++++++++++++ 49% ] 215 of 441 complete
[+++++++++++++++++ 49% ] 216 of 441 complete
[+++++++++++++++++ 49% ] 217 of 441 complete
[+++++++++++++++++ 49% ] 218 of 441 complete
[+++++++++++++++++ 50% ] 219 of 441 complete
[+++++++++++++++++ 50% ] 220 of 441 complete
[+++++++++++++++++ 50% ] 221 of 441 complete
[+++++++++++++++++ 50% ] 222 of 441 complete
[+++++++++++++++++ 51% ] 223 of 441 complete
[+++++++++++++++++ 51% ] 224 of 441 complete
[+++++++++++++++++ 51% ] 225 of 441 complete
[+++++++++++++++++ 51% ] 226 of 441 complete
[+++++++++++++++++ 51% ] 227 of 441 complete
[+++++++++++++++++ 52% ] 228 of 441 complete
[+++++++++++++++++ 52% ] 229 of 441 complete
[+++++++++++++++++ 52% ] 230 of 441 complete
[+++++++++++++++++ 52% ] 231 of 441 complete
[+++++++++++++++++ 53% ] 232 of 441 complete
[+++++++++++++++++ 53% ] 233 of 441 complete
[+++++++++++++++++ 53% ] 234 of 441 complete
[+++++++++++++++++ 53% ] 235 of 441 complete
[+++++++++++++++++ 54% ] 236 of 441 complete
[+++++++++++++++++ 54% ] 237 of 441 complete
[+++++++++++++++++ 54% ] 238 of 441 complete
[+++++++++++++++++ 54% ] 239 of 441 complete
[+++++++++++++++++ 54% ] 240 of 441 complete
[+++++++++++++++++ 55% ] 241 of 441 complete
[+++++++++++++++++ 55% ] 242 of 441 complete
[+++++++++++++++++ 55% ] 243 of 441 complete
[+++++++++++++++++ 55% ] 244 of 441 complete
[+++++++++++++++++ 56% ] 245 of 441 complete
[+++++++++++++++++ 56% ] 246 of 441 complete
[+++++++++++++++++ 56% ] 247 of 441 complete
[+++++++++++++++++ 56% ] 248 of 441 complete
[+++++++++++++++++ 56% ] 249 of 441 complete
[+++++++++++++++++ 57% ] 250 of 441 complete
[+++++++++++++++++ 57% ] 251 of 441 complete
[+++++++++++++++++ 57% ] 252 of 441 complete
[+++++++++++++++++ 57% ] 253 of 441 complete
[+++++++++++++++++ 58% ] 254 of 441 complete
[+++++++++++++++++ 58% ] 255 of 441 complete
[+++++++++++++++++ 58% ] 256 of 441 complete
[+++++++++++++++++ 58% ] 257 of 441 complete
[+++++++++++++++++ 59% ] 258 of 441 complete
[+++++++++++++++++ 59% ] 259 of 441 complete
[+++++++++++++++++ 59% ] 260 of 441 complete
[+++++++++++++++++ 59% ] 261 of 441 complete
[+++++++++++++++++ 59% ] 262 of 441 complete
[+++++++++++++++++ 60% + ] 263 of 441 complete
[+++++++++++++++++ 60% + ] 264 of 441 complete
[+++++++++++++++++ 60% + ] 265 of 441 complete
[+++++++++++++++++ 60% + ] 266 of 441 complete
[+++++++++++++++++ 61% + ] 267 of 441 complete
[+++++++++++++++++ 61% + ] 268 of 441 complete
[+++++++++++++++++ 61% + ] 269 of 441 complete
[+++++++++++++++++ 61% + ] 270 of 441 complete
[+++++++++++++++++ 61% + ] 271 of 441 complete
[+++++++++++++++++ 62% ++ ] 272 of 441 complete
[+++++++++++++++++ 62% ++ ] 273 of 441 complete
[+++++++++++++++++ 62% ++ ] 274 of 441 complete
[+++++++++++++++++ 62% ++ ] 275 of 441 complete
[+++++++++++++++++ 63% ++ ] 276 of 441 complete
[+++++++++++++++++ 63% ++ ] 277 of 441 complete
[+++++++++++++++++ 63% ++ ] 278 of 441 complete
[+++++++++++++++++ 63% ++ ] 279 of 441 complete
[+++++++++++++++++ 63% ++ ] 280 of 441 complete
[+++++++++++++++++ 64% ++ ] 281 of 441 complete
[+++++++++++++++++ 64% ++ ] 282 of 441 complete
[+++++++++++++++++ 64% ++ ] 283 of 441 complete
[+++++++++++++++++ 64% ++ ] 284 of 441 complete
[+++++++++++++++++ 65% +++ ] 285 of 441 complete
[+++++++++++++++++ 65% +++ ] 286 of 441 complete
[+++++++++++++++++ 65% +++ ] 287 of 441 complete
[+++++++++++++++++ 65% +++ ] 288 of 441 complete
[+++++++++++++++++ 66% +++ ] 289 of 441 complete
[+++++++++++++++++ 66% +++ ] 290 of 441 complete
[+++++++++++++++++ 66% +++ ] 291 of 441 complete
[+++++++++++++++++ 66% +++ ] 292 of 441 complete
[+++++++++++++++++ 66% +++ ] 293 of 441 complete
[+++++++++++++++++ 67% +++ ] 294 of 441 complete
[+++++++++++++++++ 67% +++ ] 295 of 441 complete
[+++++++++++++++++ 67% +++ ] 296 of 441 complete
[+++++++++++++++++ 67% +++ ] 297 of 441 complete
[+++++++++++++++++ 68% ++++ ] 298 of 441 complete
[+++++++++++++++++ 68% ++++ ] 299 of 441 complete
[+++++++++++++++++ 68% ++++ ] 300 of 441 complete
[+++++++++++++++++ 68% ++++ ] 301 of 441 complete
[+++++++++++++++++ 68% ++++ ] 302 of 441 complete
[+++++++++++++++++ 69% ++++ ] 303 of 441 complete
[+++++++++++++++++ 69% ++++ ] 304 of 441 complete
[+++++++++++++++++ 69% ++++ ] 305 of 441 complete
[+++++++++++++++++ 69% ++++ ] 306 of 441 complete
[+++++++++++++++++ 70% +++++ ] 307 of 441 complete
[+++++++++++++++++ 70% +++++ ] 308 of 441 complete
[+++++++++++++++++ 70% +++++ ] 309 of 441 complete
[+++++++++++++++++ 70% +++++ ] 310 of 441 complete
[+++++++++++++++++ 71% +++++ ] 311 of 441 complete
[+++++++++++++++++ 71% +++++ ] 312 of 441 complete
[+++++++++++++++++ 71% +++++ ] 313 of 441 complete
[+++++++++++++++++ 71% +++++ ] 314 of 441 complete
[+++++++++++++++++ 71% +++++ ] 315 of 441 complete
[+++++++++++++++++ 72% +++++ ] 316 of 441 complete
[+++++++++++++++++ 72% +++++ ] 317 of 441 complete
[+++++++++++++++++ 72% +++++ ] 318 of 441 complete
[+++++++++++++++++ 72% +++++ ] 319 of 441 complete
[+++++++++++++++++ 73% ++++++ ] 320 of 441 complete
[+++++++++++++++++ 73% ++++++ ] 321 of 441 complete
[+++++++++++++++++ 73% ++++++ ] 322 of 441 complete
[+++++++++++++++++ 73% ++++++ ] 323 of 441 complete
[+++++++++++++++++ 73% ++++++ ] 324 of 441 complete
[+++++++++++++++++ 74% ++++++ ] 325 of 441 complete
[+++++++++++++++++ 74% ++++++ ] 326 of 441 complete
[+++++++++++++++++ 74% ++++++ ] 327 of 441 complete
[+++++++++++++++++ 74% ++++++ ] 328 of 441 complete
[+++++++++++++++++ 75% ++++++ ] 329 of 441 complete
[+++++++++++++++++ 75% ++++++ ] 330 of 441 complete
[+++++++++++++++++ 75% ++++++ ] 331 of 441 complete
[+++++++++++++++++ 75% ++++++ ] 332 of 441 complete
[+++++++++++++++++ 76% +++++++ ] 333 of 441 complete
[+++++++++++++++++ 76% +++++++ ] 334 of 441 complete
[+++++++++++++++++ 76% +++++++ ] 335 of 441 complete
[+++++++++++++++++ 76% +++++++ ] 336 of 441 complete
[+++++++++++++++++ 76% +++++++ ] 337 of 441 complete
[+++++++++++++++++ 77% +++++++ ] 338 of 441 complete
[+++++++++++++++++ 77% +++++++ ] 339 of 441 complete
[+++++++++++++++++ 77% +++++++ ] 340 of 441 complete
[+++++++++++++++++ 77% +++++++ ] 341 of 441 complete
[+++++++++++++++++ 78% ++++++++ ] 342 of 441 complete
[+++++++++++++++++ 78% ++++++++ ] 343 of 441 complete
[+++++++++++++++++ 78% ++++++++ ] 344 of 441 complete
[+++++++++++++++++ 78% ++++++++ ] 345 of 441 complete
[+++++++++++++++++ 78% ++++++++ ] 346 of 441 complete
[+++++++++++++++++ 79% ++++++++ ] 347 of 441 complete
[+++++++++++++++++ 79% ++++++++ ] 348 of 441 complete
[+++++++++++++++++ 79% ++++++++ ] 349 of 441 complete
[+++++++++++++++++ 79% ++++++++ ] 350 of 441 complete
[+++++++++++++++++ 80% ++++++++ ] 351 of 441 complete
[+++++++++++++++++ 80% ++++++++ ] 352 of 441 complete
[+++++++++++++++++ 80% ++++++++ ] 353 of 441 complete
[+++++++++++++++++ 80% ++++++++ ] 354 of 441 complete
[+++++++++++++++++ 80% ++++++++ ] 355 of 441 complete
[+++++++++++++++++ 81% +++++++++ ] 356 of 441 complete
[+++++++++++++++++ 81% +++++++++ ] 357 of 441 complete
[+++++++++++++++++ 81% +++++++++ ] 358 of 441 complete
[+++++++++++++++++ 81% +++++++++ ] 359 of 441 complete
[+++++++++++++++++ 82% +++++++++ ] 360 of 441 complete
[+++++++++++++++++ 82% +++++++++ ] 361 of 441 complete
[+++++++++++++++++ 82% +++++++++ ] 362 of 441 complete
[+++++++++++++++++ 82% +++++++++ ] 363 of 441 complete
[+++++++++++++++++ 83% ++++++++++ ] 364 of 441 complete
[+++++++++++++++++ 83% ++++++++++ ] 365 of 441 complete
[+++++++++++++++++ 83% ++++++++++ ] 366 of 441 complete
[+++++++++++++++++ 83% ++++++++++ ] 367 of 441 complete
[+++++++++++++++++ 83% ++++++++++ ] 368 of 441 complete
[+++++++++++++++++ 84% ++++++++++ ] 369 of 441 complete
[+++++++++++++++++ 84% ++++++++++ ] 370 of 441 complete
[+++++++++++++++++ 84% ++++++++++ ] 371 of 441 complete
[+++++++++++++++++ 84% ++++++++++ ] 372 of 441 complete
[+++++++++++++++++ 85% ++++++++++ ] 373 of 441 complete
[+++++++++++++++++ 85% ++++++++++ ] 374 of 441 complete
[+++++++++++++++++ 85% ++++++++++ ] 375 of 441 complete
[+++++++++++++++++ 85% ++++++++++ ] 376 of 441 complete
[+++++++++++++++++ 85% ++++++++++ ] 377 of 441 complete
[+++++++++++++++++ 86% +++++++++++ ] 378 of 441 complete
[+++++++++++++++++ 86% +++++++++++ ] 379 of 441 complete
[+++++++++++++++++ 86% +++++++++++ ] 380 of 441 complete
[+++++++++++++++++ 86% +++++++++++ ] 381 of 441 complete
[+++++++++++++++++ 87% +++++++++++ ] 382 of 441 complete
[+++++++++++++++++ 87% +++++++++++ ] 383 of 441 complete
[+++++++++++++++++ 87% +++++++++++ ] 384 of 441 complete
[+++++++++++++++++ 87% +++++++++++ ] 385 of 441 complete
[+++++++++++++++++ 88% +++++++++++ ] 386 of 441 complete
[+++++++++++++++++ 88% +++++++++++ ] 387 of 441 complete
[+++++++++++++++++ 88% +++++++++++ ] 388 of 441 complete
[+++++++++++++++++ 88% +++++++++++ ] 389 of 441 complete
[+++++++++++++++++ 88% +++++++++++ ] 390 of 441 complete
[+++++++++++++++++ 89% ++++++++++++ ] 391 of 441 complete
[+++++++++++++++++ 89% ++++++++++++ ] 392 of 441 complete
[+++++++++++++++++ 89% ++++++++++++ ] 393 of 441 complete
[+++++++++++++++++ 89% ++++++++++++ ] 394 of 441 complete
[+++++++++++++++++ 90% ++++++++++++ ] 395 of 441 complete
[+++++++++++++++++ 90% ++++++++++++ ] 396 of 441 complete
[+++++++++++++++++ 90% ++++++++++++ ] 397 of 441 complete
[+++++++++++++++++ 90% ++++++++++++ ] 398 of 441 complete
[+++++++++++++++++ 90% ++++++++++++ ] 399 of 441 complete
[+++++++++++++++++ 91% +++++++++++++ ] 400 of 441 complete
[+++++++++++++++++ 91% +++++++++++++ ] 401 of 441 complete
[+++++++++++++++++ 91% +++++++++++++ ] 402 of 441 complete
[+++++++++++++++++ 91% +++++++++++++ ] 403 of 441 complete
[+++++++++++++++++ 92% +++++++++++++ ] 404 of 441 complete
[+++++++++++++++++ 92% +++++++++++++ ] 405 of 441 complete
[+++++++++++++++++ 92% +++++++++++++ ] 406 of 441 complete
[+++++++++++++++++ 92% +++++++++++++ ] 407 of 441 complete
[+++++++++++++++++ 93% +++++++++++++ ] 408 of 441 complete
[+++++++++++++++++ 93% +++++++++++++ ] 409 of 441 complete
[+++++++++++++++++ 93% +++++++++++++ ] 410 of 441 complete
[+++++++++++++++++ 93% +++++++++++++ ] 411 of 441 complete
[+++++++++++++++++ 93% +++++++++++++ ] 412 of 441 complete
[+++++++++++++++++ 94% ++++++++++++++ ] 413 of 441 complete
[+++++++++++++++++ 94% ++++++++++++++ ] 414 of 441 complete
[+++++++++++++++++ 94% ++++++++++++++ ] 415 of 441 complete
[+++++++++++++++++ 94% ++++++++++++++ ] 416 of 441 complete
[+++++++++++++++++ 95% ++++++++++++++ ] 417 of 441 complete
[+++++++++++++++++ 95% ++++++++++++++ ] 418 of 441 complete
[+++++++++++++++++ 95% ++++++++++++++ ] 419 of 441 complete
[+++++++++++++++++ 95% ++++++++++++++ ] 420 of 441 complete
[+++++++++++++++++ 95% ++++++++++++++ ] 421 of 441 complete
[+++++++++++++++++ 96% ++++++++++++++ ] 422 of 441 complete
[+++++++++++++++++ 96% ++++++++++++++ ] 423 of 441 complete
[+++++++++++++++++ 96% ++++++++++++++ ] 424 of 441 complete
[+++++++++++++++++ 96% ++++++++++++++ ] 425 of 441 complete
[+++++++++++++++++ 97% +++++++++++++++ ] 426 of 441 complete
[+++++++++++++++++ 97% +++++++++++++++ ] 427 of 441 complete
[+++++++++++++++++ 97% +++++++++++++++ ] 428 of 441 complete
[+++++++++++++++++ 97% +++++++++++++++ ] 429 of 441 complete
[+++++++++++++++++ 98% +++++++++++++++ ] 430 of 441 complete
[+++++++++++++++++ 98% +++++++++++++++ ] 431 of 441 complete
[+++++++++++++++++ 98% +++++++++++++++ ] 432 of 441 complete
[+++++++++++++++++ 98% +++++++++++++++ ] 433 of 441 complete
[+++++++++++++++++ 98% +++++++++++++++ ] 434 of 441 complete
[+++++++++++++++++ 99% ++++++++++++++++] 435 of 441 complete
[+++++++++++++++++ 99% ++++++++++++++++] 436 of 441 complete
[+++++++++++++++++ 99% ++++++++++++++++] 437 of 441 complete
[+++++++++++++++++ 99% ++++++++++++++++] 438 of 441 complete
[++++++++++++++++ 100% ++++++++++++++++] 439 of 441 complete
[++++++++++++++++ 100% ++++++++++++++++] 440 of 441 complete
[++++++++++++++++ 100% ++++++++++++++++] 441 of 441 complete
Just like in the paper, the data are contaminated with an error model consisting of relative and absolute noise of 2% and 1 nT, respectively.
absError = np.abs(data) * 0.02 + 1
data += np.random.randn(len(data)) * absError
Depth weighting#
In the paper of Li & Oldenburg (1996), they propose a depth weighting of the constraints with the formula
\[w_z = \frac{1}{(z+z_0)^{3/2}}\]
0.008944271909999158 0.282842712474619
Inversion#
The inversion is rather straightforward using the standard inversion
framework pygimli.Inversion.
fop: <pygimli.physics.gravimetry.MagneticsModelling.MagneticsModelling object at 0x145211c60>
Data transformation: Identity transform
Model transformation (cumulative):
0 Logarithmic LU transform, lower bound 0.0, upper bound 0.07
min/max (data): -18.31/354
min/max (error): 2.24%/618%
min/max (start model): 1.0e-03/1.0e-03
--------------------------------------------------------------------------------
inv.iter 0 ... chi² = 378.11
--------------------------------------------------------------------------------
inv.iter 1 ... chi² = 361.17 (dPhi = 4.33%) lam: 200.0
--------------------------------------------------------------------------------
inv.iter 2 ... chi² = 96.88 (dPhi = 72.71%) lam: 200.0
--------------------------------------------------------------------------------
inv.iter 3 ... chi² = 22.20 (dPhi = 75.18%) lam: 200.0
--------------------------------------------------------------------------------
inv.iter 4 ... chi² = 14.74 (dPhi = 30.63%) lam: 200.0
--------------------------------------------------------------------------------
inv.iter 5 ... chi² = 8.59 (dPhi = 34.94%) lam: 200.0
--------------------------------------------------------------------------------
inv.iter 6 ... chi² = 1.27 (dPhi = 65.88%) lam: 200.0
--------------------------------------------------------------------------------
inv.iter 7 ... chi² = 1.10 (dPhi = 3.96%) lam: 200.0
--------------------------------------------------------------------------------
inv.iter 8 ... chi² = 1.10 (dPhi = 0.64%) lam: 200.0
################################################################################
# Abort criterion reached: dPhi = 0.64 (< 2.0%) #
################################################################################
Visualization#
For showing the model, we again use the threshold filter with a value of 0.02. For comparison with the synthetic model, we plot the latter as a wireframe.
pl, _ = pg.show(grid, label="synth", style="wireframe", hold=True,
filter={"threshold": dict(value=0.025, scalars="synth")})
pv.drawMesh(pl, grid, label="inv", style="surface", cMap="Spectral_r",
filter={"threshold": dict(value=0.02, scalars="inv")})
pv.drawMesh(pl, grid, label="inv", style="surface", cMap="Spectral_r",
filter={"slice": dict(normal=[-1, 0, 0], origin=[500, 600, 250])})
pl.camera_position = "yz"
pl.camera.roll = 90
pl.camera.azimuth = 180 - 15
pl.camera.elevation = 10
pl.camera.zoom(1.2)
_ = pl.show()
The model can nicely outline the top part of the anomalous body, but not its depth extent.
We compare the data and model response by means of scatter plots:
fig, ax = plt.subplots(ncols=2, figsize=(12, 5), sharex=True, sharey=True)
vals = data
mm = np.max(np.abs(vals))
ax[0].scatter(px, py, c=vals, cmap="bwr", vmin=-mm, vmax=mm)
im = ax[1].scatter(px, py, c=inv.response, cmap="bwr", vmin=-mm, vmax=mm)
cb = plt.colorbar(im, ax=ax[1])
cb.set_label("B (nT)")
Alternatively, we can also plot the error-weighted misfit.
misfit = (inv.response - data) / absError
im = plt.scatter(py, px, c=misfit, cmap="bwr", vmin=-3, vmax=3)
cb = plt.colorbar(im)
cb.set_label("misfit / error")
References#
Li, Y. & Oldenburg, D. (1996): 3-D inversion of magnetic data. Geophysics 61(2), 394-408.
Holstein, H., Sherratt, E.M., Reid, A.B. (2007): Gravimagnetic field tensor gradiometry formulas for uniform polyhedra, SEG Ext. Abstr.
Total running time of the script: (0 minutes 10.402 seconds)