✅博主简介：擅长数据搜集与处理、建模仿真、程序设计、仿真代码、论文写作与指导，毕业论文、期刊论文经验交流。
✅成品或者定制，扫描文章底部微信二维码。
（1）二维大地电磁正演与正则化反演方法实现
大地电磁法是一种重要的地球物理勘探技术，通过测量天然电磁场在地下介质中激发的感应电磁响应来推断地下电性结构分布。在进行深度学习反演研究之前，需要首先建立完善的正演计算和传统反演方法作为基础支撑。本研究采用有限单元法实现二维大地电磁正演计算，该方法通过将研究区域离散化为大量三角形或四边形单元，将连续的偏微分方程问题转化为离散的代数方程组进行求解。
在网格剖分方面，考虑到大地电磁场在地表附近和电性界面处变化剧烈的特点，采用非均匀网格剖分策略。在测点位置和预设异常体边界附近进行网格加密，以保证计算精度；而在远离关注区域的边界附近则适当放稀网格密度，以控制计算规模和计算时间。有限单元法的核心是构建单元刚度矩阵和整体刚度矩阵，通过对麦克斯韦方程组进行弱形式变换和分部积分处理，可以得到适用于大地电磁问题的有限元方程。在数值求解过程中，采用直接法求解稀疏线性方程组，得到各节点的电磁场值，进而计算视电阻率和相位等反演所需的观测数据。
正则化反演是大地电磁数据处理的重要环节，其目的是根据地表观测数据重建地下电阻率分布。由于大地电磁反演问题具有严重的非适定性和多解性，需要引入正则化约束来稳定反演过程并获得合理的地质解释。本研究实现的正则化反演算法采用光滑约束作为模型正则化项，通过最小化目标函数来寻找最优模型参数。目标函数由数据拟合项和模型约束项两部分组成，数据拟合项衡量理论响应与观测数据之间的差异，模型约束项则限制模型参数在相邻单元之间的变化幅度，促使反演结果呈现光滑连续的特征。
在反演算法实现过程中，采用高斯牛顿法进行迭代求解。每次迭代需要计算雅可比矩阵，即观测数据对模型参数的敏感度矩阵。雅可比矩阵的计算采用伴随方程法，相比传统的扰动法具有更高的计算效率。模型测试实验表明，正则化反演算法能够有效恢复地下异常体的大致位置和形态，但在异常体边界的刻画和电阻率数值的准确性方面存在一定局限性。此外，当观测数据受到静态效应干扰时，反演结果会出现明显失真，表现为浅部出现虚假的高阻或低阻异常。这些问题为后续采用深度学习方法改进反演效果提供了明确的研究动机。
（2）大地电磁深度学习卷积神经网络设计与训练策略
针对传统正则化反演方法存在的精度和抗干扰能力不足等问题，本研究设计了一种专门用于大地电磁数据反演的深度卷积神经网络架构。该网络以大地电磁观测数据（视电阻率和相位曲线）作为输入，直接输出地下二维电阻率分布图像，实现端到端的反演映射。网络设计借鉴了计算机视觉领域图像分割任务的成熟经验，同时针对大地电磁数据的特殊性质进行了适应性改进。
网络整体采用编码器-解码器结构。编码器部分负责从输入的大地电磁响应数据中提取多尺度特征表示，采用残差网络模块作为基本构建单元。残差连接的引入有效缓解了深层网络训练过程中的梯度消失问题，使得网络能够堆叠更多的卷积层以学习更加复杂的特征映射关系。在编码器的不同层级，设置了多尺度池化模块，通过并行使用不同尺寸的池化窗口来捕获不同空间范围的上下文信息。这种设计对于大地电磁数据尤为重要，因为地下不同深度的电性结构在观测数据中的响应特征具有不同的空间尺度。
解码器部分负责将编码器提取的抽象特征逐步恢复为与输入尺寸相匹配的电阻率分布图像。在上采样过程中，通过跳跃连接将编码器各层的特征图与解码器对应层进行融合，以保留更多的空间细节信息。特征融合采用通道拼接和卷积压缩的方式，使得网络能够综合利用浅层的精细空间信息和深层的高级语义信息。在网络的输出端，采用分类与分割相结合的策略，将电阻率值域离散化为若干个区间，将连续的回归问题转化为像素级分类问题，以提高训练的稳定性。
在损失函数设计方面，本研究采用Focal Loss与L2范数相结合的复合损失函数。Focal Loss是一种改进的交叉熵损失函数，通过动态调整样本权重来缓解类别不平衡问题。在大地电磁反演任务中，背景区域的像素数量通常远多于异常体区域，传统交叉熵损失函数会导致模型过度关注易分类的背景像素而忽视对异常体边界的精确刻画。L2范数项则用于约束预测电阻率值与真实值之间的数值误差，确保模型不仅能够正确分类各像素的电阻率区间，还能够给出较为准确的电阻率估计值。在优化策略方面，采用Adam与SGD相结合的组合优化方法，训练初期使用Adam实现快速收敛，后期切换为SGD进行精细调整。
（3）深度学习反演模型在实际数据中的应用验证
为验证所设计深度学习反演模型的实用性，本研究分别利用公开的非洲南部大地电磁数据和山东齐河地区的实测数据开展了反演实验研究。在模型部署之前，首先需要构建大规模的训练样本集。由于实测数据难以获取对应的真实地下电阻率模型，本研究采用理论模型正演的方式生成训练样本。通过设计包含不同类型异常体（高阻体、低阻体、组合异常体）、不同几何形态（块状、层状、透镜状）和不同埋深的大量理论模型，利用正演程序计算对应的大地电磁响应，形成输入输出配对的训练样本。
在非洲南部大地电磁数据的反演实验中，发现深度学习反演结果整体上能够反映该区域地下电性结构的基本特征，识别出主要的高阻和低阻异常体分布。然而，反演结果中也存在较多的细碎构造和噪点，影响了地质解释的清晰性。通过深入分析，认为造成这种现象的主要原因有两方面：一是分类问题固有的不连续性导致相邻像素可能被划分到不同的电阻率区间；二是实测数据的采集点位相对稀疏，经过插值处理后引入了额外的不确定性。
针对上述问题，本研究提出了改进方案。在网络结构层面，对多尺度池化层的卷积核尺寸进行调整优化，通过增大池化窗口来增强模型对噪声的鲁棒性和输出结果的空间连续性。在后处理层面，对网络输出的分类结果进行形态学滤波，消除孤立的噪点并平滑类别边界。改进后的反演结果在保持主要地质构造特征的同时，显著减少了冗余噪声的干扰。
在山东齐河低阻覆盖区实测数据的应用中，深度学习反演方法展现出相对于传统正则化反演的优势。该区域地表存在较厚的第四系低阻覆盖层，传统方法的反演结果受覆盖层屏蔽效应影响较大，难以清晰揭示深部目标层的电性特征。深度学习反演模型由于在训练阶段学习了大量包含覆盖层干扰的理论模型，具备了一定的抗干扰能力，反演结果对深部构造的成像效果更优。
import numpy as np import torch import torch.nn as nn import torch.nn.functional as F from scipy.sparse import lil_matrix, csr_matrix from scipy.sparse.linalg import spsolve class MT2DForward: def __init__(self, nx, nz, dx, dz, freq): self.nx = nx self.nz = nz self.dx = dx self.dz = dz self.freq = freq self.omega = 2 * np.pi * freq self.mu0 = 4 * np.pi * 1e-7 def build_mesh(self, resistivity): self.sigma = 1.0 / resistivity self.x = np.arange(self.nx + 1) * self.dx self.z = np.arange(self.nz + 1) * self.dz def assemble_stiffness_matrix(self): n_nodes = (self.nx + 1) * (self.nz + 1) K = lil_matrix((n_nodes, n_nodes), dtype=complex) for i in range(self.nz): for j in range(self.nx): elem_idx = i * self.nx + j sigma_e = self.sigma[i, j] nodes = [i * (self.nx + 1) + j, i * (self.nx + 1) + j + 1, (i + 1) * (self.nx + 1) + j, (i + 1) * (self.nx + 1) + j + 1] k_local = self.compute_element_matrix(sigma_e) for ii, ni in enumerate(nodes): for jj, nj in enumerate(nodes): K[ni, nj] += k_local[ii, jj] return csr_matrix(K) def compute_element_matrix(self, sigma): k = 1j * self.omega * self.mu0 * sigma dx, dz = self.dx, self.dz a = dz / (6 * dx) b = dx / (6 * dz) c = k * dx * dz / 36 Ke = np.array([[2*a + 2*b + 4*c, -2*a + b + 2*c, a - 2*b + 2*c, -a - b + c], [-2*a + b + 2*c, 2*a + 2*b + 4*c, -a - b + c, a - 2*b + 2*c], [a - 2*b + 2*c, -a - b + c, 2*a + 2*b + 4*c, -2*a + b + 2*c], [-a - b + c, a - 2*b + 2*c, -2*a + b + 2*c, 2*a + 2*b + 4*c]]) return Ke def apply_boundary_conditions(self, K, f): n_nodes = (self.nx + 1) * (self.nz + 1) for i in range(self.nz + 1): left = i * (self.nx + 1) right = i * (self.nx + 1) + self.nx K[left, :] = 0 K[left, left] = 1 K[right, :] = 0 K[right, right] = 1 return K, f def solve(self, resistivity): self.build_mesh(resistivity) K = self.assemble_stiffness_matrix() f = np.zeros((self.nx + 1) * (self.nz + 1), dtype=complex) K, f = self.apply_boundary_conditions(K, f) E = spsolve(K, f) return E.reshape(self.nz + 1, self.nx + 1) class ResidualBlock(nn.Module): def __init__(self, in_channels, out_channels, stride=1): super(ResidualBlock, self).__init__() self.conv1 = nn.Conv2d(in_channels, out_channels, 3, stride, 1) self.bn1 = nn.BatchNorm2d(out_channels) self.conv2 = nn.Conv2d(out_channels, out_channels, 3, 1, 1) self.bn2 = nn.BatchNorm2d(out_channels) self.shortcut = nn.Sequential() if stride != 1 or in_channels != out_channels: self.shortcut = nn.Sequential( nn.Conv2d(in_channels, out_channels, 1, stride), nn.BatchNorm2d(out_channels) ) def forward(self, x): out = F.relu(self.bn1(self.conv1(x))) out = self.bn2(self.conv2(out)) out += self.shortcut(x) return F.relu(out) class MultiScalePooling(nn.Module): def __init__(self, in_channels): super(MultiScalePooling, self).__init__() self.pool1 = nn.AdaptiveAvgPool2d(1) self.pool2 = nn.AdaptiveAvgPool2d(2) self.pool3 = nn.AdaptiveAvgPool2d(4) self.conv = nn.Conv2d(in_channels * 3, in_channels, 1) def forward(self, x): size = x.size()[2:] p1 = F.interpolate(self.pool1(x), size=size, mode='bilinear', align_corners=True) p2 = F.interpolate(self.pool2(x), size=size, mode='bilinear', align_corners=True) p3 = F.interpolate(self.pool3(x), size=size, mode='bilinear', align_corners=True) return self.conv(torch.cat([p1, p2, p3], dim=1)) class MTInversionNet(nn.Module): def __init__(self, in_channels=2, num_classes=10): super(MTInversionNet, self).__init__() self.encoder1 = nn.Sequential( nn.Conv2d(in_channels, 64, 3, padding=1), nn.BatchNorm2d(64), nn.ReLU(), ResidualBlock(64, 64) ) self.encoder2 = nn.Sequential( nn.MaxPool2d(2), ResidualBlock(64, 128), ResidualBlock(128, 128) ) self.encoder3 = nn.Sequential( nn.MaxPool2d(2), ResidualBlock(128, 256), ResidualBlock(256, 256) ) self.encoder4 = nn.Sequential( nn.MaxPool2d(2), ResidualBlock(256, 512), ResidualBlock(512, 512) ) self.msp = MultiScalePooling(512) self.decoder3 = nn.Sequential( nn.ConvTranspose2d(512, 256, 2, stride=2), ResidualBlock(512, 256) ) self.decoder2 = nn.Sequential( nn.ConvTranspose2d(256, 128, 2, stride=2), ResidualBlock(256, 128) ) self.decoder1 = nn.Sequential( nn.ConvTranspose2d(128, 64, 2, stride=2), ResidualBlock(128, 64) ) self.output = nn.Conv2d(64, num_classes, 1) def forward(self, x): e1 = self.encoder1(x) e2 = self.encoder2(e1) e3 = self.encoder3(e2) e4 = self.encoder4(e3) e4 = e4 + self.msp(e4) d3 = self.decoder3[0](e4) d3 = self.decoder3[1](torch.cat([d3, e3], dim=1)) d2 = self.decoder2[0](d3) d2 = self.decoder2[1](torch.cat([d2, e2], dim=1)) d1 = self.decoder1[0](d2) d1 = self.decoder1[1](torch.cat([d1, e1], dim=1)) return self.output(d1) class FocalLoss(nn.Module): def __init__(self, alpha=0.25, gamma=2.0): super(FocalLoss, self).__init__() self.alpha = alpha self.gamma = gamma def forward(self, pred, target): ce_loss = F.cross_entropy(pred, target, reduction='none') pt = torch.exp(-ce_loss) focal_loss = self.alpha * (1 - pt) ** self.gamma * ce_loss return focal_loss.mean() class CombinedLoss(nn.Module): def __init__(self, alpha=0.25, gamma=2.0, l2_weight=0.1): super(CombinedLoss, self).__init__() self.focal = FocalLoss(alpha, gamma) self.l2_weight = l2_weight def forward(self, pred, target, pred_value, true_value): focal_loss = self.focal(pred, target) l2_loss = F.mse_loss(pred_value, true_value) return focal_loss + self.l2_weight * l2_loss def train_mt_inversion(model, train_loader, epochs=100, lr=1e-3): device = torch.device('cuda' if torch.cuda.is_available() else 'cpu') model = model.to(device) optimizer = torch.optim.Adam(model.parameters(), lr=lr) scheduler = torch.optim.lr_scheduler.StepLR(optimizer, step_size=30, gamma=0.1) criterion = FocalLoss() for epoch in range(epochs如有问题，可以直接沟通
👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇👇
