768 lines
1.0 MiB
Plaintext
768 lines
1.0 MiB
Plaintext
|
{
|
|||
|
|
"cells": [
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"# 球体重磁异常正演计算\n",
|
|||
|
|
"\n",
|
|||
|
|
"本notebook实现了球体的重力异常和磁异常正演计算，包含理论公式、数值计算和可视化展示。\n",
|
|||
|
|
"\n",
|
|||
|
|
"## 1. 理论基础\n",
|
|||
|
|
"\n",
|
|||
|
|
"### 1.1 球体重力异常\n",
|
|||
|
|
"对于半径为 $R$、密度为 $\\rho$、中心埋深为 $h$ 的球体，其在地面 $(x,y)$ 处的重力异常为：\n",
|
|||
|
|
"\n",
|
|||
|
|
"$$\\Delta g = \\frac{4\\pi G \\rho R^3}{3} \\cdot \\frac{h}{(x^2 + y^2 + h^2)^{3/2}}$$\n",
|
|||
|
|
"\n",
|
|||
|
|
"其中：\n",
|
|||
|
|
"- $G = 6.67 \\times 10^{-11} \\text{ N·m}^2/\\text{kg}^2$ 为万有引力常数\n",
|
|||
|
|
"- $\\rho$ 为球体密度\n",
|
|||
|
|
"- $R$ 为球体半径\n",
|
|||
|
|
"- $h$ 为球体中心埋深\n",
|
|||
|
|
"\n",
|
|||
|
|
"### 1.2 球体磁异常\n",
|
|||
|
|
"对于磁化强度为 $M$、磁化倾角为 $I$、磁化偏角为 $D$ 的球体，其在地面 $(x,y)$ 处的磁异常为：\n",
|
|||
|
|
"\n",
|
|||
|
|
"$$\\Delta T = \\frac{\\mu_0 M V}{4\\pi} \\cdot \\frac{1}{(x^2 + y^2 + h^2)^{5/2}} \\cdot [(2h^2 - x^2 - y^2)\\sin I + 3hx\\cos I \\cos D + 3hy\\cos I \\sin D]$$\n",
|
|||
|
|
"\n",
|
|||
|
|
"其中：\n",
|
|||
|
|
"- $\\mu_0 = 4\\pi \\times 10^{-7} \\text{ H/m}$ 为真空磁导率\n",
|
|||
|
|
"- $M$ 为磁化强度\n",
|
|||
|
|
"- $V = \\frac{4}{3}\\pi R^3$ 为球体体积\n",
|
|||
|
|
"- $I$ 为磁化倾角\n",
|
|||
|
|
"- $D$ 为磁化偏角"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 1,
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [],
|
|||
|
|
"source": [
|
|||
|
|
"# 导入必要的库\n",
|
|||
|
|
"import numpy as np\n",
|
|||
|
|
"import matplotlib.pyplot as plt\n",
|
|||
|
|
"from mpl_toolkits.mplot3d import Axes3D\n",
|
|||
|
|
"import ipywidgets as widgets\n",
|
|||
|
|
"from ipywidgets import interact, FloatSlider, IntSlider\n",
|
|||
|
|
"import warnings\n",
|
|||
|
|
"warnings.filterwarnings('ignore')\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 设置中文字体\n",
|
|||
|
|
"plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']\n",
|
|||
|
|
"plt.rcParams['axes.unicode_minus'] = False\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 物理常数\n",
|
|||
|
|
"G = 6.67e-11 # 万有引力常数 (N·m²/kg²)\n",
|
|||
|
|
"mu_0 = 4 * np.pi * 1e-7 # 真空磁导率 (H/m)"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"## 2. 正演计算函数"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 2,
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [],
|
|||
|
|
"source": [
|
|||
|
|
"def gravity_anomaly_sphere(x, y, R, rho, h):\n",
|
|||
|
|
" \"\"\"\n",
|
|||
|
|
" 计算球体重力异常\n",
|
|||
|
|
" \n",
|
|||
|
|
" 参数:\n",
|
|||
|
|
" x, y: 观测点坐标 (m)\n",
|
|||
|
|
" R: 球体半径 (m)\n",
|
|||
|
|
" rho: 球体密度 (kg/m³)\n",
|
|||
|
|
" h: 球体中心埋深 (m)\n",
|
|||
|
|
" \n",
|
|||
|
|
" 返回:\n",
|
|||
|
|
" 重力异常值 (mGal)\n",
|
|||
|
|
" \"\"\"\n",
|
|||
|
|
" r_squared = x**2 + y**2 + h**2\n",
|
|||
|
|
" r = np.sqrt(r_squared)\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 重力异常公式 (单位: m/s²)\n",
|
|||
|
|
" delta_g = (4 * np.pi * G * rho * R**3 / 3) * h / (r_squared)**(3/2)\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 转换为 mGal (1 mGal = 1e-5 m/s²)\n",
|
|||
|
|
" return delta_g * 1e5\n",
|
|||
|
|
"\n",
|
|||
|
|
"def magnetic_anomaly_sphere(x, y, R, M, h, I, D):\n",
|
|||
|
|
" \"\"\"\n",
|
|||
|
|
" 计算球体磁异常\n",
|
|||
|
|
" \n",
|
|||
|
|
" 参数:\n",
|
|||
|
|
" x, y: 观测点坐标 (m)\n",
|
|||
|
|
" R: 球体半径 (m)\n",
|
|||
|
|
" M: 磁化强度 (A/m)\n",
|
|||
|
|
" h: 球体中心埋深 (m)\n",
|
|||
|
|
" I: 磁化倾角 (度)\n",
|
|||
|
|
" D: 磁化偏角 (度)\n",
|
|||
|
|
" \n",
|
|||
|
|
" 返回:\n",
|
|||
|
|
" 磁异常值 (nT)\n",
|
|||
|
|
" \"\"\"\n",
|
|||
|
|
" # 角度转换为弧度\n",
|
|||
|
|
" I_rad = np.radians(I)\n",
|
|||
|
|
" D_rad = np.radians(D)\n",
|
|||
|
|
" \n",
|
|||
|
|
" r_squared = x**2 + y**2 + h**2\n",
|
|||
|
|
" r = np.sqrt(r_squared)\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 球体体积\n",
|
|||
|
|
" V = 4 * np.pi * R**3 / 3\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 磁异常公式\n",
|
|||
|
|
" factor = mu_0 * M * V / (4 * np.pi * r_squared**(5/2))\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 磁异常分量\n",
|
|||
|
|
" delta_T = factor * ((2*h**2 - x**2 - y**2) * np.sin(I_rad) + \n",
|
|||
|
|
" 3*h*x * np.cos(I_rad) * np.cos(D_rad) + \n",
|
|||
|
|
" 3*h*y * np.cos(I_rad) * np.sin(D_rad))\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 转换为 nT (1 nT = 1e-9 T)\n",
|
|||
|
|
" return delta_T * 1e9"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"## 3. 可视化函数"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 3,
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [],
|
|||
|
|
"source": [
|
|||
|
|
"def plot_gravity_anomaly(X, Y, gravity_data, R, rho, h):\n",
|
|||
|
|
" \"\"\"绘制重力异常等值线图和剖面图\"\"\"\n",
|
|||
|
|
" fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 等值线图\n",
|
|||
|
|
" contour = ax1.contour(X, Y, gravity_data, levels=20, colors='black', linewidths=0.5)\n",
|
|||
|
|
" ax1.clabel(contour, inline=True, fontsize=8)\n",
|
|||
|
|
" im1 = ax1.contourf(X, Y, gravity_data, levels=20, cmap='RdBu_r', alpha=0.8)\n",
|
|||
|
|
" ax1.set_xlabel('X (m)')\n",
|
|||
|
|
" ax1.set_ylabel('Y (m)')\n",
|
|||
|
|
" ax1.set_title(f'Gravity Anomaly (2D)\\nR={R}m, ρ={rho}kg/m³, h={h}m')\n",
|
|||
|
|
" plt.colorbar(im1, ax=ax1, label='Gravity Anomaly (mGal)')\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 剖面图 (y=0)\n",
|
|||
|
|
" profile_y = 0\n",
|
|||
|
|
" profile_idx = np.argmin(np.abs(Y[:, 0] - profile_y))\n",
|
|||
|
|
" profile_x = X[profile_idx, :]\n",
|
|||
|
|
" profile_gravity = gravity_data[profile_idx, :]\n",
|
|||
|
|
" \n",
|
|||
|
|
" ax2.plot(profile_x, profile_gravity, 'b-', linewidth=2, label='Gravity Anomaly (1D)')\n",
|
|||
|
|
" ax2.set_xlabel('X (m)')\n",
|
|||
|
|
" ax2.set_ylabel('Gravity Anomaly (mGal)')\n",
|
|||
|
|
" ax2.set_title(f'Gravity Anomaly (1D) (Y={profile_y}m)')\n",
|
|||
|
|
" ax2.grid(True, alpha=0.3)\n",
|
|||
|
|
" ax2.legend()\n",
|
|||
|
|
" \n",
|
|||
|
|
" plt.tight_layout()\n",
|
|||
|
|
" plt.show()\n",
|
|||
|
|
"\n",
|
|||
|
|
"def plot_magnetic_anomaly(X, Y, magnetic_data, R, M, h, I, D):\n",
|
|||
|
|
" \"\"\"绘制磁异常等值线图和剖面图\"\"\"\n",
|
|||
|
|
" fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 等值线图\n",
|
|||
|
|
" contour = ax1.contour(X, Y, magnetic_data, levels=20, colors='black', linewidths=0.5)\n",
|
|||
|
|
" ax1.clabel(contour, inline=True, fontsize=8)\n",
|
|||
|
|
" im1 = ax1.contourf(X, Y, magnetic_data, levels=20, cmap='RdBu_r', alpha=0.8)\n",
|
|||
|
|
" ax1.set_xlabel('X (m)')\n",
|
|||
|
|
" ax1.set_ylabel('Y (m)')\n",
|
|||
|
|
" ax1.set_title(f'磁异常等值线图\\nR={R}m, M={M}A/m, h={h}m, I={I}°, D={D}°')\n",
|
|||
|
|
" plt.colorbar(im1, ax=ax1, label='磁异常 (nT)')\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 剖面图 (y=0)\n",
|
|||
|
|
" profile_y = 0\n",
|
|||
|
|
" profile_idx = np.argmin(np.abs(Y[:, 0] - profile_y))\n",
|
|||
|
|
" profile_x = X[profile_idx, :]\n",
|
|||
|
|
" profile_magnetic = magnetic_data[profile_idx, :]\n",
|
|||
|
|
" \n",
|
|||
|
|
" ax2.plot(profile_x, profile_magnetic, 'r-', linewidth=2, label='磁异常')\n",
|
|||
|
|
" ax2.set_xlabel('X (m)')\n",
|
|||
|
|
" ax2.set_ylabel('磁异常 (nT)')\n",
|
|||
|
|
" ax2.set_title(f'磁异常剖面图 (Y={profile_y}m)')\n",
|
|||
|
|
" ax2.grid(True, alpha=0.3)\n",
|
|||
|
|
" ax2.legend()\n",
|
|||
|
|
" \n",
|
|||
|
|
" plt.tight_layout()\n",
|
|||
|
|
" plt.show()\n",
|
|||
|
|
"\n",
|
|||
|
|
"def plot_3d_anomaly(X, Y, data, title, zlabel):\n",
|
|||
|
|
" \"\"\"绘制3D异常图\"\"\"\n",
|
|||
|
|
" fig = plt.figure(figsize=(12, 8))\n",
|
|||
|
|
" ax = fig.add_subplot(111, projection='3d')\n",
|
|||
|
|
" \n",
|
|||
|
|
" surf = ax.plot_surface(X, Y, data, cmap='RdBu_r', alpha=0.8)\n",
|
|||
|
|
" ax.set_xlabel('X (m)')\n",
|
|||
|
|
" ax.set_ylabel('Y (m)')\n",
|
|||
|
|
" ax.set_zlabel(zlabel)\n",
|
|||
|
|
" ax.set_title(title)\n",
|
|||
|
|
" \n",
|
|||
|
|
" plt.colorbar(surf, ax=ax, label=zlabel)\n",
|
|||
|
|
" plt.show()"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"## 4. 交互式参数调整"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 4,
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [],
|
|||
|
|
"source": [
|
|||
|
|
"def interactive_gravity_anomaly():\n",
|
|||
|
|
" \"\"\"交互式重力异常计算\"\"\"\n",
|
|||
|
|
" \n",
|
|||
|
|
" @interact\n",
|
|||
|
|
" def gravity_plot(R=FloatSlider(min=10, max=200, step=10, value=50, description='半径 (m)'),\n",
|
|||
|
|
" rho=FloatSlider(min=1000, max=5000, step=100, value=3000, description='密度 (kg/m³)'),\n",
|
|||
|
|
" h=FloatSlider(min=20, max=200, step=10, value=100, description='埋深 (m)'),\n",
|
|||
|
|
" x_range=FloatSlider(min=100, max=1000, step=50, value=400, description='X范围 (m)'),\n",
|
|||
|
|
" y_range=FloatSlider(min=100, max=1000, step=50, value=400, description='Y范围 (m)')):\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 创建网格\n",
|
|||
|
|
" x = np.linspace(-x_range, x_range, 100)\n",
|
|||
|
|
" y = np.linspace(-y_range, y_range, 100)\n",
|
|||
|
|
" X, Y = np.meshgrid(x, y)\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 计算重力异常\n",
|
|||
|
|
" gravity_data = gravity_anomaly_sphere(X, Y, R, rho, h)\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 绘制结果\n",
|
|||
|
|
" plot_gravity_anomaly(X, Y, gravity_data, R, rho, h)\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 显示参数信息\n",
|
|||
|
|
" print(f\"球体参数:\")\n",
|
|||
|
|
" print(f\"半径: {R} m\")\n",
|
|||
|
|
" print(f\"密度: {rho} kg/m³\")\n",
|
|||
|
|
" print(f\"埋深: {h} m\")\n",
|
|||
|
|
" print(f\"最大重力异常: {np.max(gravity_data):.2f} mGal\")\n",
|
|||
|
|
" print(f\"最小重力异常: {np.min(gravity_data):.2f} mGal\")\n",
|
|||
|
|
"\n",
|
|||
|
|
"def interactive_magnetic_anomaly():\n",
|
|||
|
|
" \"\"\"交互式磁异常计算\"\"\"\n",
|
|||
|
|
" \n",
|
|||
|
|
" @interact\n",
|
|||
|
|
" def magnetic_plot(R=FloatSlider(min=10, max=200, step=10, value=50, description='半径 (m)'),\n",
|
|||
|
|
" M=FloatSlider(min=100, max=5000, step=100, value=1000, description='磁化强度 (A/m)'),\n",
|
|||
|
|
" h=FloatSlider(min=20, max=200, step=10, value=100, description='埋深 (m)'),\n",
|
|||
|
|
" I=FloatSlider(min=-90, max=90, step=5, value=60, description='磁化倾角 (°)'),\n",
|
|||
|
|
" D=FloatSlider(min=0, max=360, step=10, value=0, description='磁化偏角 (°)'),\n",
|
|||
|
|
" x_range=FloatSlider(min=100, max=1000, step=50, value=400, description='X范围 (m)'),\n",
|
|||
|
|
" y_range=FloatSlider(min=100, max=1000, step=50, value=400, description='Y范围 (m)')):\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 创建网格\n",
|
|||
|
|
" x = np.linspace(-x_range, x_range, 100)\n",
|
|||
|
|
" y = np.linspace(-y_range, y_range, 100)\n",
|
|||
|
|
" X, Y = np.meshgrid(x, y)\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 计算磁异常\n",
|
|||
|
|
" magnetic_data = magnetic_anomaly_sphere(X, Y, R, M, h, I, D)\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 绘制结果\n",
|
|||
|
|
" plot_magnetic_anomaly(X, Y, magnetic_data, R, M, h, I, D)\n",
|
|||
|
|
" \n",
|
|||
|
|
" # 显示参数信息\n",
|
|||
|
|
" print(f\"球体参数:\")\n",
|
|||
|
|
" print(f\"半径: {R} m\")\n",
|
|||
|
|
" print(f\"磁化强度: {M} A/m\")\n",
|
|||
|
|
" print(f\"埋深: {h} m\")\n",
|
|||
|
|
" print(f\"磁化倾角: {I}°\")\n",
|
|||
|
|
" print(f\"磁化偏角: {D}°\")\n",
|
|||
|
|
" print(f\"最大磁异常: {np.max(magnetic_data):.2f} nT\")\n",
|
|||
|
|
" print(f\"最小磁异常: {np.min(magnetic_data):.2f} nT\")"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"## 5. 使用示例"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 5,
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"name": "stdout",
|
|||
|
|
"output_type": "stream",
|
|||
|
|
"text": [
|
|||
|
|
"=== 示例1: 固定参数重力异常 ===\n"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"data": {
|
|||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAABdIAAAJNCAYAAAA1YYe/AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzs3XdcU9f7B/BPFknYS5YDHLgAF6CCs6JY96p1171rtdVaaW3VDm31W7XW0aXgHq117z1xoeAWN1ZBFNmb5P7+yI9IZCsQxuf9euXlzcm55z73gJo8Ofe5IkEQBBARERERERERERERUY7E+g6AiIiIiIiIiIiIiKg0YyKdiIiIiIiIiIiIiCgPTKQTEREREREREREREeWBiXQiIiIiIiIiIiIiojwwkU5ERERERERERERElAcm0omIiIiIiIiIiIiI8sBEOhERERERERERERFRHphIJyIiIiIiIiIiIiLKAxPpRERERERERERERER5YCKdqIK4evUqRo4ciZo1a0KpVEKpVMLZ2Rljx47FpUuX9BbX7NmzIRKJdNqWL1+OgICAYjvmy5cvIZfLIRKJ9HruxSWnOS0KPj4+GDdunPZ5UFAQJk6cCDc3N5iYmMDW1hbt27fH0aNHs+07bNgwiEQi7cPIyAhOTk7o3r07/P39kZqamm2f1q1bY8qUKUV+HkREREREREREhcVEOlEF8Pvvv8Pd3R3nz5/H5MmTsXv3buzZswdTpkzBjRs34Onpifv37+sltlGjRiEwMFCnrbgT6WvXrkVaWhoAYOXKlcV2nPJkx44dOHPmDL7++mtt28aNG3HhwgWMGDECO3bswF9//QW5XA4fHx+sWbMm2xhKpRKBgYEIDAzE7t278e2338LIyAijR4+Gu7s7/vvvP53+3333HZYvX447d+4U+/kREREREREREeVFJAiCoO8giKj4nDlzBq1bt0aXLl3wzz//wMDAIFufv//+Gy1atICDg0Ou4yQlJcHQ0LA4Q9VydXWFtbU1jh8/Xizju7m5ITIyEo6OjggNDUV4eDiUSmWxHEsfZs+ejTlz5qAo/3lv1qwZatSogY0bN2rbIiMjYWNjo9NPpVKhSZMmSExMxL1797Ttw4YNwz///IOEhIRsYx88eBBdu3ZFkyZNcO7cOZ3X3Nzc4OXlhT/++KPIzoWIiIiIiIiIqLC4Ip2onJs7dy4kEgl+//33HJPoANC3b1+dJPqwYcNgbGyMa9euwdfXFyYmJvDx8QEAHDp0CD169ECVKlWgUChQq1YtjB07Fi9fvtTuv337dohEIhw5ciTbsVasWAGRSISrV68CyF6GxMnJCTdu3MCJEye0ZUCcnJyQkJAAc3NzjB07NtuYjx49gkQiwYIFC/Kdj/Pnz+P69esYMmQIRo8ejdjYWGzdujVbv7Zt28LV1RUXL15Eq1atYGhoiBo1auDHH3+EWq3W6RsWFobBgwfDxsYGcrkc9erVw88//6zT79GjRxCJRFiwYAF++uknODk5QalUom3btggNDUV6ejpmzJgBBwcHmJmZoVevXoiMjNQ5zubNm+Hr6wt7e3solUrUq1cPM2bMQGJiYp7nPHLkSFhaWiIpKSnba+3atYOLi0ue+1+5cgUXLlzAkCFDdNrfTKIDgEQigbu7O548eZLnmFn5+vpi9OjROH/+PE6ePKnz2pAhQ7BhwwbEx8cXeDwiIiIiIiIioqLGRDpROaZSqXDs2DF4eHjA3t6+UPumpaWhe/fuaNeuHXbs2IE5c+YAAO7fvw8vLy+sWLECBw8exDfffIPz58+jZcuWSE9PBwB07doVNjY28Pf3zzZuQEAAmjRpggYNGuR43G3btqFGjRpo3LixtgzItm3bYGxsjBEjRmD9+vWIjY3V2Wf58uUwMDDAiBEj8j2vzFIuI0aMQP/+/WFoaJhreZeIiAgMGjQIgwcPxs6dO9GpUyf4+flh3bp12j4vXryAt7c3Dh48iO+++w47d+5E+/btMW3aNHz88cfZxly2bBnOnDmDZcuW4a+//sLt27fRrVs3jBw5Ei9evMCqVaswf/58HD58GKNGjdLZ9+7du+jcuTNWrlyJ/fv3Y8qUKdiyZQu6deuW5zlPnjwZ0dHR2LBhg077zZs3cezYMUycODHP/Xfv3g2JRILWrVvn2Q8AMjIycOrUqXyT82/q3r07AGRLpLdt2xaJiYnFdnUCEREREREREVFBSPUdABEVn5cvXyI5ORmOjo7ZXlOpVDqlPyQSic7K8PT0dHzzzTcYPny4zn5ZbzYpCAK8vb3Rtm1bODo6Yt++fejevTukUikGDx6MFStWIDY2FmZmZgCAW7du4cKFC/j1119zjblx48ZQKpUwNTVF8+bNdV77+OOP8csvv8Df3197E8qUlBSsWrUKAwYMgJWVVZ7zkZSUhM2bN6N58+aoX78+AM1q/DVr1uD+/fuoWbOmTv+oqCjs3bsXTZs2BQC0b98ex48fx4YNG/DRRx8BABYuXIinT5/i/Pnz2n4dO3aESqXCb7/9hilTpqB27draMc3NzbF9+3aIxZrvMV++fIkpU6agbt262LFjh7bf7du3sXjxYsTFxcHU1BQAMHPmTJ25b9GiBerVq4c2bdrg6tWruX450aBBA7Rp0wbLli3TSc4vXboUpqam2nPJTWBgIJydnWFsbJxnP0BzhcG9e/ewffv2fPtmlfk7+uzZM532xo0bQyQS4cyZM/l+YUBEREREREREVFy4Ip2ognJ3d4dMJtM+fv7552x9+vTpk60tMjIS48aNQ9WqVSGVSiGTybRJ0Fu3bmn7jRgxAsnJydi8ebO2zd/fH3K5HAMHDnyrmGvUqIGuXbti+fLl2i8BNmzYgKioqBxXf79py5YtiIuL01m5PmLECAiCkOPqeTs7O21yPFODBg3w+PFj7fOjR4+ifv362foNGzYMgiDg6NGjOu2dO3fWJtEBoF69egCALl266PTLbA8LC9O2PXjwAAMHDoSdnR0kEglkMhnatGkDQHfuczJ58mQEBwfjzJkzAIC4uDisXbsWQ4cOzTdB/uzZsxzLuLzpr7/+wg8//ICpU6eiR48e+fbPKrd67jKZDObm5nj69GmhxiMiIiIiIiIiKkpMpBOVY9bW1lAqlTqJ30wbNmzAxYsXsXPnzhz3NTQ01K6EzqRWq+Hr64t///0X06dPx5EjR3DhwgXtDSKTk5O1fV1cXODp6alNUKtUKqxbtw49evSApaXlW5/T5MmTcffuXRw6dAiAplSKl5cXmjRpku++K1euhEKhwPvvv4+YmBjExMSgQYMGcHJyQkBAAFQqlU7/nFa4y+VynfOMiorKsWxOZs35qKgonfY3zz2zbn1u7SkpKQCAhIQEtGrVCufPn8f333+P48eP4+LFi/j3338B6M59Tnr06AEnJycsW7YMgKbETmJiYr5lXTLHVigUefbx9/fH2LFjMWbMmALVqn9T5u9oTje8VSgU+Z4fEREREREREVFxYmkXonJMIpGgXbt2OHjwIMLDw3USvpmlTR49epTjvlnLvGS6fv06QkJCEBAQgKFDh2rb7927l+MYw4cPx4QJE3Dr1i08ePAA4eHh2UrFFFa7du3g6uqKpUuXwtjYGJcvX9apWZ6b0NBQnD59GgBQrVq1HPscOHAAnTt3LlQ8VlZWCA8Pz9aeWaLE2tq6UOPl5ujRo3j27BmOHz+uXYUOADExMQXaXywWY+LEifjyyy/x888/Y/ny5fDx8UGdOnXy3dfa2hqvXr3K9XV/f3+MGjUKQ4cOxW+//Zbj705+Mr/Qadu2bbbXoqOji2weiYiIiIiIiIjeBlekE5Vzfn5+UKlUGDdunPZmoG8rM0Eql8t12n///fcc+w8YMAAKhQIBAQEICAhA5cqV4evrm+9x3lz1/aZPPvkEe/bsgZ+fH2xtbdG3b998x8y8oeiff/6JY8eO6Tz27t0LmUyGVatW5TvOm3x8fHDz5k1cvnxZp33NmjUQiUR47733Cj1mTgo79zkZNWoUDAwMMGjQINy5c6dA5XAAoG7dunjw4EGOrwUEBGDUqFEYPHgw/vrrr7dKoh86dAh
|
|||
|
|
"text/plain": [
|
|||
|
|
"<Figure size 1500x600 with 3 Axes>"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
"metadata": {},
|
|||
|
|
"output_type": "display_data"
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"name": "stdout",
|
|||
|
|
"output_type": "stream",
|
|||
|
|
"text": [
|
|||
|
|
"最大重力异常: 1.04 mGal\n",
|
|||
|
|
"最小重力异常: 0.01 mGal\n"
|
|||
|
|
]
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"# 示例1: 固定参数的重力异常计算\n",
|
|||
|
|
"print(\"=== 示例1: 固定参数重力异常 ===\")\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 参数设置\n",
|
|||
|
|
"R = 50 # 半径 (m)\n",
|
|||
|
|
"rho = 3000 # 密度 (kg/m³)\n",
|
|||
|
|
"h = 100 # 埋深 (m)\n",
|
|||
|
|
"x_range = 400 # X范围 (m)\n",
|
|||
|
|
"y_range = 400 # Y范围 (m)\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 创建网格\n",
|
|||
|
|
"x = np.linspace(-x_range, x_range, 100)\n",
|
|||
|
|
"y = np.linspace(-y_range, y_range, 100)\n",
|
|||
|
|
"X, Y = np.meshgrid(x, y)\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 计算重力异常\n",
|
|||
|
|
"gravity_data = gravity_anomaly_sphere(X, Y, R, rho, h)\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 绘制结果\n",
|
|||
|
|
"plot_gravity_anomaly(X, Y, gravity_data, R, rho, h)\n",
|
|||
|
|
"\n",
|
|||
|
|
"print(f\"最大重力异常: {np.max(gravity_data):.2f} mGal\")\n",
|
|||
|
|
"print(f\"最小重力异常: {np.min(gravity_data):.2f} mGal\")"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 6,
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"name": "stdout",
|
|||
|
|
"output_type": "stream",
|
|||
|
|
"text": [
|
|||
|
|
"=== 示例2: 固定参数磁异常 ===\n"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"data": {
|
|||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAABdEAAAJOCAYAAABYwk4SAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzs3Xd4VNXWx/HvZJJMem+EXkINKFIDShEIIEVERUW5IIrei8JFwYJYggooNnzBfhW8IKJXxAKCgFJEepVeBOkhlJBCeua8f8QMGUJIQkImCb/P88zjyTlrztlnJ8TMmjVrmwzDMBARERERERERERERkXycHD0AEREREREREREREZHySkl0EREREREREREREZECKIkuIiIiIiIiIiIiIlIAJdFFRERERERERERERAqgJLqIiIiIiIiIiIiISAGURBcRERERERERERERKYCS6CIiIiIiIiIiIiIiBVASXURERERERERERESkAEqii4iIiIiIiIiIiIgUwNnRAxARuZzdu3cTFRV1xZjffvuNW2655YoxGzZsICsrq9BzVYS4iIiIK8aIiIiIiIiIiEjpUxJdRMql7OxsIiMjWbVq1WWP33zzzZhMpkJjsrOzi3SuihAnIiIiIiIiIiJlT+1cREREREREREREREQKoCS6iIiIiIiIiIiIiEgBlEQXERERERERERERESmAkugiIiIiIiIiIiIiIgVQEl1EREREREREREREpABKoouIiIiIiIiIiIiIFEBJdBERERERERERERGRAiiJLiIiIiIiIiIiIiJSACXRRUREREREREREREQKoCS6iIiIiIiIiIiIiEgBlEQXERERERERERERESmAkugiIiIiIiIiIiIiIgVwdvQAREQux2w2s23bNvz8/C57PDs7m+zs7EJjnJycinSuihAnIiIiIiIiIiJlz2QYhuHoQYiIiIiIiIiIiIiIlEcqbRQRERERERERERERKYCS6HLVZsyYgclksj2cnZ2pUqUK9957L/v37y+163Tq1MnuOrmPHj165IvNzMxk/Pjx1KpVC4vFQsOGDZk6dWqpjeVq5I6/Tp06XO6DHytXrrTd04wZM0p8vWPHjjFq1Cg6duyIn59foeddunQpUVFReHh4EBQUxJAhQ4iLi8sXV5y5PXjwIP3798fPzw8vLy+6devG5s2bCxzDmTNnsFgsmEwmNm7ceMX7++OPPzCZTGzZsuWKcaVt+fLlmEwmvvnmm2t2jVWrVvHwww/TokUL23z89ddfBcZPnTqVhg0bYrFYqF27NuPHjyczMzNfXFxcHEOGDCEoKAgPDw+ioqL45Zdfrtl9FIXJZOLxxx8vtfPt2LGDu+++m+DgYCwWC7Vq1WL48OH54or6s3ns2DGio6Px8PCgffv27N69+6rHlvf3ltlsxt/fnxtuuIFHH32UtWvXXvV5L6c4/05feOEFAgMDqVq1Kp988kmpjkNERERERESkMlESXUps+vTprFmzhqVLl/L444/zww8/cPPNNxMfH19q16hTpw5r1qyxe0yZMiVf3PDhw5k0aRKPPfYYP//8M3fccQf//ve/mThxYqmN5Wp4e3tz6NAhfv3113zHPvvsM3x8fErtWgcOHOCLL77A1dWV22677YqxK1asoGfPnoSGhvL999/z7rvvsnTpUrp06UJ6erpdbFHn9vTp09xyyy3s27ePzz77jK+//pq0tDQ6derE3r17LzuOmTNnkpGRAcCnn356xTHPnTuX2rVr07x588KmosL55ZdfWLp0KTVq1KBdu3ZXjJ0wYQL//ve/6d+/Pz///DPDhw9n4sSJPPbYY3Zx6enpdOnShV9++YV3332X77//ntDQUHr06MGKFSuu5e2UmWXLltG6dWsSExP58MMPWbx4Ma+88gpubm52ccX52Rw2bBihoaEsXLiQtm3bcvfdd5dojHfddRdr1qxh1apVzJkzh3/84x+sXbuWqKgo/v3vf5fo3HkV9d/pV199xezZs5kxYwZvv/02zz33HBs2bCi1cYiIiIhUBLt378bPz++Kj/3791eauML89ttvWCwWDh8+zKFDh/D29ubOO++8bOzs2bMxmUx89NFHJf4+bN68ma5du+Ll5YWfnx/9+/fn4MGDJT5vXpmZmdStW/eyeQQRkSIxRK7S9OnTDcDYsGGD3f7x48cbgPHZZ5+VynU6duxoNGnSpNC4HTt2GCaTyZg4caLd/mHDhhnu7u7G2bNnS2U8xZU7/rZt2xoDBw60O5aYmGh4eHgYw4YNMwBj+vTpJb5edna2bXvDhg1XPG+rVq2Mxo0bG5mZmbZ9v//+uwEY77//vm1fceb2qaeeMlxcXIy//vrLti8hIcEICgoyBgwYcNlxREZGGiEhIUarVq0MX19fIyUlpcD7a9y4sTFmzJgCj18ry5YtMwDjf//73zW7Rt7v3RtvvGEAxqFDh/LFnTlzxnBzczMeeeQRu/0TJkwwTCaTsXPnTtu+9957zwCM1atX2/ZlZmYajRs3Nlq3bl36N1FEgPHYY4+V+DwXLlwwqlSpYvTq1cuwWq1XjC3Oz6anp6eRkJBg+9rf3984c+bMVY2xoHvNysoyhg4dmu/f29Uqzr/Txx57zHj33XdtXz/xxBPGm2++WeIxiIiIiFQk27dvN9q3b1/g8fbt2xu7d++uNHFXYrVajZtuusnu79aPP/7YAIwvvvjCLvbkyZNGQECA0b179yuesyh2795teHt7G7fccouxYMECY+7cuUaTJk2M8PBwIy4ursTnz2vGjBkl+rteRK5vqkSXUteyZUsATp06VabX/e677zAMgwcffNBu/4MPPkhqaiqLFi2y7evUqRORkZGsWbOGdu3a4e7uTq1atZg+fToACxYs4KabbsLDw4OmTZvaPfdqDR06lG+//Zbz58/b9s2ZMweAe++9t8Tnz+XkVLR/1sePH2fDhg0MGjQIZ2dn2/527dpRv3595s2bZ9tXnLmdN28et956KzVr1rTt8/HxoX///vz4449kZWXZnWPdunXs2LGDQYMGMWzYMBISEpg7d+5lx7xnzx527dplq4b466+/MJlMvPHGG7z++uvUqlULd3d3OnXqxL59+8jMzOTZZ58lPDwcX19f7rjjjsu2qimOzMxMxo0bR3h4OD4+PnTt2rXACvviKur3btGiRaSlpV32+2EYBt99951t37x582jQoAFRUVG2fc7OzjzwwAOsX7+e48eP2/bntliZPn06DRo0wN3dnZYtW7J27VoMw+CNN96gdu3aeHl5ceutt3LgwIGS3XAp+N///sfJkyd56qmnMJlMV4wtzs9m3bp1ef/990lMTGTWrFm4uLgQEBBQqmM3m81MmzaNoKAg3njjjRKfrzj/TuvWrcvXX3/N0aNH2bt3LwsWLKBBgwYlHoOIiIiIVEyLFi1i8+bNjBgxwrZv2LBh9OzZkxEjRnDy5Enb/kcffRTDMAr9FHFRvPjii1gsFubPn89tt91G//79WbBgAadPn+bNN98s8fnzuu+++0qtel5Erj9KokupO3ToEAD169e322+1WsnKyir0kZ2dne+cf/75JwEBATg7O1O3bl3GjRtHamqqXcyOHTsIDg4mLCzMbn+zZs1sx/OKjY3lwQcf5OGHH+b777+nadOmDB06lJdffpmxY8fy9NNPM3fuXLy8vOjXrx8nTpwo0bzce++9mM1mvvzyS9u+Tz/9lLvuuuuy7VwMwyjSfF2alC6q3PnInZ+8mjVrZjdfRZ3b1NRU/vzzzwLPmZqamu9jebl/eA0dOpR7770XDw+PAv8Ymzt3LlWrVqVNmzZ2+9977z1+//133nvvPf7zn/+wZ88e+vTpw0MPPcTp06f57LPPmDx5MkuXLuXhhx++4rwU5rnnnuPw4cP85z//4eOPP2b//v306dPH7ue2JD/rRZE7302bNrXbX6VKFYKCgvJ97wr6fgDs3LnTbv/8+fP5z3/+w2uvvcaXX35JUlISvXr1YvTo0fz+++9MmzaNjz/+2PZmhnGZPv9Xq6g/73mvuXLlSgCys7O5+eabcXV1xd/fn/vuu8/u32xxfzbff/9
|
|||
|
|
"text/plain": [
|
|||
|
|
"<Figure size 1500x600 with 3 Axes>"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
"metadata": {},
|
|||
|
|
"output_type": "display_data"
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"name": "stdout",
|
|||
|
|
"output_type": "stream",
|
|||
|
|
"text": [
|
|||
|
|
"最大磁异常: 95523.32 nT\n",
|
|||
|
|
"最小磁异常: -7447.79 nT\n"
|
|||
|
|
]
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"# 示例2: 固定参数的磁异常计算\n",
|
|||
|
|
"print(\"=== 示例2: 固定参数磁异常 ===\")\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 参数设置\n",
|
|||
|
|
"R = 50 # 半径 (m)\n",
|
|||
|
|
"M = 1000 # 磁化强度 (A/m)\n",
|
|||
|
|
"h = 100 # 埋深 (m)\n",
|
|||
|
|
"I = 60 # 磁化倾角 (°)\n",
|
|||
|
|
"D = 0 # 磁化偏角 (°)\n",
|
|||
|
|
"x_range = 400 # X范围 (m)\n",
|
|||
|
|
"y_range = 400 # Y范围 (m)\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 创建网格\n",
|
|||
|
|
"x = np.linspace(-x_range, x_range, 100)\n",
|
|||
|
|
"y = np.linspace(-y_range, y_range, 100)\n",
|
|||
|
|
"X, Y = np.meshgrid(x, y)\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 计算磁异常\n",
|
|||
|
|
"magnetic_data = magnetic_anomaly_sphere(X, Y, R, M, h, I, D)\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 绘制结果\n",
|
|||
|
|
"plot_magnetic_anomaly(X, Y, magnetic_data, R, M, h, I, D)\n",
|
|||
|
|
"\n",
|
|||
|
|
"print(f\"最大磁异常: {np.max(magnetic_data):.2f} nT\")\n",
|
|||
|
|
"print(f\"最小磁异常: {np.min(magnetic_data):.2f} nT\")"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 7,
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"name": "stdout",
|
|||
|
|
"output_type": "stream",
|
|||
|
|
"text": [
|
|||
|
|
"=== 示例3: 3D重力异常可视化 ===\n"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"data": {
|
|||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAvwAAAKQCAYAAADjbq6uAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzs/Xl4XGl5541/n3NO7ZsW27JkS/K+tHu13e6WTIclpEl3AhlgMj28E7YAoV8SGELCO8MSkitX3uEXwiSEENYQCNBJ+IVAZzrTgTRcTNMN3dBNW5JlW7ZkybIsybIslbbazva8f5TOcVWpJNVyTtUj6f5cl2isOnrOc06d5Xvfz70wzjkHQRAEQRAEQRCbEqneEyAIgiAIgiAIwj1I8BMEQRAEQRDEJoYEP0EQBEEQBEFsYkjwEwRBEARBEMQmhgQ/QRAEQRAEQWxiSPATBEEQBEEQxCaGBD9BEARBEARBbGJI8BMEQRAEQRDEJoYEP0EQBEEQBEFsYkjwEwRBEARBEMQmRqn3BAiC2LhcuHABXV1da27zzDPP4IEHHlhzmxdeeAG6rq87Vj23O3jwIHp6evCRj3wEZ8+exfT0NAKBAA4fPozf/u3fxm/8xm/k/c0rXvEKPP300wAAxhhCoRBaWlpwzz334JFHHsEb3vAGSNItn8tb3/pW/Mu//Muqc3jjG9+IL3/5y2vOkyAIgiCKQYKfIIiKMQwDt99+O5599tmin7/sZS8DY2zdbQzDKGmsem4HAHNzc2hvb8eb3vQm7Nq1C4lEAo899hje/OY348qVK/joRz+a97f79u3DY489BgBIJBIYGRnB448/jl//9V/HAw88gCeeeAKxWMz+/Fvf+hZe/epXr5jD97//ffzN3/xN0fkRBEEQxHqQ4CcIgiiRV7ziFXjFK16R97tf/dVfxcjICL74xS+uEPyBQAD3339/3u/e+c534itf+Qp+8zd/E7/1W7+Fb37zm25PmyAIgtjiUAw/QRBElWzbtg2KUrr/5O1vfzsefvhh/NM//RNGR0ddnBlBEARBkOAnCIIoG9M0oes6pqen8dnPfhbf+9738N/+238ra4zXve514JzjmWeecWmWBEEQBJGFQnoIgiDK5D3veQ++8IUvAAC8Xi8+/elP493vfndZY3R2dgIAJiYmHJ8fQRAEQeRCgp8gCKJMPvzhD+Od73wnbty4gSeeeAK/8zu/g0Qigd///d8veQzOuYszJAiCIIhbkOAnCIIok46ODnR0dAAAHn74YQDAhz70Ibz1rW/F9u3bSxrDit1va2tzZ5IEQRAEsQzF8BMEQVTJqVOnoOs6hoeHS/6b//W//hcYY/iFX/gFF2dGEARBECT4CYIgquaHP/whJEnCvn37Str+K1/5Cv7t3/4Nb3rTm+yVAoIgCIJwCwrpIQiCKJHf+q3fQjQaxalTp9DS0oKbN2/in/7pn/DNb34TH/zgB1eE86RSKTz//PP2/x8eHsbjjz+Of/3Xf8XLX/5yfP7zn6/HYRAEQRBbDBL8BEEQJdLV1YWvfOUr+Lu/+zvMzc0hHA7jrrvuwte//nX8xm/8xorth4eH0dXVBQAIhUJoaWnB8ePH8U//9E94wxveAEmiRVaCIAjCfUjwEwRBlMjb3/52vP3tby9p2//zf/6Pu5MhCIIgiBIh9xJBEARBEARBbGLIw08QRMXIsoze3l40NDQU/dwwDBiGse42kiSVNFY9t3ObQCCAN77xjWCMrfiMc47Xv/71rs+BIAiC2JwwTt1fCIIgCIIgCGLTQiE9BEEQBEEQBLGJIcFPEARBEARBEJsYiuEnCIIgCIIgak46nYaqqjXfr9frhd/vr/l+6wkJfoIgCIIgCKKmpNNpbI80YEnP1HzfO3fuxMjIyJYS/ST4CYIgCIIgiJqiqiqW9Az+nzt+ET65dnI0Y+j4xNkfQFVVEvwEQRAEQRAE4TZBRYFf9tRsf/LKysdbAkraJQiCIAiCIIhNDAl+giAIgiAIgtjEUEgPQRAEQRAEURckxiAV6TDu5v62IuThJwiCIAiCIIhNDHn4CYIgCIIgiLogsexPLfe3FSEPP0EQBEEQBEFsYkjwEwRBEARBEMQmhkJ6CIIgCIIgiLogL//Ucn9bEfLwEwRBEARBEMQmhjz8BEEQBEEQRF2QUOOk3drtSii26nETBEEQBEEQxJaAPPwEQRAEQRBEXZAYg0yNt1yHPPwEQRAEQRAEsYkhwU8QBEEQBEEQmxgK6SEIgiAIgiDqAnXarQ3k4ScIgiAIgiCITQx5+AmCIAiCIIi6IDFW00RaStolCIIgCIIgCGLTQYKfIAiCIAiCIDYxJPgJgiAIgiCIuiDV4adcfvSjH+G1r30t2trawBjD448/vu7fPP300zhx4gT8fj/27duHz3/+8xXs2TlI8BMEQRAEQRDEKiQSCdx11134zGc+U9L2IyMjePjhh/HAAw/gzJkz+PCHP4z3ve99+Od//meXZ7o6lLRLEARBEARB1AUJDBJqmLRbwb4eeughPPTQQyVv//nPfx4dHR341Kc+BQA4evQoXnzxRXzyk5/EG9/4xrL37wTk4ScIgiAIgiC2FAsLC3k/mUzGsbGfe+45PPjgg3m/e81rXoMXX3wRmqY5tp9yIMFPEARBEARB1AWZ1f4HANrb2xGLxeyfj3/8444d0/Xr19HS0pL3u5aWFui6jps3bzq2n3KgkB6CIAiCIAhiSzE2NoZoNGr/2+fzOTo+K6j3zzkv+vtaQYKfIAiCIAiC2FJEo9E8we8kO3fuxPXr1/N+d+PGDSiKgubmZlf2uR4k+AmCIAiCIIi6sBk77XZ1deGJJ57I+92///u/4+TJk/B4PK7vvxgUw08QBEEQBEEQq7C0tISenh709PQAyJbd7OnpwdWrVwEAH/rQh/CWt7zF3v7RRx/F6OgoPvCBD+DChQv427/9W3z5y1/G7//+79dj+gDIw08QBEEQBEHUCSknkbZW+yuXF198Ea985Svtf3/gAx8AALz1rW/FV7/6VUxOTtriHwD27t2LJ598Er/7u7+Lv/7rv0ZbWxs+/elP160kJwAwbmUREARBEARBEEQNWFhYQCwWwxe7fhVBpXZhLkldw28996+Yn593LYZfRCikhyAIgiAIgiA2MRTSQxAEQRAEQdQFiVUWZlPN/rYi5OEnCIIgCIIgiE0MefgJgiAIgiCIuiChtt7nrerp3qrHTRAEQRAEQRBbAvLwEwRBEARBEHVBZgxyDRtv1XJfIkEefoIgCIIgCILYxJDgJwiCIAiCIIhNDIX0EARBEARBEHWBynLWBvLwEwRBEARBEMQmhjz8BEEQBEEQRF2QWfanlvvbipCHnyAIgiAIgiA2MST4CYIgCIIgCGITQyE9BEEQBEEQRF2QwCChdnE2tdyXSJDgJ4gtCuccmqYhnU5DURQoigJZliFJEtgWbUxCEARBEJsREvwEsQUxTROqqsIwDOi6DsMwkMlkwBiDJEnweDyQZRmKooAxRgYAQRAE4QpUlrM2kOAniC0E5xyGYUDTNHDOIUmS7dXnnNuf67puC31L+FsrAGQAEARBEMTGggQ/QWwRrBAewzAAwBb5FpaQlyTJ3p5zDl3XoWnaqgaAtT1BEARBEGJCgp8gtgCWV980zbwY/VzBX0g5BoAVAkQGAEEQBFEOEmOQa7hqLG3RFWoS/ASxibEEuq7rAFBVQu56BoA1PgB4vV7bCCADgCAIgiDqCwl+gtikmKZpe/UBOB57v5oBcObMGezcuRMtLS1gjOWF/5ABQBAEQeRCSbu1gQQ/QWwyOOd5Yr9WZTYtA8AyAmRZtudhrQAUGgBWFSCCIAiCINyDBD9BbCKKJebWWlBb+7Pi+3PnZhkAqqraxoFlAORWASIIgiC2BnKNY/hruS+RIMFPEJsE0zQxMzMDv98Pj8dTV+FcLBm4FAPAWhnITQImA4AgCIIgqoMEP0FscHJr57/wwgs4efIkvF5v3eZTqkAv1QAozAEgA4AgCIIgyoMEP0FsYApDeBhja5barBWVzCHXALD+3uoInNsFmAwAgiCIzYO0/FPL/W1FSPA
|
|||
|
|
"text/plain": [
|
|||
|
|
"<Figure size 1200x800 with 2 Axes>"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
"metadata": {},
|
|||
|
|
"output_type": "display_data"
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"# 示例3: 3D可视化\n",
|
|||
|
|
"print(\"=== 示例3: 3D重力异常可视化 ===\")\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 使用较小的网格以提高性能\n",
|
|||
|
|
"x = np.linspace(-200, 200, 50)\n",
|
|||
|
|
"y = np.linspace(-200, 200, 50)\n",
|
|||
|
|
"X, Y = np.meshgrid(x, y)\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 计算重力异常\n",
|
|||
|
|
"gravity_data = gravity_anomaly_sphere(X, Y, R=50, rho=3000, h=100)\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 3D可视化\n",
|
|||
|
|
"plot_3d_anomaly(X, Y, gravity_data, '球体重力异常3D图', '重力异常 (mGal)')"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"## 6. 交互式计算\n",
|
|||
|
|
"\n",
|
|||
|
|
"运行下面的代码块，使用滑块调整参数来实时查看重磁异常的变化："
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 8,
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"name": "stdout",
|
|||
|
|
"output_type": "stream",
|
|||
|
|
"text": [
|
|||
|
|
"=== 交互式重力异常计算 ===\n"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"data": {
|
|||
|
|
"application/vnd.jupyter.widget-view+json": {
|
|||
|
|
"model_id": "0f183af415864da1af4c7d3abaf4bfc1",
|
|||
|
|
"version_major": 2,
|
|||
|
|
"version_minor": 0
|
|||
|
|
},
|
|||
|
|
"text/plain": [
|
|||
|
|
"interactive(children=(FloatSlider(value=50.0, description='半径 (m)', max=200.0, min=10.0, step=10.0), FloatSlid…"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
"metadata": {},
|
|||
|
|
"output_type": "display_data"
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"# 交互式重力异常计算\n",
|
|||
|
|
"print(\"=== 交互式重力异常计算 ===\")\n",
|
|||
|
|
"interactive_gravity_anomaly()"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 9,
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"name": "stdout",
|
|||
|
|
"output_type": "stream",
|
|||
|
|
"text": [
|
|||
|
|
"=== 交互式磁异常计算 ===\n"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"data": {
|
|||
|
|
"application/vnd.jupyter.widget-view+json": {
|
|||
|
|
"model_id": "ff8fcbe1939e4af5899f75843faa2090",
|
|||
|
|
"version_major": 2,
|
|||
|
|
"version_minor": 0
|
|||
|
|
},
|
|||
|
|
"text/plain": [
|
|||
|
|
"interactive(children=(FloatSlider(value=50.0, description='半径 (m)', max=200.0, min=10.0, step=10.0), FloatSlid…"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
"metadata": {},
|
|||
|
|
"output_type": "display_data"
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"# 交互式磁异常计算\n",
|
|||
|
|
"print(\"=== 交互式磁异常计算 ===\")\n",
|
|||
|
|
"interactive_magnetic_anomaly()"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"## 7. 参数敏感性分析"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 10,
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"name": "stdout",
|
|||
|
|
"output_type": "stream",
|
|||
|
|
"text": [
|
|||
|
|
"=== 重力异常参数敏感性分析 ===\n"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"data": {
|
|||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAABdEAAASmCAYAAADBBeLHAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzs3XmczWX/x/HXWWbfDWMMxtIiURHJ0EKWIbtEy62EihZu0kIbQj8lKd2lLCntGvtSQ6KUJUKiyBJiZJt9P8vvjzEn01hmmDnfOTPv5+Mxj7v5nu/5Xu9zzbnHNZ9zfa/L5HQ6nYiIiIiIiIiIiIiISCFmowOIiIiIiIiIiIiIiJRVKqKLiIiIiIiIiIiIiJyDiugiIiIiIiIiIiIiIuegIrqIiIiIiIiIiIiIyDmoiC4iIiIiIiIiIiIicg4qoouIiIiIiIiIiIiInIOK6CIiIiIiIiIiIiIi56AiuoiIiIiIiIiIiIjIOaiILiIiIiIiIiIiIiJyDlajA4iIlFW//fYbMTEx5z3n+++/5+abbz7vOT/99BM2m+2C19J57jvviiuuOO85IiIiIuK5NI4vv+dpHC8iRlERXUTkHOx2Ow0bNmTt2rVnffymm27CZDJd8By73V6ka+k8950nIiIiIuWXxvHl9zwREaNoORcRERERERERERERkXNQEV1ERERERERERERE5BxURBcREREREREREREROQcV0UVEREREREREREREzkFFdBERERERERERERGRc1ARXURERERERERERETkHFREFxERERERERERERE5BxXRRURERERERERERETOQUV0EREREREREREREZFzUBFdREREREREREREROQcVEQXERERERERERERETkHFdFFRERERERERERERM7BanQAEZGyymKxsG3bNkJDQ8/6uN1ux263X/Acs9lcpGvpPPedJyIiIiLll8bx5fc8ERGjmJxOp9PoECIiIiIiIiIiIiIiZZE+xhMREREREREREREROQcV0UVEREREREREREREzkFFdBERERERERERERGRc/DojUUdDgdHjhwhKCgIk8lkdBwRERERqaCcTiepqalERUVp47Mi0DheRERERMqCoo7jPbqIfuTIEWrWrGl0DBERERERAA4dOkSNGjWMjlHmaRwvIiIiImXJhcbxHl1EDwoKAvJeZHBwsFvbdjgcHD9+nCpVqmi2URGpz4pPfVZ86rPiU58Vn/qs+NRnF0f9VnxG9VlKSgo1a9Z0jU/l/DSO9yzqs+JTnxWf+qz41GfFpz4rPvVZ8anPis/IPivqON6ji+j5t34GBwcbMvjOysoiODhY/4coIvVZ8anPik99Vnzqs+JTnxWf+uziqN+Kz+g+09IkRaNxvGdRnxWf+qz41GfFpz4rPvVZ8anPik99Vnxloc8uNI7XT1JERERERERERERE5BxURBcREREREREREREROQcV0UVEREREREREREREzsGj10QXERER8WR2u53c3NxzPu5wOMjNzSUrK0vrKRZRafaZt7e3fg4iIiIiBrjQuLks0Ri++Eqzz7y8vLBYLJd8HRXRRURERNzM6XRy9OhRkpKSLniew+EgNTVVG1YWUWn2mdlspk6dOnh7e5fodUVERETk7Io6bi5LNIYvvtLus9DQUCIjIy/p2iqii4iIiLhZ/h8CERER+Pv7n3Mw53Q6sdlsWK1WDcCLqLT6zOFwcOTIERISEoiOjtbPQ0RERMQNijpuLks0hi++0uozp9NJRkYGx44dA6BatWoXfS0V0UVERETcyG63u/4QCA8PP++5GoAXX2n2WZUqVThy5Ag2mw0vL68SvbaIiIiIFFSccXNZojF88ZVmn/n5+QFw7NgxIiIiLnppFy3MIyIiIuJG+Ws5+vv7G5xEiit/GRe73W5wEhEREZHyT+NmKSn576FLWVdfRXQRERERA2hWiufRz0xERETE/TQGk0tVEu8hFdFFRERERERERERERM5BRXQRERERERERERERuaB+/frRvXt3o2O4nYroIiIiIlIk/fr1w2QyYTKZsFqtREdHM3jwYBITEy/qeidPnqRDhw5ERUXh4+NDzZo1eeyxx0hJSSlw3vbt27n11lvx8/OjevXqjB07FqfTWRIvSURERESkTPvzzz8ZMGAAderUwc/Pj8suu4wXX3yRnJycAucdPHiQLl26EBAQQOXKlRkyZEihc4oyrl6zZg1NmjTB19eXunXrMm3atFJ7baNHj+auu+666OePGjWK2NhY2rdvz3PPPVeCyQpTEV1EREREiqxDhw4kJCTw559/MmPGDBYvXswjjzxyUdcym81069aNRYsWsXv3bmbPns3KlSsZNGiQ65yUlBTatWtHVFQUP/30E1OnTmXSpElMnjy5pF6SiIiIiEiZ9fvvv+NwOHj33XfZsWMHr7/+OtOmTWPUqFGuc+x2O506dSI9PZ21a9fy2WefERcXxxNPPOE6pyjj6v3793P77bdz8803s2XLFkaNGsWQIUOIi4srlde2aNEiunXrdtHPnzBhAhaLhR07djB79uySC3YWKqKLiIiISJH5+PgQGRlJjRo1aN++PX369CE+Pv6irhUWFsbgwYNp2rQptWrVok2bNjzyyCN8//33rnM+/vhjsrKymD17Ng0bNqRnz56MGjWKyZMnu2bNjB49mkaNGjFr1ixq1arluq7dbueVV14hMjKSiIgIxo8fXyJ9ICIiIiJyIdu2baNZs2b4+vpSpUoVRo8efVHX6dChA++//z7t27enbt26dO3alREjRjBv3jzXOfHx8ezcuZOPPvqIxo0b07ZtW1577TWmT5/uusuzKOPqadOmER0dzZQpU6hfvz4DBw6kf//+TJo06Zz5Nm/eXGisPW7cOCIiIggKCmLgwIE888wzNGrUqMDzDh06xK+//krHjh0B8Pb25t1336Vz5874+/tTv3591q1bx549e2jVqhUBAQHExMSwd+/eAtdZtmwZGzduZMiQIRfVv0WlIrqIiIiIgZxOJxk5NkO+LnVJlH379vHVV1/h5eXlOtaxY0cCAwPP+3UuR44cYd68edx6662uY+vWrePWW2/Fx8fHdSw2NpYjR47w559/uo7t3buX5cuXs3z5cubMmcP7779Pp06d+Ouvv1izZg0TJ07kueeeY/369Zf0mkVERETEOJ4ydj58+DC33XYbLVq0YNu2bUyZMoVXX32VWbNmATBo0KALjpkPHjx4zusnJydTqVIl1/fr1q2jYcOGREVFuY7FxsaSnZ3N5s2bXedcaFy9bt062rdvX6Ct2NhYNm3aRG5ubqEcq1evpk2bNowZM4Znn30WyCvWjx8/nokTJ7J582aio6N55513Cj130aJF3HLLLYSGhrqOjRs3jvvuu4+tW7dy1VVXcc899/Dwww8zcuRINm3aBMBjjz0G5L0XsrOzAfDz82Ps2LEXvcxkUVhL7coiIiIickGZuXaufuFrQ9reOTYWf+/iDQeXLFlCYGAgdrudrKwsgAK3gM6YMYPMzMxiXfPuu+9m4cKFZGZm0qVLF2bMmOF67OjRo9SuXbvA+VWrVnU9VqdOHQAcDgezZs0iMDCQK6+8ktatW7Nr1y6WLVuG2WymXr16TJw4kdWrV9O8efNi5fNk3333Ha+++iqbN28mISGB+fPnn3cjqH79+vHBBx8UOn711VezY8cOAGbPns0DDzxQ6JzMzEx8fX1LLLuIiIjIv3nK2HnWrFlERkby+uuvYzKZqFevHtu3b2fy5Mn079+fsWPHMmLEiPNe48yC+Jn27t3L1KlTee2111zHjh496hoj5wsLC8Pb25ujR4+6zrnQuPps16latSo2m40TJ05QrVo11/GFCxfSt29f3n33Xe6++27X8alTpzJgwADXePGFF14gPj6etLS0AtdduHBhoaVc+vXrR+/evQF4+umniYmJ4fnnnyc2NhaAoUOHuq5rs9m47bbbsFqtHD9+nKFDhxIWFnbWPisJKqKLiIiISJG1bt2ad955h4yMDGbMmMHu3bt5/PHHXY9Xr1692Nd8/fXXefHFF9m1axejRo1i+PDhvP3
|
|||
|
|
"text/plain": [
|
|||
|
|
"<Figure size 1500x1200 with 4 Axes>"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
"metadata": {},
|
|||
|
|
"output_type": "display_data"
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"# 重力异常参数敏感性分析\n",
|
|||
|
|
"print(\"=== 重力异常参数敏感性分析 ===\")\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 基准参数\n",
|
|||
|
|
"base_R, base_rho, base_h = 50, 3000, 100\n",
|
|||
|
|
"x_profile = np.linspace(-200, 200, 100)\n",
|
|||
|
|
"y_profile = 0\n",
|
|||
|
|
"\n",
|
|||
|
|
"fig, axes = plt.subplots(2, 2, figsize=(15, 12))\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 半径变化的影响\n",
|
|||
|
|
"radii = [30, 50, 70, 100]\n",
|
|||
|
|
"for R in radii:\n",
|
|||
|
|
" gravity = gravity_anomaly_sphere(x_profile, y_profile, R, base_rho, base_h)\n",
|
|||
|
|
" axes[0,0].plot(x_profile, gravity, label=f'R={R}m')\n",
|
|||
|
|
"axes[0,0].set_xlabel('X (m)')\n",
|
|||
|
|
"axes[0,0].set_ylabel('重力异常 (mGal)')\n",
|
|||
|
|
"axes[0,0].set_title('半径变化对重力异常的影响')\n",
|
|||
|
|
"axes[0,0].legend()\n",
|
|||
|
|
"axes[0,0].grid(True, alpha=0.3)\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 密度变化的影响\n",
|
|||
|
|
"densities = [2000, 3000, 4000, 5000]\n",
|
|||
|
|
"for rho in densities:\n",
|
|||
|
|
" gravity = gravity_anomaly_sphere(x_profile, y_profile, base_R, rho, base_h)\n",
|
|||
|
|
" axes[0,1].plot(x_profile, gravity, label=f'ρ={rho}kg/m³')\n",
|
|||
|
|
"axes[0,1].set_xlabel('X (m)')\n",
|
|||
|
|
"axes[0,1].set_ylabel('重力异常 (mGal)')\n",
|
|||
|
|
"axes[0,1].set_title('密度变化对重力异常的影响')\n",
|
|||
|
|
"axes[0,1].legend()\n",
|
|||
|
|
"axes[0,1].grid(True, alpha=0.3)\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 埋深变化的影响\n",
|
|||
|
|
"depths = [50, 100, 150, 200]\n",
|
|||
|
|
"for h in depths:\n",
|
|||
|
|
" gravity = gravity_anomaly_sphere(x_profile, y_profile, base_R, base_rho, h)\n",
|
|||
|
|
" axes[1,0].plot(x_profile, gravity, label=f'h={h}m')\n",
|
|||
|
|
"axes[1,0].set_xlabel('X (m)')\n",
|
|||
|
|
"axes[1,0].set_ylabel('重力异常 (mGal)')\n",
|
|||
|
|
"axes[1,0].set_title('埋深变化对重力异常的影响')\n",
|
|||
|
|
"axes[1,0].legend()\n",
|
|||
|
|
"axes[1,0].grid(True, alpha=0.3)\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 综合对比\n",
|
|||
|
|
"gravity_base = gravity_anomaly_sphere(x_profile, y_profile, base_R, base_rho, base_h)\n",
|
|||
|
|
"axes[1,1].plot(x_profile, gravity_base, 'k-', linewidth=2, label='基准参数')\n",
|
|||
|
|
"axes[1,1].set_xlabel('X (m)')\n",
|
|||
|
|
"axes[1,1].set_ylabel('重力异常 (mGal)')\n",
|
|||
|
|
"axes[1,1].set_title('基准参数重力异常剖面')\n",
|
|||
|
|
"axes[1,1].legend()\n",
|
|||
|
|
"axes[1,1].grid(True, alpha=0.3)\n",
|
|||
|
|
"\n",
|
|||
|
|
"plt.tight_layout()\n",
|
|||
|
|
"plt.show()"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "code",
|
|||
|
|
"execution_count": 11,
|
|||
|
|
"metadata": {},
|
|||
|
|
"outputs": [
|
|||
|
|
{
|
|||
|
|
"name": "stdout",
|
|||
|
|
"output_type": "stream",
|
|||
|
|
"text": [
|
|||
|
|
"=== 磁异常参数敏感性分析 ===\n"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"data": {
|
|||
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAABUoAAAJNCAYAAADnHU+iAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQABAABJREFUeJzs3XdYU9f/B/B3CCEJG9k4wIGKgnvhrrPu1ap1rzpbtWpbUevostVfq62rtVVpnbWODkfduAeiqCjiFgeIIrJXkvP7g28ikbAUCOP9ep48vTk599zPuV4ofDhDIoQQICIiIiIiIiIiIirDTIwdABEREREREREREZGxMVFKREREREREREREZR4TpURERERERERERFTmMVFKREREREREREREZR4TpURERERERERERFTmMVFKREREREREREREZR4TpURERERERERERFTmMVFKREREREREREREZR4TpURERERERERERFTmMVFKREREREREREREZR4TpURERERERERERFTmMVFKREREREREREREZR4TpURERERERERERFTmMVFKREREREREREREZR4TpURERERERERERFTmMVFKREREREREREREZR4TpURERERERERERFTmMVFKREREREREREREZR4TpVTmXb58GaNHj0bVqlWhVCqhVCrh6emJcePG4fz580aLa/78+ZBIJHplK1euhL+/f6Fd89mzZ5DL5ZBIJEbte2ExdE8LQvv27TF+/Hjd+6CgIEyaNAk+Pj6wsrKCs7MzOnTogMOHD2c5d8SIEZBIJLqXhYUFPDw80LNnT6xbtw6pqalZzmndujWmTp1a4P0gIiIiIiIiKsuYKKUy7eeff0bDhg1x9uxZTJkyBbt27cLu3bsxdepUXL16FY0bN8bt27eNEtuYMWNw+vRpvbLCTpSuX78eaWlpAIA1a9YU2nVKk7///hsnT57EZ599pivbvHkzzp07h1GjRuHvv//Gr7/+Crlcjvbt2+P333/P0oZSqcTp06dx+vRp7Nq1C59//jksLCzw/vvvo2HDhnj48KFe/S+++AIrV65EWFhYofePiIiIiIiIqKyQCCGEsYMgMoaTJ0+idevW6NatG7Zt2wYzM7Msdf7880+0aNECbm5u2baTlJQEc3PzwgxVx9vbGw4ODggICCiU9n18fBAVFQV3d3fcuHEDERERUCqVhXItY5g/fz4WLFiAgvy217RpU1SpUgWbN2/WlUVFRcHJyUmvnlqtRoMGDZCYmIhbt27pykeMGIFt27YhISEhS9v79+9H9+7d0aBBA5w5c0bvMx8fH/j6+mL16tUF1hciIiIiIiKisowjSqnM+vrrryGVSvHzzz8bTJICwLvvvquXJB0xYgQsLS1x5coVdOrUCVZWVmjfvj0A4MCBA+jVqxcqVKgAhUKBatWqYdy4cXj27Jnu/L/++gsSiQSHDh3Kcq1Vq1ZBIpHg8uXLALJOE/fw8MDVq1dx9OhR3TRtDw8PJCQkwNbWFuPGjcvS5r179yCVSrF48eJc78fZs2cREhKCoUOH4v3330dsbCy2b9+epV7btm3h7e2NwMBAtGrVCubm5qhSpQq++eYbaDQavbrh4eEYMmQInJycIJfL4eXlhe+++06v3r179yCRSLB48WJ8++238PDwgFKpRNu2bXHjxg2kp6dj5syZcHNzg42NDfr06YOoqCi96/zxxx/o1KkTXF1doVQq4eXlhZkzZyIxMTHHPo8ePRrlypVDUlJSls/atWuH2rVr53j+xYsXce7cOQwdOlSv/NUkKQBIpVI0bNgQDx48yLHNzDp16oT3338fZ8+exbFjx/Q+Gzp0KDZt2oT4+Pg8t0dERERERERE2WOilMoktVqNI0eOoFGjRnB1dc3XuWlpaejZsyfatWuHv//+GwsWLAAA3L59G76+vli1ahX279+PuXPn4uzZs2jZsiXS09MBAN27d4eTkxPWrVuXpV1/f380aNAAderUMXjdnTt3okqVKqhfv75umvbOnTthaWmJUaNGYePGjYiNjdU7Z+XKlTAzM8OoUaNy7Zd2qv2oUaMwcOBAmJubZzv9PjIyEoMHD8aQIUPwzz//oEuXLvDz88OGDRt0dZ4+fYrmzZtj//79+OKLL/DPP/+gQ4cOmDFjBj744IMsba5YsQInT57EihUr8Ouvv+L69evo0aMHRo8ejadPn2Lt2rVYtGgRDh48iDFjxuide/PmTXTt2hVr1qzBf//9h6lTp2Lr1q3o0aNHjn2eMmUKYmJisGnTJr3ya9eu4ciRI5g0aVKO5+/atQtSqRStW7fOsR4AqFQqHD9+PNfk66t69uwJAFkSpW3btkViYmKhjS4mIiIiIiIiKnMEURkUGRkpAIiBAwdm+UylUon09HTdS6PR6D4bPny4ACDWrl2bY/sajUakp6eL+/fvCwDi77//1n02bdo0oVQqxYsXL3Rl165dEwDEsmXLdGXz5s0Tr36J1q5dW7Rp0ybL9W7fvi1MTEzEkiVLdGXJycnC3t5ejBw5MsdYhRAiMTFRWFtbi2bNmun1VSKRiFu3bunVbdOmjQAgzp49q1deq1Yt0blzZ937mTNnGqw3YcIEIZFIRFhYmBBCiLt37woAom7dukKtVuvqLV26VAAQPXv21Dt/6tSpAoCIjY012BftvT969KgAIC5duqT7zNA9bdOmjahXr16WGK2trUV8fLzBa2h16dJF1KxZM8c6WrNnzxYAxF9//aVXPnz4cGFhYZHteaGhoQKAmDBhgl55WlqakEgk4tNPP83T9YmIiIiM7dKlS2LUqFGiSpUqQqFQCIVCIapVqybGjh0rAgMDjRaXoZ8RV6xYIdatW1do13z69KkwMzMTAIza98Ji6J4WhHbt2olx48bplc2ePVt069ZNuLm5CQBi+PDhBs/V/i6nfZmbmwt3d3fRo0cPsXbtWpGSkpLlnFatWokpU6bkK8Zbt24JMzMzcerUKRETEyPKly8vmjRpIlQqVZa6x48fFyYmJmLmzJn5uoYht2/fFn369BE2NjbCwsJCdOjQQQQFBb1xu696nXtCVJJwRCnRKxo2bAiZTKZ7fffdd1nq9OvXL0tZVFQUxo8fj4oVK8LU1BQymQzu7u4AgNDQUF29UaNGITk5GX/88YeubN26dZDL5Rg0aNBrxVylShV0794dK1eu1K2/uWnTJkRHRxscvfmqrVu3Ii4uTm/k6ahRoyCEMDj61cXFBU2aNNErq1OnDu7fv697f/jwYdSqVStLvREjRkAIkWUH+K5du8LE5OW3JC8vLwBAt27d9Oppy8PDw3Vld+7cwaBBg+Di4gKpVAqZTIY2bdoA0L/3hkyZMgXBwcE4efIkACAuLg7r16/H8OHDYWlpmeO5jx8/NjjN/lW//vorvvrqK0yfPh29evXKtX5mIpv1VGUyGWxtbfHo0aN8tUdERERkDNxEVR83Uc0/Q5uoAsCSJUsQHR2Nnj17ZrukmlZRbKI6Y8YMdOzYEb6+vrC1tcXatWtx7tw5fPvtt3r1kpKSMHLkSNSuXVs3S/F1PX36FK1atcKNGzewdu1abN26FSkpKWjbtm2BbwDLjWWptGOilMokBwcHKJVKvcSe1qZNmxAYGIh//vnH4Lnm5uawtrbWK9NoNOjUqRN27NiBTz75BIcOHcK5c+d0G/AkJyfr6tauXRuNGzfWJSDVajU2bNiAXr16oVy5cq/dpylTpuDmzZs4cOAAgIyp7L6+vmjQoEGu565ZswYKhQJvv/02Xrx4gRcvXqBOnTrw8PCAv78/1Gq1Xn17e/ssbcjlcr1+RkdHG1zWQLvma3R0tF75q33X/pCTXXlKSgoAICEhAa1atcLZs2fx5ZdfIiAgAIGBgdixYwcA/XtvSK9eveDh4YEVK1YAyFgCITExMddp99q2FQpFjnXWrVuHcePGYezYsXlaK/ZV2mfU0IZiCoUi1/4RERERGdvJkycxceJEdOnSBRcuXMDkyZPRvn17tGvXDpMmTcKJEyewdevWXDcRNbSufEGoUKECmjVrVihtZ2ft2rVwcnJC48aNsXnzZv5Mlwdff/01+vT
|
|||
|
|
"text/plain": [
|
|||
|
|
"<Figure size 1500x600 with 3 Axes>"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
"metadata": {},
|
|||
|
|
"output_type": "display_data"
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"source": [
|
|||
|
|
"# 磁异常参数敏感性分析\n",
|
|||
|
|
"print(\"=== 磁异常参数敏感性分析 ===\")\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 基准参数\n",
|
|||
|
|
"base_R, base_M, base_h, base_I, base_D = 50, 1000, 100, 60, 0\n",
|
|||
|
|
"\n",
|
|||
|
|
"fig, axes = plt.subplots(2, 3, figsize=(18, 12))\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 磁化强度变化的影响\n",
|
|||
|
|
"magnetizations = [500, 1000, 2000, 3000]\n",
|
|||
|
|
"for M in magnetizations:\n",
|
|||
|
|
" magnetic = magnetic_anomaly_sphere(x_profile, y_profile, base_R, M, base_h, base_I, base_D)\n",
|
|||
|
|
" axes[0,0].plot(x_profile, magnetic, label=f'M={M}A/m')\n",
|
|||
|
|
"axes[0,0].set_xlabel('X (m)')\n",
|
|||
|
|
"axes[0,0].set_ylabel('磁异常 (nT)')\n",
|
|||
|
|
"axes[0,0].set_title('磁化强度对磁异常的影响')\n",
|
|||
|
|
"axes[0,0].legend()\n",
|
|||
|
|
"axes[0,0].grid(True, alpha=0.3)\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 磁化倾角变化的影响\n",
|
|||
|
|
"inclinations = [0, 30, 60, 90]\n",
|
|||
|
|
"for I in inclinations:\n",
|
|||
|
|
" magnetic = magnetic_anomaly_sphere(x_profile, y_profile, base_R, base_M, base_h, I, base_D)\n",
|
|||
|
|
" axes[0,1].plot(x_profile, magnetic, label=f'I={I}°')\n",
|
|||
|
|
"axes[0,1].set_xlabel('X (m)')\n",
|
|||
|
|
"axes[0,1].set_ylabel('磁异常 (nT)')\n",
|
|||
|
|
"axes[0,1].set_title('磁化倾角对磁异常的影响')\n",
|
|||
|
|
"axes[0,1].legend()\n",
|
|||
|
|
"axes[0,1].grid(True, alpha=0.3)\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 磁化偏角变化的影响\n",
|
|||
|
|
"declinations = [0, 45, 90, 135]\n",
|
|||
|
|
"for D in declinations:\n",
|
|||
|
|
" magnetic = magnetic_anomaly_sphere(x_profile, y_profile, base_R, base_M, base_h, base_I, D)\n",
|
|||
|
|
" axes[0,2].plot(x_profile, magnetic, label=f'D={D}°')\n",
|
|||
|
|
"axes[0,2].set_xlabel('X (m)')\n",
|
|||
|
|
"axes[0,2].set_ylabel('磁异常 (nT)')\n",
|
|||
|
|
"axes[0,2].set_title('磁化偏角对磁异常的影响')\n",
|
|||
|
|
"axes[0,2].legend()\n",
|
|||
|
|
"axes[0,2].grid(True, alpha=0.3)\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 半径变化的影响\n",
|
|||
|
|
"radii = [30, 50, 70, 100]\n",
|
|||
|
|
"for R in radii:\n",
|
|||
|
|
" magnetic = magnetic_anomaly_sphere(x_profile, y_profile, R, base_M, base_h, base_I, base_D)\n",
|
|||
|
|
" axes[1,0].plot(x_profile, magnetic, label=f'R={R}m')\n",
|
|||
|
|
"axes[1,0].set_xlabel('X (m)')\n",
|
|||
|
|
"axes[1,0].set_ylabel('磁异常 (nT)')\n",
|
|||
|
|
"axes[1,0].set_title('半径对磁异常的影响')\n",
|
|||
|
|
"axes[1,0].legend()\n",
|
|||
|
|
"axes[1,0].grid(True, alpha=0.3)\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 埋深变化的影响\n",
|
|||
|
|
"depths = [50, 100, 150, 200]\n",
|
|||
|
|
"for h in depths:\n",
|
|||
|
|
" magnetic = magnetic_anomaly_sphere(x_profile, y_profile, base_R, base_M, h, base_I, base_D)\n",
|
|||
|
|
" axes[1,1].plot(x_profile, magnetic, label=f'h={h}m')\n",
|
|||
|
|
"axes[1,1].set_xlabel('X (m)')\n",
|
|||
|
|
"axes[1,1].set_ylabel('磁异常 (nT)')\n",
|
|||
|
|
"axes[1,1].set_title('埋深对磁异常的影响')\n",
|
|||
|
|
"axes[1,1].legend()\n",
|
|||
|
|
"axes[1,1].grid(True, alpha=0.3)\n",
|
|||
|
|
"\n",
|
|||
|
|
"# 综合对比\n",
|
|||
|
|
"magnetic_base = magnetic_anomaly_sphere(x_profile, y_profile, base_R, base_M, base_h, base_I, base_D)\n",
|
|||
|
|
"axes[1,2].plot(x_profile, magnetic_base, 'k-', linewidth=2, label='基准参数')\n",
|
|||
|
|
"axes[1,2].set_xlabel('X (m)')\n",
|
|||
|
|
"axes[1,2].set_ylabel('磁异常 (nT)')\n",
|
|||
|
|
"axes[1,2].set_title('基准参数磁异常剖面')\n",
|
|||
|
|
"axes[1,2].legend()\n",
|
|||
|
|
"axes[1,2].grid(True, alpha=0.3)\n",
|
|||
|
|
"\n",
|
|||
|
|
"plt.tight_layout()\n",
|
|||
|
|
"plt.show()"
|
|||
|
|
]
|
|||
|
|
},
|
|||
|
|
{
|
|||
|
|
"cell_type": "markdown",
|
|||
|
|
"metadata": {},
|
|||
|
|
"source": [
|
|||
|
|
"## 8. 总结\n",
|
|||
|
|
"\n",
|
|||
|
|
"本notebook实现了球体的重磁异常正演计算，主要特点包括：\n",
|
|||
|
|
"\n",
|
|||
|
|
"1. **理论基础完整**：包含了重磁异常正演的理论公式和参数说明\n",
|
|||
|
|
"2. **计算功能完善**：实现了球体重力异常和磁异常的数值计算\n",
|
|||
|
|
"3. **可视化丰富**：提供了等值线图、剖面图和3D图等多种可视化方式\n",
|
|||
|
|
"4. **交互性强**：支持通过滑块实时调整参数并查看结果变化\n",
|
|||
|
|
"5. **分析深入**：包含了参数敏感性分析，帮助理解各参数对异常的影响\n",
|
|||
|
|
"\n",
|
|||
|
|
"### 使用建议：\n",
|
|||
|
|
"- 对于重力异常，主要关注球体的密度、半径和埋深参数\n",
|
|||
|
|
"- 对于磁异常，除了几何参数外，磁化强度和方向参数也很重要\n",
|
|||
|
|
"- 可以通过交互式界面快速探索不同参数组合的效果\n",
|
|||
|
|
"- 参数敏感性分析有助于理解各参数的相对重要性\n",
|
|||
|
|
"\n",
|
|||
|
|
"### 注意事项：\n",
|
|||
|
|
"- 计算结果基于理论公式，实际应用中需要考虑地质复杂性\n",
|
|||
|
|
"- 网格分辨率会影响计算精度和显示效果\n",
|
|||
|
|
"- 3D可视化对计算资源要求较高，建议使用较小的网格"
|
|||
|
|
]
|
|||
|
|
}
|
|||
|
|
],
|
|||
|
|
"metadata": {
|
|||
|
|
"kernelspec": {
|
|||
|
|
"display_name": "Python (forward_modeling)",
|
|||
|
|
"language": "python",
|
|||
|
|
"name": "forward_modeling"
|
|||
|
|
},
|
|||
|
|
"language_info": {
|
|||
|
|
"codemirror_mode": {
|
|||
|
|
"name": "ipython",
|
|||
|
|
"version": 3
|
|||
|
|
},
|
|||
|
|
"file_extension": ".py",
|
|||
|
|
"mimetype": "text/x-python",
|
|||
|
|
"name": "python",
|
|||
|
|
"nbconvert_exporter": "python",
|
|||
|
|
"pygments_lexer": "ipython3",
|
|||
|
|
"version": "3.12.12"
|
|||
|
|
}
|
|||
|
|
},
|
|||
|
|
"nbformat": 4,
|
|||
|
|
"nbformat_minor": 4
|
|||
|
|
}
|