# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
from SimPEG import maps
from SimPEG.potential_fields import magnetics
from discretize import TensorMesh
import os
print("=== Magnetic Forward Modeling for Coal Fire Detection ===")
output_dir = "Magnetic_Results"
os.makedirs(output_dir, exist_ok=True)
print(f"Results will be saved to: {output_dir}")
def create_model1_horizontal_slab():
model = np.zeros(mesh.nC)
centers = mesh.cell_centers
indices = (
(centers[:, 0] >= -80) & (centers[:, 0] <= 80) &
(centers[:, 2] >= -60) & (centers[:, 2] <= -40) &
(centers[:, 1] >= -5) & (centers[:, 1] <= 5)
)
model[indices] = 0.08
return model, "01_Horizontal_Slab"
def create_model2_dipping_slab():
model = np.zeros(mesh.nC)
centers = mesh.cell_centers
angle = np.radians(30)
rotated_x = centers[:, 0] * np.cos(angle) - centers[:, 2] * np.sin(angle)
rotated_z = centers[:, 0] * np.sin(angle) + centers[:, 2] * np.cos(angle)
indices = (
(rotated_x >= -60) & (rotated_x <= 60) &
(rotated_z >= -10) & (rotated_z <= 10) &
(centers[:, 1] >= -5) & (centers[:, 1] <= 5)
)
model[indices] = 0.07
return model, "02_Dipping_Slab"
def create_model3_rectangular_prism():
model = np.zeros(mesh.nC)
centers = mesh.cell_centers
indices = (
(centers[:, 0] >= -40) & (centers[:, 0] <= 40) &
(centers[:, 2] >= -100) & (centers[:, 2] <= -60) &
(centers[:, 1] >= -5) & (centers[:, 1] <= 5)
)
model[indices] = 0.06
return model, "03_Rectangular_Prism"
def create_model4_two_anomalies():
model = np.zeros(mesh.nC)
centers = mesh.cell_centers
indices1 = (
(centers[:, 0] >= -100) & (centers[:, 0] <= -50) &
(centers[:, 2] >= -50) & (centers[:, 2] <= -30) &
(centers[:, 1] >= -5) & (centers[:, 1] <= 5)
)
indices2 = (
(centers[:, 0] >= 30) & (centers[:, 0] <= 90) &
(centers[:, 2] >= -90) & (centers[:, 2] <= -50) &
(centers[:, 1] >= -5) & (centers[:, 1] <= 5)
)
model[indices1] = 0.09
model[indices2] = 0.05
return model, "04_Two_Anomalies"
def create_model5_complex():
model = np.zeros(mesh.nC)
centers = mesh.cell_centers
indices1 = (
(centers[:, 0] >= -120) & (centers[:, 0] <= -60) &
(centers[:, 2] >= -40) & (centers[:, 2] <= -20) &
(centers[:, 1] >= -5) & (centers[:, 1] <= 5)
)
angle = np.radians(45)
rotated_x = (centers[:, 0] - 50) * np.cos(angle) - (centers[:, 2] + 70) * np.sin(angle)
rotated_z = (centers[:, 0] - 50) * np.sin(angle) + (centers[:, 2] + 70) * np.cos(angle)
indices2 = (
(rotated_x >= -30) & (rotated_x <= 30) &
(rotated_z >= -8) & (rotated_z <= 8) &
(centers[:, 1] >= -5) & (centers[:, 1] <= 5)
)
model[indices1] = 0.08
model[indices2] = 0.065
return model, "05_Complex_Model"
print("Initializing...")
mesh = TensorMesh([
np.ones(100) * 10,
np.ones(1) * 10,
np.ones(80) * 5
], x0=['C', 'C', 0])
print(f"Mesh created: {mesh.nC} cells")
x_locations = np.linspace(-150, 150, 61)
y_locations = np.zeros_like(x_locations)
z_locations = np.ones_like(x_locations) * 10
receiver_locations = np.c_[x_locations, y_locations, z_locations]
H0 = 50000
incl = 50
decl = -5
receivers = magnetics.receivers.Point(receiver_locations, components="tmi")
source_field = magnetics.sources.UniformBackgroundField(
receiver_list=[receivers],
amplitude=H0,
inclination=incl,
declination=decl
)
survey = magnetics.survey.Survey(source_field)
print("Survey system setup completed: 61 stations")
all_models = [
create_model1_horizontal_slab,
create_model2_dipping_slab,
create_model3_rectangular_prism,
create_model4_two_anomalies,
create_model5_complex
]
print("\nStarting forward modeling...")
results = {}
for i, model_func in enumerate(all_models):
print(f"Progress: {i+1}/{len(all_models)} - Calculating model {i+1}...")
model, model_name = model_func()
simulation = magnetics.simulation.Simulation3DIntegral(
survey=survey,
mesh=mesh,
modelType='scalar',
chiMap=maps.IdentityMap(nP=mesh.nC)
)
data = simulation.dpred(model)
results[model_name] = {
'model': model,
'data': data,
'receiver_locations': receiver_locations
}
fig, axes = plt.subplots(2, 2, figsize=(15, 12))
fig.suptitle(f'Magnetic Modeling - {model_name}', fontsize=16)
model_2d = model.reshape(mesh.shape[0], mesh.shape[1], mesh.shape[2])[:, 0, :]
x_centers = mesh.cell_centers_x
z_centers = mesh.cell_centers_z
X, Z = np.meshgrid(x_centers, z_centers, indexing='ij')
im1 = axes[0,0].pcolormesh(X, Z, model_2d.T, cmap='jet', shading='auto')
axes[0,0].set_xlabel('Distance (m)')
axes[0,0].set_ylabel('Depth (m)')
axes[0,0].set_title('Underground Model')
axes[0,0].set_aspect('equal')
plt.colorbar(im1, ax=axes[0,0], label='Susceptibility')
axes[0,0].axhline(y=0, color='k', linestyle='-', linewidth=1)
axes[0,1].plot(receiver_locations[:,0], data, 'b-', linewidth=2, marker='o', markersize=3)
axes[0,1].set_xlabel('Station Position (m)')
axes[0,1].set_ylabel('Magnetic Anomaly (nT)')
axes[0,1].set_title('Surface Magnetic Data')
axes[0,1].grid(True, alpha=0.3)
axes[0,1].fill_between(receiver_locations[:,0], data, alpha=0.3, color='blue')
ax3d = axes[1,0]
centers = mesh.cell_centers
mask = model > 0.01
if np.any(mask):
scatter = ax3d.scatter(centers[mask, 0], centers[mask, 1], centers[mask, 2],
c=model[mask], cmap='jet', s=10)
ax3d.set_xlabel('X (m)')
ax3d.set_ylabel('Y (m)')
ax3d.set_zlabel('Depth (m)')
ax3d.set_title('3D Model View')
axes[1,1].axis('off')
info_text = f"""Model Info:
Name: {model_name}
Data Range: {np.min(data):.1f} ~ {np.max(data):.1f} nT
Peak Anomaly: {np.max(np.abs(data)):.1f} nT
Stations: {len(data)}
Anomaly Cells: {np.sum(model > 0.01)}
Max Susceptibility: {np.max(model):.3f} SI"""
axes[1,1].text(0.1, 0.9, info_text, fontsize=12, verticalalignment='top', linespacing=1.5)
plt.tight_layout()
plt.savefig(f'{output_dir}/{model_name}.png', dpi=300, bbox_inches='tight')
plt.close(fig)
print(f" Saved: {model_name}.png")
print("\nSaving data files...")
np.save(f'{output_dir}/all_models_data.npy', results)
print("Generating report...")
with open(f'{output_dir}/data_report.txt', 'w') as f:
f.write("Magnetic Forward Modeling Dataset Summary\n")
f.write("="*50 + "\n\n")
f.write(f"Generated: {np.datetime64('now')}\n")
f.write(f"Total Models: {len(all_models)}\n")
f.write(f"Stations: {len(receiver_locations)}\n")
f.write(f"Mesh Size: {mesh.nC} cells\n\n")
f.write("Magnetic Anomaly Statistics:\n")
f.write("-"*30 + "\n")
for name in results:
data = results[name]['data']
f.write(f"{name}: {np.min(data):7.1f} ~ {np.max(data):7.1f} nT (Peak: {np.max(np.abs(data)):7.1f} nT)\n")
print("="*50)
print("🎉 Project Completed!")
print(f"✅ Generated {len(all_models)} typical fire zone models")
print(f"✅ All results saved to '{output_dir}' folder")
print(f"✅ Data file: all_models_data.npy")
print(f"✅ Report: data_report.txt")
print("\nYou now have a complete dataset for deep learning!")
print("- Underground models are your [LABELS]")
print("- Surface magnetic data are your [FEATURES]")
最新发布
对象缓存机制
扫一扫
906
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
