帮我写出实验报告的python代码，要求在doc文档中，其中有五个实验要求，你用一些全局变量表示总体...
Created on: June 17, 2026
Answered using GPT-5.5 Thinking by Chat01
Question
重磁勘探学-实验指导书2026.docx
帮我写出实验报告的python代码，要求在doc文档中，其中有五个实验要求，你用一些全局变量表示总体条件，用五个函数完成五个实验要求
Python
pythonimport numpy as np MU0=4*np.pi*1e-7 T0_NT=50000.0 K=0.2 def magnetization_vector(inc_deg=60,a_deg=0): I=np.deg2rad(inc_deg); A=np.deg2rad(a_deg) M=K*(T0_NT*1e-9)/MU0 return np.array([M*np.cos(I)*np.cos(A),M*np.cos(I)*np.sin(A),M*np.sin(I)]) def unit(inc_deg=60,a_deg=0): I=np.deg2rad(inc_deg); A=np.deg2rad(a_deg) return np.array([np.cos(I)*np.cos(A),np.cos(I)*np.sin(A),np.sin(I)]) def sphere(x,y,h=30,radius=20,inc=60,a=0): rx,ry=np.broadcast_arrays(np.asarray(x,float),np.asarray(y,float)) rz=-h*np.ones_like(rx) r2=rx**2+ry**2+rz**2 r=np.sqrt(r2) m=magnetization_vector(inc,a)*(4/3*np.pi*radius**3) dot=m[0]*rx+m[1]*ry+m[2]*rz coeff=MU0/(4*np.pi) bx=coeff*(3*rx*dot/r**5 - m[0]/r**3)*1e9 by=coeff*(3*ry*dot/r**5 - m[1]/r**3)*1e9 bz=coeff*(3*rz*dot/r**5 - m[2]/r**3)*1e9 dt=bx*unit(inc,a)[0]+by*unit(inc,a)[1]+bz*unit(inc,a)[2] return bx,by,bz,dt def za_to_hax(za,dx,mult=1): F=np.fft.fft(za-np.mean(za)) k=2*np.pi*np.fft.fftfreq(len(za),dx) return np.real(np.fft.ifft(F*(mult*1j*np.sign(k)))) x=np.arange(-300,302,2) bx,by,za,dt=sphere(x,np.zeros_like(x)) for mult in [1,-1]: h=za_to_hax(za,2,mult) corr=np.corrcoef(h,bx)[0,1] scale=np.dot(h,bx)/np.dot(h,h) rmse=np.sqrt(np.mean((h*scale-bx)**2)) print(mult,corr,scale,rmse,np.max(np.abs(bx)))
Python
pythondef cyl(x,h=100,radius=30,inc=60,a=0): rx=np.asarray(x,float) rz=-h*np.ones_like(rx) r2=rx**2+rz**2 M=magnetization_vector(inc,a) mline=np.array([M[0],M[2]])*np.pi*radius**2 dot=mline[0]*rx+mline[1]*rz coeff=MU0/(2*np.pi) bx=coeff*(2*rx*dot/r2**2-mline[0]/r2)*1e9 bz=coeff*(2*rz*dot/r2**2-mline[1]/r2)*1e9 dt=bx*unit(inc,a)[0]+bz*unit(inc,a)[2] return bx,bz,dt x=np.arange(-1000,1000+2,2) bx,za,dt=cyl(x) for mult in [1,-1]: h=za_to_hax(za,2,mult) corr=np.corrcoef(h,bx)[0,1] scale=np.dot(h,bx)/np.dot(h,h) rmse=np.sqrt(np.mean((h-bx)**2)) print(mult,corr,scale,rmse,np.max(np.abs(bx)))
Thought
Python
pythonfrom pathlib import Path code = r'''# -*- coding: utf-8 -*- """ 《重磁勘探学原理》实验五：磁异常正演、异常处理及转换 说明： 1. 全局变量表示总体条件：地磁场、磁倾角、磁化率、球体/水平圆柱体几何参数等。 2. 五个 experiment_* 函数分别对应实验指导书中的五个实验要求。 3. 运行本文件后，图片会自动保存到 OUT_DIR 文件夹，可直接放入实验报告。 4. 坐标约定：x 为剖面方向，z 轴向下为正；Za 为向下垂直分量，Hax 为 x 方向水平分量。 """ import os from dataclasses import dataclass from typing import Dict, Tuple, List import numpy as np import matplotlib.pyplot as plt # ========================= # 一、全局变量：总体实验条件 # ========================= # 输出文件夹 OUT_DIR = "实验五_磁异常输出图件" # 物理常数 MU0 = 4.0 * np.pi * 1e-7 # 真空磁导率，H/m # 地磁场和磁化条件 T0_NT = 50000.0 # 地磁场强度，nT KAPPA = 0.2 # 磁化率，SI INC_DEG = 60.0 # 磁倾角 I，degree APRIME_DEG = 0.0 # 剖面与磁化强度水平投影夹角 A'，degree # 球体参数 SPHERE_H = 30.0 # 球心埋深，m SPHERE_R = 20.0 # 球半径，m # 水平圆柱体参数 CYL_H = 100.0 # 水平圆柱体中心埋深，m CYL_R = 30.0 # 水平圆柱体半径，m # 剖面与网格参数 PROFILE_X_MIN = -300.0 PROFILE_X_MAX = 300.0 PROFILE_DX = 2.0 GRID_X_MIN = -150.0 GRID_X_MAX = 150.0 GRID_Y_MIN = -150.0 GRID_Y_MAX = 150.0 GRID_STEP = 3.0 # 向上延拓高度 UPWARD_HEIGHTS = [20.0, 40.0, 80.0] # ========================= # 二、基础工具函数 # ========================= def setup_matplotlib() -> None: """设置中文字体和负号显示。""" plt.rcParams["font.sans-serif"] = [ "SimHei", "Microsoft YaHei", "Arial Unicode MS", "DejaVu Sans" ] plt.rcParams["axes.unicode_minus"] = False def ensure_out_dir() -> None: """创建输出文件夹。""" os.makedirs(OUT_DIR, exist_ok=True) def deg2rad(deg: float) -> float: return np.deg2rad(deg) def magnetization_amplitude(kappa: float = KAPPA, t0_nt: float = T0_NT) -> float: """ 由 J = k * T / μ0 计算磁化强度，单位 A/m。 注意 T 需由 nT 转成 T。 """ t0_tesla = t0_nt * 1e-9 return kappa * t0_tesla / MU0 def unit_vector_by_inc_aprime(inc_deg: float = INC_DEG, aprime_deg: float = APRIME_DEG) -> np.ndarray: """ 由磁倾角 I 和水平投影夹角 A' 得到方向单位向量。 z 轴向下为正，因此垂向分量为 sin(I)。 """ inc = deg2rad(inc_deg) a = deg2rad(aprime_deg) return np.array([ np.cos(inc) * np.cos(a), np.cos(inc) * np.sin(a), np.sin(inc) ], dtype=float) def magnetization_vector(inc_deg: float = INC_DEG, aprime_deg: float = APRIME_DEG, kappa: float = KAPPA, t0_nt: float = T0_NT) -> np.ndarray: """磁化强度矢量，单位 A/m。""" return magnetization_amplitude(kappa, t0_nt) * unit_vector_by_inc_aprime( inc_deg, aprime_deg ) def make_profile_x() -> np.ndarray: """生成主剖面 x 坐标。""" return np.arange(PROFILE_X_MIN, PROFILE_X_MAX + PROFILE_DX, PROFILE_DX) def make_grid_xy() -> Tuple[np.ndarray, np.ndarray]: """生成平面网格坐标。""" x = np.arange(GRID_X_MIN, GRID_X_MAX + GRID_STEP, GRID_STEP) y = np.arange(GRID_Y_MIN, GRID_Y_MAX + GRID_STEP, GRID_STEP) return np.meshgrid(x, y) def savefig(filename: str, dpi: int = 300) -> None: """保存当前图像。""" ensure_out_dir() path = os.path.join(OUT_DIR, filename) plt.tight_layout() plt.savefig(path, dpi=dpi) plt.close() print(f"已保存：{path}") def plot_profile(x: np.ndarray, curves: Dict[str, np.ndarray], title: str, filename: str, ylabel: str = "异常值 / nT") -> None: """绘制剖面曲线图。""" plt.figure(figsize=(8.5, 5.2)) for name, y in curves.items(): plt.plot(x, y, linewidth=1.6, label=name) plt.axhline(0, linewidth=0.8) plt.xlabel("x / m") plt.ylabel(ylabel) plt.title(title) plt.grid(True, linestyle="--", alpha=0.45) plt.legend() savefig(filename) def plot_contour(X: np.ndarray, Y: np.ndarray, Z: np.ndarray, title: str, filename: str, cbar_label: str = "异常值 / nT") -> None: """绘制平面等值线图。""" plt.figure(figsize=(7.2, 6.0)) cf = plt.contourf(X, Y, Z, levels=30) cs = plt.contour(X, Y, Z, levels=12, linewidths=0.5) plt.clabel(cs, inline=True, fontsize=7, fmt="%.0f") plt.xlabel("x / m") plt.ylabel("y / m") plt.title(title) plt.axis("equal") cbar = plt.colorbar(cf) cbar.set_label(cbar_label) savefig(filename) def plot_surface(X: np.ndarray, Y: np.ndarray, Z: np.ndarray, title: str, filename: str, zlabel: str = "异常值 / nT") -> None: """绘制三维曲面图。""" fig = plt.figure(figsize=(8.2, 6.2)) ax = fig.add_subplot(111, projection="3d") surf = ax.plot_surface(X, Y, Z, cmap="viridis", linewidth=0, antialiased=True) ax.set_xlabel("x / m") ax.set_ylabel("y / m") ax.set_zlabel(zlabel) ax.set_title(title) fig.colorbar(surf, ax=ax, shrink=0.70, pad=0.10, label=zlabel) savefig(filename) def normalize_for_compare(y: np.ndarray) -> np.ndarray: """归一化到最大绝对值为 1，便于形态对比。""" y = np.asarray(y, dtype=float) m = np.nanmax(np.abs(y)) if m == 0: return y return y / m # ========================= # 三、正演计算公式 # ========================= def sphere_field(x, y, h: float = SPHERE_H, radius: float = SPHERE_R, inc_deg: float = INC_DEG, aprime_deg: float = APRIME_DEG, kappa: float = KAPPA) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: """ 均匀磁化球体磁异常正演。 返回： Bx, By, Bz, deltaT，单位均为 nT。 其中 Bx 即 Hax，Bz 即 Za，deltaT 为总场异常 ΔT 的投影近似值。 球体外部磁场等效于位于球心的磁偶极子： B = μ0/(4π) * [3r(m·r)/r^5 - m/r^3] """ rx, ry = np.broadcast_arrays(np.asarray(x, dtype=float), np.asarray(y, dtype=float)) rz = -h * np.ones_like(rx) # 观测面 z=0，球心 z=h，故源到测点矢量的 z 分量为 -h r2 = rx**2 + ry**2 + rz**2 r = np.sqrt(r2) volume = 4.0 / 3.0 * np.pi * radius**3 m = magnetization_vector(inc_deg, aprime_deg, kappa) * volume mdotr = m[0] * rx + m[1] * ry + m[2] * rz coeff = MU0 / (4.0 * np.pi) bx = coeff * (3.0 * rx * mdotr / r**5 - m[0] / r**3) * 1e9 by = coeff * (3.0 * ry * mdotr / r**5 - m[1] / r**3) * 1e9 bz = coeff * (3.0 * rz * mdotr / r**5 - m[2] / r**3) * 1e9 fhat = unit_vector_by_inc_aprime(inc_deg, aprime_deg) delta_t = bx * fhat[0] + by * fhat[1] + bz * fhat[2] return bx, by, bz, delta_t def cylinder_field(x, h: float = CYL_H, radius: float = CYL_R, inc_deg: float = INC_DEG, aprime_deg: float = APRIME_DEG, kappa: float = KAPPA) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: """ 无限长水平圆柱体磁异常正演，剖面垂直于圆柱轴向。 返回： Bx, Bz, deltaT，单位均为 nT。 其中 Bx 即 Hax，Bz 即 Za。 圆柱体外部磁场可用二维线偶极子表示： B = μ0/(2π) * [2r(m_l·r)/r^4 - m_l/r^2] m_l = πR^2 M，为单位长度磁矩。 """ rx = np.asarray(x, dtype=float) rz = -h * np.ones_like(rx) r2 = rx**2 + rz**2 m3d = magnetization_vector(inc_deg, aprime_deg, kappa) # 水平圆柱体沿 y 方向无限延伸，外部异常主要由 x-z 平面分量产生 mline_x = np.pi * radius**2 * m3d[0] mline_z = np.pi * radius**2 * m3d[2] mdotr = mline_x * rx + mline_z * rz coeff = MU0 / (2.0 * np.pi) bx = coeff * (2.0 * rx * mdotr / r2**2 - mline_x / r2) * 1e9 bz = coeff * (2.0 * rz * mdotr / r2**2 - mline_z / r2) * 1e9 fhat = unit_vector_by_inc_aprime(inc_deg, aprime_deg) delta_t = bx * fhat[0] + bz * fhat[2] return bx, bz, delta_t # ========================= # 四、异常处理公式 # ========================= def upward_continue_profile(data: np.ndarray, dx: float, height: float) -> np.ndarray: """ 一维剖面异常向上延拓，FFT 形式： F_up(k) = F(k) * exp(-|k| * height) 注意： 1. FFT 默认存在周期延拓，边部会有一定影响； 2. 分析时建议重点看剖面中部。 """ data = np.asarray(data, dtype=float) k = 2.0 * np.pi * np.fft.fftfreq(data.size, d=dx) f = np.fft.fft(data) f_up = f * np.exp(-np.abs(k) * height) return np.real(np.fft.ifft(f_up)) def upward_continue_grid(data: np.ndarray, dx: float, dy: float, height: float) -> np.ndarray: """ 二维平面异常向上延拓，FFT 形式： F_up(kx, ky) = F(kx, ky) * exp(-sqrt(kx^2 + ky^2) * height) """ data = np.asarray(data, dtype=float) ny, nx = data.shape kx = 2.0 * np.pi * np.fft.fftfreq(nx, d=dx) ky = 2.0 * np.pi * np.fft.fftfreq(ny, d=dy) KX, KY = np.meshgrid(kx, ky) kk = np.sqrt(KX**2 + KY**2) f = np.fft.fft2(data) f_up = f * np.exp(-kk * height) return np.real(np.fft.ifft2(f_up)) def za_to_hax_fft(za: np.ndarray, dx: float) -> np.ndarray: """ 由 Za 转换为 Hax 的频率域算法。 对二维位场剖面，水平导数分量与垂向导数分量近似满足希尔伯特变换关系： Hax(k) = i * sign(k) * Za(k) 该公式最适合二维体异常，例如水平圆柱体剖面。 """ za = np.asarray(za, dtype=float) za0 = za - np.mean(za) k = 2.0 * np.pi * np.fft.fftfreq(za0.size, d=dx) f = np.fft.fft(za0) hax = np.real(np.fft.ifft(f * 1j * np.sign(k))) return hax # ========================= # 五、五个实验函数 # ========================= def experiment1_sphere_forward() -> Dict[str, np.ndarray]: """ 实验要求 1： 正演计算球体磁场 Za、Hax、ΔT； 绘制平面等值线图、曲面图、主剖面异常图。 """ print("\n========== 实验 1：球体磁异常正演 ==========") X, Y = make_grid_xy() bx, by, za, dt = sphere_field(X, Y) field_map = { "Za": za, "Hax": bx, "DeltaT": dt } chinese_name = { "Za": "球体 Za 垂直磁异常", "Hax": "球体 Hax 水平磁异常", "DeltaT": "球体 ΔT 总场异常" } for key, value in field_map.items(): plot_contour( X, Y, value, title=f"{chinese_name[key]}平面等值线图", filename=f"实验1_{key}_平面等值线图.png" ) plot_surface( X, Y, value, title=f"{chinese_name[key]}三维曲面图", filename=f"实验1_{key}_三维曲面图.png" ) x = make_profile_x() y0 = np.zeros_like(x) bx_p, by_p, za_p, dt_p = sphere_field(x, y0) plot_profile( x, {"Za": za_p, "Hax": bx_p, "ΔT": dt_p}, title="实验1：球体主剖面磁异常", filename="实验1_球体主剖面_Za_Hax_DeltaT.png" ) return { "X": X, "Y": Y, "Za_grid": za, "Hax_grid": bx, "DeltaT_grid": dt, "x": x, "Za_profile": za_p, "Hax_profile": bx_p, "DeltaT_profile": dt_p } def experiment2_cylinder_forward() -> Dict[str, np.ndarray]: """ 实验要求 2： 正演计算水平圆柱体磁场 Za、Hax、ΔT； 绘制主剖面异常图。 """ print("\n========== 实验 2：水平圆柱体磁异常正演 ==========") x = make_profile_x() bx, za, dt = cylinder_field(x) plot_profile( x, {"Za": za, "Hax": bx, "ΔT": dt}, title="实验2：水平圆柱体主剖面磁异常", filename="实验2_水平圆柱体主剖面_Za_Hax_DeltaT.png" ) # 分量单独图，便于报告排版 plot_profile(x, {"Za": za}, "实验2：水平圆柱体 Za 主剖面", "实验2_圆柱体_Za主剖面.png") plot_profile(x, {"Hax": bx}, "实验2：水平圆柱体 Hax 主剖面", "实验2_圆柱体_Hax主剖面.png") plot_profile(x, {"ΔT": dt}, "实验2：水平圆柱体 ΔT 主剖面", "实验2_圆柱体_DeltaT主剖面.png") return { "x": x, "Za_profile": za, "Hax_profile": bx, "DeltaT_profile": dt } def experiment3_parameter_analysis() -> Dict[str, Dict[str, np.ndarray]]: """ 实验要求 3： 改变球体与水平圆柱体的几何参数、磁化强度方向 I， 观察主剖面 Za 的变化，分析磁化方向与计算剖面对磁性体磁场特征的影响。 """ print("\n========== 实验 3：参数变化对 Za 的影响 ==========") x = make_profile_x() y0 = np.zeros_like(x) results = {} # 1）球体埋深变化 sphere_depth_cases = [20.0, 30.0, 50.0, 80.0] curves = {} for h in sphere_depth_cases: _, _, za, _ = sphere_field(x, y0, h=h, radius=SPHERE_R, inc_deg=INC_DEG) curves[f"h={h:.0f} m"] = za plot_profile( x, curves, "实验3：球体埋深变化对 Za 的影响", "实验3_球体埋深变化_Za.png" ) results["sphere_depth"] = curves # 2）球体半径变化 sphere_radius_cases = [10.0, 20.0, 30.0, 40.0] curves = {} for r in sphere_radius_cases: _, _, za, _ = sphere_field(x, y0, h=SPHERE_H, radius=r, inc_deg=INC_DEG) curves[f"R={r:.0f} m"] = za plot_profile( x, curves, "实验3：球体半径变化对 Za 的影响", "实验3_球体半径变化_Za.png" ) results["sphere_radius"] = curves # 3）球体磁倾角变化 inc_cases = [0.0, 30.0, 60.0, 90.0] curves = {} for inc in inc_cases: _, _, za, _ = sphere_field(x, y0, h=SPHERE_H, radius=SPHERE_R, inc_deg=inc) curves[f"I={inc:.0f}°"] = za plot_profile( x, curves, "实验3：球体磁倾角 I 变化对 Za 的影响", "实验3_球体磁倾角变化_Za.png" ) results["sphere_inc"] = curves # 4）圆柱体埋深变化 cyl_depth_cases = [60.0, 100.0, 150.0, 200.0] curves = {} for h in cyl_depth_cases: _, za, _ = cylinder_field(x, h=h, radius=CYL_R, inc_deg=INC_DEG) curves[f"h={h:.0f} m"] = za plot_profile( x, curves, "实验3：水平圆柱体埋深变化对 Za 的影响", "实验3_圆柱体埋深变化_Za.png" ) results["cyl_depth"] = curves # 5）圆柱体半径变化 cyl_radius_cases = [15.0, 30.0, 45.0, 60.0] curves = {} for r in cyl_radius_cases: _, za, _ = cylinder_field(x, h=CYL_H, radius=r, inc_deg=INC_DEG) curves[f"R={r:.0f} m"] = za plot_profile( x, curves, "实验3：水平圆柱体半径变化对 Za 的影响", "实验3_圆柱体半径变化_Za.png" ) results["cyl_radius"] = curves # 6）圆柱体磁倾角变化 curves = {} for inc in inc_cases: _, za, _ = cylinder_field(x, h=CYL_H, radius=CYL_R, inc_deg=inc) curves[f"I={inc:.0f}°"] = za plot_profile( x, curves, "实验3：水平圆柱体磁倾角 I 变化对 Za 的影响", "实验3_圆柱体磁倾角变化_Za.png" ) results["cyl_inc"] = curves # 7）A' 变化，体现计算剖面与磁化水平投影夹角影响 aprime_cases = [0.0, 30.0, 60.0, 90.0] curves = {} for a in aprime_cases: _, _, za, _ = sphere_field(x, y0, h=SPHERE_H, radius=SPHERE_R, inc_deg=INC_DEG, aprime_deg=a) curves[f"A'={a:.0f}°"] = za plot_profile( x, curves, "实验3：球体 A' 变化对 Za 的影响", "实验3_球体Aprime变化_Za.png" ) results["sphere_aprime"] = curves curves = {} for a in aprime_cases: _, za, _ = cylinder_field(x, h=CYL_H, radius=CYL_R, inc_deg=INC_DEG, aprime_deg=a) curves[f"A'={a:.0f}°"] = za plot_profile( x, curves, "实验3：水平圆柱体 A' 变化对 Za 的影响", "实验3_圆柱体Aprime变化_Za.png" ) results["cyl_aprime"] = curves print("实验3分析提示：") print("1. 埋深增大时，异常幅值减小，异常曲线变宽。") print("2. 半径增大时，异常幅值增大；球体约受体积 R^3 控制，圆柱体约受截面积 R^2 控制。") print("3. I 和 A' 改变会改变正负异常分布、极值位置和对称性。") return results def sphere_vzz_kernel(x, y, h: float = SPHERE_H) -> np.ndarray: """ 球体重力位二阶垂向导数核函数形态： ∂²(1/r)/∂z² = (3z² - r²) / r^5 只用于形态对比，不含万有引力常数和密度。 """ rx, ry = np.broadcast_arrays(np.asarray(x, dtype=float), np.asarray(y, dtype=float)) rz = -h * np.ones_like(rx) r2 = rx**2 + ry**2 + rz**2 return (3.0 * rz**2 - r2) / (r2 ** 2.5) def cylinder_vzz_kernel(x, h: float = CYL_H) -> np.ndarray: """ 水平圆柱体重力位二阶垂向导数核函数形态。 对二维 ln(r) 位函数，∂²ln(r)/∂z² = (x² - z²)/r^4。 磁场符号可能因位函数和坐标正方向约定不同而相差一个负号。 """ rx = np.asarray(x, dtype=float) rz = -h * np.ones_like(rx) r2 = rx**2 + rz**2 return (rx**2 - rz**2) / r2**2 def experiment4_vertical_superposition() -> Dict[str, np.ndarray]: """ 实验要求 4： 垂直磁化时，分析球体磁场与水平圆柱体磁场 Za、Hax、ΔT 的叠加异常； 并与重力位二阶偏导数形态进行对比。 """ print("\n========== 实验 4：垂直磁化叠加异常与重力位二阶导数对比 ==========") x = make_profile_x() y0 = np.zeros_like(x) # 垂直磁化：I = 90° bx_s, _, za_s, dt_s = sphere_field(x, y0, inc_deg=90.0) bx_c, za_c, dt_c = cylinder_field(x, inc_deg=90.0) bx_sum = bx_s + bx_c za_sum = za_s + za_c dt_sum = dt_s + dt_c plot_profile( x, {"球体 Za": za_s, "圆柱体 Za": za_c, "叠加 Za": za_sum}, "实验4：垂直磁化 Za 叠加异常", "实验4_垂直磁化_Za叠加异常.png" ) plot_profile( x, {"球体 Hax": bx_s, "圆柱体 Hax": bx_c, "叠加 Hax": bx_sum}, "实验4：垂直磁化 Hax 叠加异常", "实验4_垂直磁化_Hax叠加异常.png" ) plot_profile( x, {"球体 ΔT": dt_s, "圆柱体 ΔT": dt_c, "叠加 ΔT": dt_sum}, "实验4：垂直磁化 ΔT 叠加异常", "实验4_垂直磁化_DeltaT叠加异常.png" ) # 与重力位二阶垂向导数核函数形态对比：归一化后比较 vzz_s = sphere_vzz_kernel(x, y0) vzz_c = cylinder_vzz_kernel(x) vzz_sum = normalize_for_compare(vzz_s) - normalize_for_compare(vzz_c) # 圆柱体部分用负号，是为了与当前磁场正方向约定下的 Za 形态更一致 plot_profile( x, { "归一化叠加 Za": normalize_for_compare(za_sum), "归一化二阶导数形态": normalize_for_compare(vzz_sum) }, "实验4：垂直磁化 Za 与重力位二阶导数形态对比", "实验4_Za与重力位二阶导数形态对比.png", ylabel="归一化幅值" ) print("实验4分析提示：") print("垂直磁化时，磁异常与重力位二阶偏导数具有相同或相近的核函数形态，") print("这体现了重磁位场之间的泊松公式联系；不同分量之间主要表现为比例系数、方向和符号约定差异。") return { "x": x, "Za_sphere": za_s, "Za_cylinder": za_c, "Za_sum": za_sum, "Hax_sum": bx_sum, "DeltaT_sum": dt_sum, "Vzz_shape_sum": vzz_sum } def experiment5_processing() -> Dict[str, np.ndarray]: """ 实验要求 5： 对磁异常进行向上延拓和 Za 求 Hax 的异常处理。 """ print("\n========== 实验 5：磁异常向上延拓与 Za 转 Hax ==========") # A. 球体平面 Za 做二维向上延拓，给出等值线图 X, Y = make_grid_xy() _, _, za_grid, _ = sphere_field(X, Y) for height in UPWARD_HEIGHTS: za_up = upward_continue_grid(za_grid, GRID_STEP, GRID_STEP, height) plot_contour( X, Y, za_up, title=f"实验5：球体 Za 向上延拓 {height:.0f} m 平面等值线图", filename=f"实验5_球体Za向上延拓_{height:.0f}m_等值线图.png" ) # B. 水平圆柱体剖面 Za 做一维向上延拓 x = make_profile_x() hax_true, za, dt = cylinder_field(x) curves = {"原始 Za": za} for height in UPWARD_HEIGHTS: curves[f"上延 {height:.0f} m"] = upward_continue_profile(za, PROFILE_DX, height) plot_profile( x, curves, "实验5：水平圆柱体 Za 剖面向上延拓", "实验5_圆柱体Za剖面向上延拓.png" ) # C. 由 Za 转换为 Hax，并与正演 Hax 对比 hax_from_za = za_to_hax_fft(za, PROFILE_DX) plot_profile( x, { "正演 Hax": hax_true, "Za 转换 Hax": hax_from_za }, "实验5：Za 求 Hax 处理结果与正演结果对比", "实验5_Za求Hax_与正演对比.png" ) # 误差指标，仅作程序检查 rmse = float(np.sqrt(np.mean((hax_from_za - hax_true) ** 2))) corr = float(np.corrcoef(hax_from_za, hax_true)[0, 1]) print(f"Za 转 Hax 与正演 Hax 对比：RMSE = {rmse:.3f} nT，相关系数 = {corr:.5f}") print("注意：Za->Hax 的 FFT 转换适用于二维位场剖面，边部可能有周期延拓误差。") return { "x": x, "Za": za, "Hax_true": hax_true, "Hax_from_Za": hax_from_za, "DeltaT": dt } # ========================= # 六、主程序 # ========================= def main() -> None: setup_matplotlib() ensure_out_dir() print("开始运行《重磁勘探学原理》实验五 Python 程序") print(f"磁化强度 J = {magnetization_amplitude():.4f} A/m") print(f"输出文件夹：{OUT_DIR}") exp1 = experiment1_sphere_forward() exp2 = experiment2_cylinder_forward() exp3 = experiment3_parameter_analysis() exp4 = experiment4_vertical_superposition() exp5 = experiment5_processing() print("\n全部实验图件已生成。") print("建议在实验报告中按：实验目的—原理—代码—结果图—结果分析 的结构整理。") if __name__ == "__main__": main() ''' path = Path("/mnt/data/magnetic_experiment5_report_code.py") path.write_text(code, encoding="utf-8") # 语法检查 import py_compile py_compile.compile(str(path), doraise=True) print(f"已生成并通过语法检查：{path}")
Answer
已按实验指导书里的五个实验要求整理成一个完整 Python 脚本：用全局变量表示总体实验条件，并用五个 experiment_* 函数分别完成五个实验。指导书要求包括球体正演、水平圆柱体正演、参数变化分析、垂直磁化叠加异常、向上延拓与 Za 求 Hax 处理等内容。(重磁勘探学-实验指导书2026.docx)
你可以下载完整代码：
运行方式：
bashpython magnetic_experiment5_report_code.py
运行后会自动生成文件夹：
text实验五_磁异常输出图件
里面包含实验报告可用的平面等值线图、三维曲面图和主剖面异常图。
代码结构如下：
python# 全局变量：T0_NT、INC_DEG、APRIME_DEG、KAPPA、SPHERE_H、SPHERE_R、CYL_H、CYL_R 等 def experiment1_sphere_forward(): # 球体 Za、Hax、ΔT 正演 # 输出平面等值线图、曲面图、主剖面图 def experiment2_cylinder_forward(): # 水平圆柱体 Za、Hax、ΔT 正演 # 输出主剖面图 def experiment3_parameter_analysis(): # 改变埋深、半径、磁倾角 I、A' # 分析 Za 主剖面变化 def experiment4_vertical_superposition(): # 垂直磁化 I=90° # 球体与圆柱体 Za、Hax、ΔT 叠加 # 与重力位二阶导数形态对比 def experiment5_processing(): # 磁异常向上延拓 # Za 转 Hax
代码里已经写了中文注释，图名、坐标轴、单位和色标也都按实验报告需要设置好了。