拓冰网站定制 
简介这份资源围绕重力测量中的地形改正技术展开面向地球科学、遥感与GIS领域的学习者和研究人员帮助解决地形起伏导致重力异常难以精确解析的问题。压缩包共4个文件以txt数据文件和m脚本为主包含高程与重力观测数据及地形改正计算程序整体约1KB结构精简便于直接上手运行与验证。已有403人学习下载说明该方向具备一定关注度。读者可借助其中的数据与脚本理解从DEM高程数据准备、理论重力计算到地形改正项求解的完整流程并对照积分法、差分法及Kolosov、Harmon等近似模型的实现思路掌握滤波平滑、改正项合并与结果验证等关键环节。资源虽小但覆盖了地形改正的核心数据与代码骨架适合作为重力数据处理入门练习或课程实验的参考素材也可在此基础上扩展大气压力改正、仪器零点改正等后续处理步骤。1. 地形改正到底在改什么从重力异常里的“假信号”说起做过重力测量的人都有个体会同一台相对重力仪同一个测点隔天复测读数能对上但把测点从山脚挪到山顶布格异常却像换了个地方。这不是仪器漂了是脚下那堆岩石在捣乱。地形改正terrain correction要干的事就是把测点周围高低起伏的地形对重力观测值的引力影响算出来从观测重力里扣掉让剩下的异常真正反映地下密度不均匀而不是被地表形状带偏。它和中间层改正、自由空气改正一起构成布格改正的完整链条缺了它山区重力资料基本没法用。这套流程适合做区域重力调查、矿产勘查、地热与工程物探的人尤其是测区落差超过几十米、地形破碎的场景——平原区可以偷懒山区偷懒就是给自己埋雷。2. 地形改正的物理模型与网格离散从积分公式到可算的扇形环2.1 从牛顿引力出发把地形切成能算的块地形改正的物理起点很朴素测点周围每一块地形质量都对测点有引力把所有块的垂直分量加起来就是地形对重力的影响。设测点位于原点周围地形用高程表示经典做法是把水平面按极坐标分成若干扇形环每个环再按方位角分成若干扇形块每块用平均高程近似成一个直立棱柱或圆柱体。对第 i 个扇形块其垂直引力分量近似为Δg_i G·ρ·Δθ·(r₂ - r₁)·[1/√(r² h²)] 形式的高程相关项实际工程里更常用的是扇形分区模板把每个扇形块的高程差代入预先算好的系数表。这里的关键参数有三个扇形环数、方位角分块数、最大影响半径。环数决定近区精度方位分块决定方向分辨率最大半径决定你截断到多远就不再算。常见做法是近区用密环、远区用疏环半径从几米到几十公里分几段处理。提示地形改正的符号约定必须和你的布格改正流程统一。有的软件定义地形改正为“扣除地形引力”有的定义为“加回”混用会让布格异常整体偏移。2.2 网格离散DEM 分辨率直接决定改正量的可信度现代地形改正基本不再手工查表而是用 DEM 做数值积分。流程是以测点为中心按扇形环和方位角生成采样点从 DEM 读取每个采样点高程计算该块与测点的高差和水平距离代入棱柱公式求和。DEM 分辨率是第一个要盯的参数。如果 DEM 网格比你的近区扇形块还粗近区改正就是插值出来的假数。经验上近区050 米最好用 1 米或 5 米 DEM中区50500 米用 1030 米 DEM远区可以用 90 米 DEM。第二个参数是密度。地形改正对密度是线性敏感的密度取 2.67 g/cm³ 还是 2.4 g/cm³改正量能差 10% 以上。山区没有密度实测时常用 2.67 g/cm³ 作为地壳平均密度但沉积区要小心。2.3 用 Python 跑通一个扇形环地形改正的最小实现下面这段代码演示单点地形改正的核心计算给定测点坐标、DEM 数组、密度和扇形参数输出改正量。它不依赖 GIS 库方便你嵌进自己的处理链。import numpy as np def terrain_correction(x0, y0, dem, dx, rho2.67, n_ring8, n_sector16, r_max5000.0): 单点扇形环地形改正 x0, y0: 测点在地面坐标系中的位置米 dem: 二维高程数组单位米 dx: DEM 网格间距单位米 rho: 地形密度g/cm³ n_ring: 扇形环数对数间隔 n_sector: 方位角分块数 r_max: 最大影响半径单位米 返回地形改正量单位 mGal G 6.67430e-11 # 引力常数 rho_si rho * 1000.0 # g/cm³ - kg/m³ # 对数间隔的环半径近密远疏 r_edges np.logspace(np.log10(1.0), np.log10(r_max), n_ring 1) dtheta 2 * np.pi / n_sector total 0.0 for i in range(n_ring): r1, r2 r_edges[i], r_edges[i 1] r_mid 0.5 * (r1 r2) for j in range(n_sector): theta (j 0.5) * dtheta # 采样点水平坐标 xs x0 r_mid * np.cos(theta) ys y0 r_mid * np.sin(theta) # 最近邻取高程实际项目建议双线性插值 ix int(round(xs / dx)) iy int(round(ys / dx)) if ix 0 or iy 0 or ix dem.shape[1] or iy dem.shape[0]: continue h dem[iy, ix] - dem[int(round(y0/dx)), int(round(x0/dx))] # 棱柱垂直引力分量近似 dr r2 - r1 dg G * rho_si * dtheta * dr * h / np.sqrt(r_mid**2 h**2) total dg return total * 1e5 # SI - mGal逻辑说明外层循环遍历扇形环内层遍历方位角块每个块取环中点作为代表距离用最近邻从 DEM 取高程。h是该块相对测点的高差dg是垂直引力分量的近似。最后乘 1e5 把 SI 单位转成 mGal。参数上n_ring和n_sector越大越精细但计算量按乘积增长r_max取 5000 米是常见中远区截断实际要看测区地形和精度要求。rho默认 2.67沉积区要换成实测值。这段代码是教学级实现生产环境要把最近邻换成双线性插值并处理测点本身高程。3. 从单点到整幅图批量地形改正的工程化流程3.1 数据准备DEM 拼接、投影统一与测点高程校正单点算通了整幅图还有一堆工程问题。第一步是 DEM 拼接。测区跨多幅 DEM 时先做无缝拼接再统一投影到以测区中心为原点的平面坐标系单位用米。第二步是测点高程校正。重力测点的高程通常来自 GPS 或水准和 DEM 在测点位置的高程可能有几米差异。这个差异会直接进入近区改正因为近区对高差最敏感。常见做法是用测点实测高程替换 DEM 在测点位置的值或者对近区采样点做高程改正。第三步是密度模型。如果测区有密度测井或岩石标本数据按地质单元分区给密度没有就统一用 2.67但在报告里写明。3.2 批量计算用多进程把整幅图跑完整幅图可能有几千到几万个测点单点循环太慢。下面用multiprocessing做并行把每个测点的地形改正独立计算。import numpy as np from multiprocessing import Pool def worker(args): x0, y0, dem, dx, rho, n_ring, n_sector, r_max args return terrain_correction(x0, y0, dem, dx, rho, n_ring, n_sector, r_max) def batch_terrain_correction(points, dem, dx, rho2.67, n_ring8, n_sector16, r_max5000.0, n_proc8): points: [(x, y), ...] 测点平面坐标 返回: 每个测点的地形改正量列表 tasks [(x, y, dem, dx, rho, n_ring, n_sector, r_max) for x, y in points] with Pool(n_proc) as p: results p.map(worker, tasks) return results逻辑说明worker把单点计算包装成可序列化的任务batch_terrain_correction用进程池分发。n_proc按机器核数设一般取核数的 70%80%留出内存带宽。注意 DEM 数组会被每个进程复制如果 DEM 很大改用共享内存或把 DEM 切成测区子块。参数上n_ring和n_sector在批量阶段可以先小后大先用粗参数跑一遍看量级再对重点测点加密。3.3 结果检查用剖面和统计量抓出异常改正量批量跑完不能直接信。先做三件事一是画地形改正量等值线图看是否和地形起伏正相关如果出现和地形无关的条带多半是 DEM 拼接缝或投影问题二是沿一条已知剖面检查改正量随距离的衰减正常应该随半径增大快速衰减三是统计改正量的均值和标准差山区典型值在 0.110 mGal 量级如果出现几十 mGal 的点回去查该点近区 DEM 是否有建筑或陡崖。下面这段代码做基本统计和异常点标记。import numpy as np def qc_terrain_correction(values, points, threshold3.0): 对地形改正量做统计质检 values: 改正量列表 points: 对应测点坐标 threshold: 标准差倍数阈值 返回: 异常点索引列表 arr np.array(values) mu, sigma arr.mean(), arr.std() print(f均值{mu:.3f} mGal, 标准差{sigma:.3f} mGal) outliers np.where(np.abs(arr - mu) threshold * sigma)[0] for idx in outliers: print(f异常点 {points[idx]}, 改正量{arr[idx]:.3f} mGal) return outliers逻辑说明用均值和标准差做粗筛超过 3 倍标准差的点标出来人工复核。threshold可以按测区经验调地形特别破碎时放宽到 4。这个质检不能替代逐点检查但能快速定位明显错误。4. 地形改正的避坑与排查五个让结果翻车的细节4.1 近区高程用错测点高程和 DEM 不一致现象近区改正量异常大或异常小和相邻测点对不上。原因测点实测高程和 DEM 在测点位置的高程差了几米近区扇形块的高差被放大。解决批量计算前先算每个测点实测高程与 DEM 高程的差值超过 2 米的点单独处理用实测高程替换 DEM 局部值或在近区计算时对高差做校正。4.2 密度取错沉积区用 2.67 导致系统偏差现象整幅图的地形改正量整体偏大布格异常出现和地形相关的长波长假异常。原因沉积盆地地形密度实际在 2.22.4 g/cm³用 2.67 会高估改正量。解决有密度测井就用测井值没有就按地质图分区给密度至少把沉积区和基岩区分开。报告里写明密度来源和不确定性。4.3 最大半径截断太早远区地形贡献被丢掉现象山区大范围地形起伏的测点改正量偏小布格异常里残留和区域地形相关的趋势。原因r_max只取了几公里远区地形引力没算进去。解决做一次截断试验把r_max从 5 公里逐步加到 20 公里、50 公里看改正量收敛情况。一般到 2030 公里就基本收敛但具体看地形。4.4 DEM 拼接缝和投影变形改正量出现条带现象改正量等值线图上出现直线条带和地形无关。原因多幅 DEM 拼接时高程基准不一致或投影变形导致距离计算偏差。解决拼接前统一高程基准投影用测区中心经线的高斯投影或 UTM检查拼接缝两侧的高程连续性。投影变形大的测区距离计算要用实际地面距离而不是平面距离。4.5 符号约定混乱改正量加反了现象布格异常和自由空气异常的趋势完全相反或者改正后异常更乱。原因地形改正的符号和中间层改正、自由空气改正的符号约定不统一。解决在流程文档里写死符号约定用一个人工算例验证一个孤立山包测点在山脚地形改正应该是把山包的引力扣掉改正量符号要能让布格异常回到区域背景。5. 进阶用自适应分区和交叉验证把地形改正做到可信地形改正做到能跑不难做到可信要花功夫。我一般会加两步自适应分区和交叉验证。自适应分区是根据测点周围地形起伏动态调整扇形环和方位角分块。地形平坦的方向少分块陡崖方向多分块。实现上先算测点周围 8 个方向的坡度坡度大的方向把n_sector加密坡度小的方向减半。这样在保证精度的同时省计算量。交叉验证是换一套 DEM 或换一组扇形参数重算比较两次改正量的差异。差异小于 0.05 mGal 的测点可以放心用差异大的点回去查 DEM 和参数。下面是一个自适应分块的简化实现按方向坡度调整方位角采样密度。import numpy as np def adaptive_sectors(x0, y0, dem, dx, base_sector16, slope_thresh15.0): 根据8方向坡度自适应确定方位角分块数 返回: 每个方向的分块数列表 dirs np.linspace(0, 2*np.pi, 8, endpointFalse) sectors [] for theta in dirs: # 沿该方向取3个采样点估坡度 hs [] for r in [50, 100, 200]: ix int(round((x0 r*np.cos(theta)) / dx)) iy int(round((y0 r*np.sin(theta)) / dx)) if 0 ix dem.shape[1] and 0 iy dem.shape[0]: hs.append(dem[iy, ix]) if len(hs) 2: sectors.append(base_sector) continue slope abs(hs[-1] - hs[0]) / 150.0 * 100 # 百分比坡度 if slope slope_thresh: sectors.append(base_sector * 2) else: sectors.append(base_sector) return sectors逻辑说明沿 8 个方向各取 3 个距离的高程用首尾高差估坡度坡度超过阈值就把该方向的分块数加倍。slope_thresh默认 15%山区可以降到 10%。这个函数返回每个方向的分块数传给主计算函数替换固定的n_sector。注意这只是一个方向性加密实际项目还要考虑环半径的自适应。交叉验证的表格可以这样组织验证项方案 A方案 B差异阈值处理DEM 来源1m 近区 30m 远区5m 全区0.05 mGal超阈值查近区密度2.67 统一分区密度0.1 mGal超阈值查地质图最大半径5 km20 km0.05 mGal超阈值扩半径扇形参数8 环 16 扇16 环 32 扇0.03 mGal超阈值加密这张表是我自己跑项目时的检查清单每次换测区先跑一遍心里有数再批量。地形改正没有绝对正确的值只有在你给定的 DEM、密度和参数下自洽的值。把不确定性量化出来比追求一个“精确”数字更有意义。希望帮到你。本文还有配套的精品资源点击获取 重力地形改正原理与Python实现：从扇形环离散到批量工程化
简介这份资源围绕重力测量中的地形改正技术展开面向地球科学、遥感与GIS领域的学习者和研究人员帮助解决地形起伏导致重力异常难以精确解析的问题。压缩包共4个文件以txt数据文件和m脚本为主包含高程与重力观测数据及地形改正计算程序整体约1KB结构精简便于直接上手运行与验证。已有403人学习下载说明该方向具备一定关注度。读者可借助其中的数据与脚本理解从DEM高程数据准备、理论重力计算到地形改正项求解的完整流程并对照积分法、差分法及Kolosov、Harmon等近似模型的实现思路掌握滤波平滑、改正项合并与结果验证等关键环节。资源虽小但覆盖了地形改正的核心数据与代码骨架适合作为重力数据处理入门练习或课程实验的参考素材也可在此基础上扩展大气压力改正、仪器零点改正等后续处理步骤。1. 地形改正到底在改什么从重力异常里的“假信号”说起做过重力测量的人都有个体会同一台相对重力仪同一个测点隔天复测读数能对上但把测点从山脚挪到山顶布格异常却像换了个地方。这不是仪器漂了是脚下那堆岩石在捣乱。地形改正terrain correction要干的事就是把测点周围高低起伏的地形对重力观测值的引力影响算出来从观测重力里扣掉让剩下的异常真正反映地下密度不均匀而不是被地表形状带偏。它和中间层改正、自由空气改正一起构成布格改正的完整链条缺了它山区重力资料基本没法用。这套流程适合做区域重力调查、矿产勘查、地热与工程物探的人尤其是测区落差超过几十米、地形破碎的场景——平原区可以偷懒山区偷懒就是给自己埋雷。2. 地形改正的物理模型与网格离散从积分公式到可算的扇形环2.1 从牛顿引力出发把地形切成能算的块地形改正的物理起点很朴素测点周围每一块地形质量都对测点有引力把所有块的垂直分量加起来就是地形对重力的影响。设测点位于原点周围地形用高程表示经典做法是把水平面按极坐标分成若干扇形环每个环再按方位角分成若干扇形块每块用平均高程近似成一个直立棱柱或圆柱体。对第 i 个扇形块其垂直引力分量近似为Δg_i G·ρ·Δθ·(r₂ - r₁)·[1/√(r² h²)] 形式的高程相关项实际工程里更常用的是扇形分区模板把每个扇形块的高程差代入预先算好的系数表。这里的关键参数有三个扇形环数、方位角分块数、最大影响半径。环数决定近区精度方位分块决定方向分辨率最大半径决定你截断到多远就不再算。常见做法是近区用密环、远区用疏环半径从几米到几十公里分几段处理。提示地形改正的符号约定必须和你的布格改正流程统一。有的软件定义地形改正为“扣除地形引力”有的定义为“加回”混用会让布格异常整体偏移。2.2 网格离散DEM 分辨率直接决定改正量的可信度现代地形改正基本不再手工查表而是用 DEM 做数值积分。流程是以测点为中心按扇形环和方位角生成采样点从 DEM 读取每个采样点高程计算该块与测点的高差和水平距离代入棱柱公式求和。DEM 分辨率是第一个要盯的参数。如果 DEM 网格比你的近区扇形块还粗近区改正就是插值出来的假数。经验上近区050 米最好用 1 米或 5 米 DEM中区50500 米用 1030 米 DEM远区可以用 90 米 DEM。第二个参数是密度。地形改正对密度是线性敏感的密度取 2.67 g/cm³ 还是 2.4 g/cm³改正量能差 10% 以上。山区没有密度实测时常用 2.67 g/cm³ 作为地壳平均密度但沉积区要小心。2.3 用 Python 跑通一个扇形环地形改正的最小实现下面这段代码演示单点地形改正的核心计算给定测点坐标、DEM 数组、密度和扇形参数输出改正量。它不依赖 GIS 库方便你嵌进自己的处理链。import numpy as np def terrain_correction(x0, y0, dem, dx, rho2.67, n_ring8, n_sector16, r_max5000.0): 单点扇形环地形改正 x0, y0: 测点在地面坐标系中的位置米 dem: 二维高程数组单位米 dx: DEM 网格间距单位米 rho: 地形密度g/cm³ n_ring: 扇形环数对数间隔 n_sector: 方位角分块数 r_max: 最大影响半径单位米 返回地形改正量单位 mGal G 6.67430e-11 # 引力常数 rho_si rho * 1000.0 # g/cm³ - kg/m³ # 对数间隔的环半径近密远疏 r_edges np.logspace(np.log10(1.0), np.log10(r_max), n_ring 1) dtheta 2 * np.pi / n_sector total 0.0 for i in range(n_ring): r1, r2 r_edges[i], r_edges[i 1] r_mid 0.5 * (r1 r2) for j in range(n_sector): theta (j 0.5) * dtheta # 采样点水平坐标 xs x0 r_mid * np.cos(theta) ys y0 r_mid * np.sin(theta) # 最近邻取高程实际项目建议双线性插值 ix int(round(xs / dx)) iy int(round(ys / dx)) if ix 0 or iy 0 or ix dem.shape[1] or iy dem.shape[0]: continue h dem[iy, ix] - dem[int(round(y0/dx)), int(round(x0/dx))] # 棱柱垂直引力分量近似 dr r2 - r1 dg G * rho_si * dtheta * dr * h / np.sqrt(r_mid**2 h**2) total dg return total * 1e5 # SI - mGal逻辑说明外层循环遍历扇形环内层遍历方位角块每个块取环中点作为代表距离用最近邻从 DEM 取高程。h是该块相对测点的高差dg是垂直引力分量的近似。最后乘 1e5 把 SI 单位转成 mGal。参数上n_ring和n_sector越大越精细但计算量按乘积增长r_max取 5000 米是常见中远区截断实际要看测区地形和精度要求。rho默认 2.67沉积区要换成实测值。这段代码是教学级实现生产环境要把最近邻换成双线性插值并处理测点本身高程。3. 从单点到整幅图批量地形改正的工程化流程3.1 数据准备DEM 拼接、投影统一与测点高程校正单点算通了整幅图还有一堆工程问题。第一步是 DEM 拼接。测区跨多幅 DEM 时先做无缝拼接再统一投影到以测区中心为原点的平面坐标系单位用米。第二步是测点高程校正。重力测点的高程通常来自 GPS 或水准和 DEM 在测点位置的高程可能有几米差异。这个差异会直接进入近区改正因为近区对高差最敏感。常见做法是用测点实测高程替换 DEM 在测点位置的值或者对近区采样点做高程改正。第三步是密度模型。如果测区有密度测井或岩石标本数据按地质单元分区给密度没有就统一用 2.67但在报告里写明。3.2 批量计算用多进程把整幅图跑完整幅图可能有几千到几万个测点单点循环太慢。下面用multiprocessing做并行把每个测点的地形改正独立计算。import numpy as np from multiprocessing import Pool def worker(args): x0, y0, dem, dx, rho, n_ring, n_sector, r_max args return terrain_correction(x0, y0, dem, dx, rho, n_ring, n_sector, r_max) def batch_terrain_correction(points, dem, dx, rho2.67, n_ring8, n_sector16, r_max5000.0, n_proc8): points: [(x, y), ...] 测点平面坐标 返回: 每个测点的地形改正量列表 tasks [(x, y, dem, dx, rho, n_ring, n_sector, r_max) for x, y in points] with Pool(n_proc) as p: results p.map(worker, tasks) return results逻辑说明worker把单点计算包装成可序列化的任务batch_terrain_correction用进程池分发。n_proc按机器核数设一般取核数的 70%80%留出内存带宽。注意 DEM 数组会被每个进程复制如果 DEM 很大改用共享内存或把 DEM 切成测区子块。参数上n_ring和n_sector在批量阶段可以先小后大先用粗参数跑一遍看量级再对重点测点加密。3.3 结果检查用剖面和统计量抓出异常改正量批量跑完不能直接信。先做三件事一是画地形改正量等值线图看是否和地形起伏正相关如果出现和地形无关的条带多半是 DEM 拼接缝或投影问题二是沿一条已知剖面检查改正量随距离的衰减正常应该随半径增大快速衰减三是统计改正量的均值和标准差山区典型值在 0.110 mGal 量级如果出现几十 mGal 的点回去查该点近区 DEM 是否有建筑或陡崖。下面这段代码做基本统计和异常点标记。import numpy as np def qc_terrain_correction(values, points, threshold3.0): 对地形改正量做统计质检 values: 改正量列表 points: 对应测点坐标 threshold: 标准差倍数阈值 返回: 异常点索引列表 arr np.array(values) mu, sigma arr.mean(), arr.std() print(f均值{mu:.3f} mGal, 标准差{sigma:.3f} mGal) outliers np.where(np.abs(arr - mu) threshold * sigma)[0] for idx in outliers: print(f异常点 {points[idx]}, 改正量{arr[idx]:.3f} mGal) return outliers逻辑说明用均值和标准差做粗筛超过 3 倍标准差的点标出来人工复核。threshold可以按测区经验调地形特别破碎时放宽到 4。这个质检不能替代逐点检查但能快速定位明显错误。4. 地形改正的避坑与排查五个让结果翻车的细节4.1 近区高程用错测点高程和 DEM 不一致现象近区改正量异常大或异常小和相邻测点对不上。原因测点实测高程和 DEM 在测点位置的高程差了几米近区扇形块的高差被放大。解决批量计算前先算每个测点实测高程与 DEM 高程的差值超过 2 米的点单独处理用实测高程替换 DEM 局部值或在近区计算时对高差做校正。4.2 密度取错沉积区用 2.67 导致系统偏差现象整幅图的地形改正量整体偏大布格异常出现和地形相关的长波长假异常。原因沉积盆地地形密度实际在 2.22.4 g/cm³用 2.67 会高估改正量。解决有密度测井就用测井值没有就按地质图分区给密度至少把沉积区和基岩区分开。报告里写明密度来源和不确定性。4.3 最大半径截断太早远区地形贡献被丢掉现象山区大范围地形起伏的测点改正量偏小布格异常里残留和区域地形相关的趋势。原因r_max只取了几公里远区地形引力没算进去。解决做一次截断试验把r_max从 5 公里逐步加到 20 公里、50 公里看改正量收敛情况。一般到 2030 公里就基本收敛但具体看地形。4.4 DEM 拼接缝和投影变形改正量出现条带现象改正量等值线图上出现直线条带和地形无关。原因多幅 DEM 拼接时高程基准不一致或投影变形导致距离计算偏差。解决拼接前统一高程基准投影用测区中心经线的高斯投影或 UTM检查拼接缝两侧的高程连续性。投影变形大的测区距离计算要用实际地面距离而不是平面距离。4.5 符号约定混乱改正量加反了现象布格异常和自由空气异常的趋势完全相反或者改正后异常更乱。原因地形改正的符号和中间层改正、自由空气改正的符号约定不统一。解决在流程文档里写死符号约定用一个人工算例验证一个孤立山包测点在山脚地形改正应该是把山包的引力扣掉改正量符号要能让布格异常回到区域背景。5. 进阶用自适应分区和交叉验证把地形改正做到可信地形改正做到能跑不难做到可信要花功夫。我一般会加两步自适应分区和交叉验证。自适应分区是根据测点周围地形起伏动态调整扇形环和方位角分块。地形平坦的方向少分块陡崖方向多分块。实现上先算测点周围 8 个方向的坡度坡度大的方向把n_sector加密坡度小的方向减半。这样在保证精度的同时省计算量。交叉验证是换一套 DEM 或换一组扇形参数重算比较两次改正量的差异。差异小于 0.05 mGal 的测点可以放心用差异大的点回去查 DEM 和参数。下面是一个自适应分块的简化实现按方向坡度调整方位角采样密度。import numpy as np def adaptive_sectors(x0, y0, dem, dx, base_sector16, slope_thresh15.0): 根据8方向坡度自适应确定方位角分块数 返回: 每个方向的分块数列表 dirs np.linspace(0, 2*np.pi, 8, endpointFalse) sectors [] for theta in dirs: # 沿该方向取3个采样点估坡度 hs [] for r in [50, 100, 200]: ix int(round((x0 r*np.cos(theta)) / dx)) iy int(round((y0 r*np.sin(theta)) / dx)) if 0 ix dem.shape[1] and 0 iy dem.shape[0]: hs.append(dem[iy, ix]) if len(hs) 2: sectors.append(base_sector) continue slope abs(hs[-1] - hs[0]) / 150.0 * 100 # 百分比坡度 if slope slope_thresh: sectors.append(base_sector * 2) else: sectors.append(base_sector) return sectors逻辑说明沿 8 个方向各取 3 个距离的高程用首尾高差估坡度坡度超过阈值就把该方向的分块数加倍。slope_thresh默认 15%山区可以降到 10%。这个函数返回每个方向的分块数传给主计算函数替换固定的n_sector。注意这只是一个方向性加密实际项目还要考虑环半径的自适应。交叉验证的表格可以这样组织验证项方案 A方案 B差异阈值处理DEM 来源1m 近区 30m 远区5m 全区0.05 mGal超阈值查近区密度2.67 统一分区密度0.1 mGal超阈值查地质图最大半径5 km20 km0.05 mGal超阈值扩半径扇形参数8 环 16 扇16 环 32 扇0.03 mGal超阈值加密这张表是我自己跑项目时的检查清单每次换测区先跑一遍心里有数再批量。地形改正没有绝对正确的值只有在你给定的 DEM、密度和参数下自洽的值。把不确定性量化出来比追求一个“精确”数字更有意义。希望帮到你。本文还有配套的精品资源点击获取