Create a 3D layer model

[zhichaoluan]

from pygimli.physics import ert
import pygimli as pg
import numpy as np
import matplotlib.pyplot as plt
import pygimli.meshtools as mt

定义模型层和标记

top_layer = mt.createWorld(start=[-10, -10, 0], end=[140, 10, -10], marker=1)
middle_layer = mt.createWorld(start=[-10, -10, -10], end=[140, 10, -20], marker=2)
bottom_layer = mt.createWorld(start=[-10, -10, -20], end=[140, 10, -40], marker=3)
cube_region = mt.createCube(size=[128, 0.4, 40], pos=[128/2, 0, -20], marker=4)

合并模型区域

start_model = mt.merge([top_layer, middle_layer, bottom_layer, cube_region])

展示初始模型

pg.show(start_model, alpha=0.3, markers=True)

加载数据

filename = "cc.shm"
cc.zip

shm = pg.DataContainerERT(filename)
data = ert.load(filename)
data['k'] = ert.geometricFactors(data)

计算视电阻率

data['rhoa'] = data['k'] * data['u'] / data['i'] * 1000
data.remove(data['rhoa'] < 0)
data['err'] = ert.estimateError(data, relativeError=0.02)

在模型中添加传感器位置

for s in shm.sensors():
start_model.createNode(s)
for s in shm.sensorPositions():
start_model.createNode(s - [0, 0, 1e-2 / 2])

创建反演网格

inversion_mesh = mt.createMesh(start_model, quality=1.4)
pg.show(inversion_mesh, markers=True, showMesh=True)

配置ERT管理器

mgr = ert.ERTManager()
mgr.setData(data)
mgr.setMesh(inversion_mesh)

设置正则化参数

mgr.inv.setRegularization(1, fix=100) # 顶层
mgr.inv.setRegularization(2, fix=200) # 中层
mgr.inv.setRegularization(3, fix=300) # 底层
mgr.inv.setRegularization(4, startModel=1e4, limits=[100, 1e4 + 20]) # 立方体区域

执行反演

mgr.invert(lam=10, zWeight=0.1, verbose=True)

保存结果并展示误差

mgr.saveResult()
mgr.showMisfit()

I want to create a layer model and assign different resistivity values to each for constraint inversion, but I have some problems to help me solve.

Traceback (most recent call last):
File "D:/wangkangbo/test/11111.py", line 36, in
mgr.invert(lam=10,zWeight=0.1,verbose=True)
File "D:\Anaconda3\envs\pg1\lib\site-packages\pygimli\frameworks\methodManager.py", line 777, in invert
if self.fop.mesh() is None:
File "D:\Anaconda3\envs\pg1\lib\site-packages\pygimli\frameworks\modelling.py", line 665, in mesh
self._applyRegionProperties()
File "D:\Anaconda3\envs\pg1\lib\site-packages\pygimli\frameworks\modelling.py", line 325, in _applyRegionProperties
if rMgr.region(rID).fixValue() != vals['fix']:
RuntimeError: C:/msys64/home/halbm/gimli/gimli/core/src/regionManager.cpp:593 GIMLI::Region* GIMLI::RegionManager::region(GIMLI::SIndex) no region with marker 1

[halbmy]

Dear colleague,
when raising an issue, there is a template that should be filled describing some background, what you did, what you expected and what actually happened. You are copying your code (not even into the part where it should go) along with chinese symbols, no images, but an error message that is clear: You do not have a region with marker 1.

[halbmy]

You should use only one call of mt.createWorld(). There you can specify layers by the layers argument, see
https://www.pygimli.org/pygimliapi/_generated/pygimli.meshtools.html#pygimli.meshtools.createWorld
and the below linked examples (mostly 2D) using this function.

[zhichaoluan]

"Sorry, my explanation was not very clear."

"I want to generate a 3D layered model and set an anomaly in it to depict the vertical membrane structures in a landfill. During inversion, I want to keep the resistivity of the background model constant and only invert the central vertical membrane area. However, I'm not sure how to set up the 3D layered model, as it seems the 'layer' parameter can't be used in three dimensions. I'm encountering the following error messages in my code."

[halbmy]

Have you tried mt.createWorld(start=[-10, -10, -40], end=[140, 10, 0], layers=[-10, -30])?

Note that you will need a much larger mesh for accurate simulations.

[zhichaoluan]

"It displays an error like this,"

[halbmy]

Yes, I am sorry, obviously this feature does not exist yet. I would suggest the following:

• creating a world as written but without layers (marker 1)
• adding a cube as a middle layer (marker 2)
• setting a region marker (marker 3) into the first layer

[zhichaoluan]

from pygimli.physics import ert
import pygimli as pg
import numpy as np
import matplotlib.pyplot as plt
import pygimli.meshtools as mt
world=mt.createWorld(start=[-10, -10,-40], end=[140, 10,0],worldMarker=True,marker=1)
cube3=mt.createCube(size=[128, 0.4, 40], pos=[128/2, 0, -40/2],marker=4)
start_model = world +cube3
pg.show(start_model,alpha=0.3,markers=True)

filename = "cc.shm"
shm = pg.DataContainerERT(filename)
data=ert.load(filename)
data['k']=ert.geometricFactors(data)
data['rhoa']=data['k']*data['u']/data['i']*1000
data.remove(data['rhoa']<0)
data['err'] = ert.estimateError(data,relativeError=0.02)
for s in shm.sensors():
start_model.createNode(s)
for s in shm.sensorPositions():
start_model.createNode(s - [0,0,1e-2/2])
inversion_mesh= mt.createMesh(start_model,quality=1.4)
pg.show(inversion_mesh, markers=True, showMesh=True)
mgr=ert.ERTManager()
mgr.setData(data)
mgr.setMesh(inversion_mesh)
mgr.inv.setRegularization(1, background=True)
mgr.inv.setRegularization(1,fix=200)
mgr.inv.setRegularization(2,startModel=1e4,limits=[100,1e4+20])
mgr.invert(lam=10,zWeight=0.1,verbose=True)
mgr.saveResult()
mgr.showMisfit()

I have created a model that includes an anomaly. In the inversion process, I want to constrain the resistivity of the background to remain constant, optimizing only the resistivity of the anomalous cube. Now, I want to change the uniform background model to a layered model, assigning different resistivities to each layer, and similarly constrain the resistivities of the layered background to remain unchanged during inversion

When I haven't established a layered structure, my forward modeling mesh is correct, as shown below

However, when I try to use a layered model, there is a problem with establishing the mesh; the mesh is not divided correctly

[zhichaoluan]

cc.zip

[halbmy]

Please paste code and messages as text instead of screenshots.

In your case the layers (cubes) intersect with the vertical sheet. Whereas such things are automatically resolved in 2D, in 3D this needs to be solved by splitting every cube into two, one on either side of the sheet.