MATLAB地形改正实战：基于dixing.m的重力异常处理

简介面向地球科学、遥感与GIS领域学习者的地形改正代码包用于消除地表起伏对重力观测数据的影响适用于布格重力异常计算、地质构造研究与矿产资源勘查等场景。资源内含MATLAB实现算法与配套样例数据压缩包仅1KB共4个文件其中3个txt文件提供测区的高程与重力测量数据1个m文件为地形改正核心脚本结构精简可直接运行或修改参数。目前已有402人学习下载。脚本基于数字高程模型完成地形改正项计算覆盖数据准备、理论重力参照、积分/差分模型选择以及滤波平滑等关键步骤帮助用户快速掌握从原始观测到改正后重力异常的完整处理流程。通过对比改正前后的重力异常图还可辅助识别地下密度异常体适合作为教学演示或科研参考工具尤其适合地球物理、测绘专业的中初级研究者。1. 地形改正到底在改什么先搞清楚这份MATLAB资源解决的事做重力测量的人都知道地形起伏是重力异常分析里最头疼的干扰项。一个山谷和一个山脊哪怕地下构造完全一样地表测到的重力值也能差出好几个毫伽。这导致你辛辛苦苦采回来的数据如果直接画等值线看到的全是地形的影子而不是真正的地质体响应。地形改正就是为了把这种地形影响剥掉。我这次拆的资源里有一个dixing.m脚本和三个文本数据文件专门用来做这件事。它适合正在做重力数据处理、或者刚接触大地测量但被地形改正公式绕晕的从业者。你不用重新发明轮子把数据格式对齐跑通脚本就能输出改正后的重力异常省掉手动积分那堆破事。2. 地形改正原理与模型选择为什么不能靠“减一个常数”打发2.1 地形对重力场的影响机制地形对重力观测值的影响不是单一方向的。站在山顶时脚下多出来的岩体对重力仪有一个向上吸引的分量会抵消一部分正常重力站在山谷时周围高处的地形又产生一个水平方向的引力分量同样改变测量结果。这种影响随着测点与周围地形的高差变化而复杂分布如果简单地把测区平均高程当成参考面再统一减去一个固定值那等于无视了地形的空间变化改正后的异常依然残留大量地形信号。更麻烦的是地形改正的半径范围不是随便定的。理论上的地形影响范围无限远但实际计算时通常只取一个有限半径比如 20 公里、50 公里。半径取小了远区地形的影响被忽略取大了计算量成倍上升。dixing.m里应该就内置了这类半径参数你需要根据自己测区的地形起伏程度来调整。2.2 常用模型对比积分法、差分法与近似模型摘要里提到了积分法、差分法还有 Kolosov 模型、Harmon 模型或 Laplace 模型。实际工程里积分法是最常用的——把 DEM 划分成一个个柱体或棱柱体计算每个柱体对测点的引力然后求和。这个过程在 MATLAB 里就是用循环或矩阵运算实现。差分法则是通过比较不同高度上的重力值来估算地形影响适用于地形平缓的区域但起伏大的山区误差明显。选择哪个模型主要看测区的地形特征和你的计算资源。我一般会先用积分法因为它物理意义清楚代码写起来也不复杂。Kolosov 这类近似模型适合在缺少高精度 DEM 时使用用简单的几何形状代替真实地形但代价是精度受限。如果你手里有 SRTM 或 ASTER 这类全球 DEM 数据直接走积分法是最稳的。下面是dixing.m里常见的一段地形改正核心计算逻辑我按自己的习惯标注了关键参数% 地形改正计算核心段 % data_a.txt: 高程网格数据单位为米 % data_h.txt: 重力观测值单位为毫伽 % data_r.txt: 测点半径范围内的高程参考值 dem load(data_a.txt); % 加载DEM二维矩阵 g_obs load(data_h.txt); % 加载实测重力值 r_max 50; % 地形改正最大半径单位公里可调 rho 2.67; % 地层密度单位g/cm^3常用2.67 G 6.674e-11; % 万有引力常数单位N·m^2/kg^2 % 网格分辨率假设DEM是等间距网格 dx 100; % 网格间距单位米 dy 100; % 初始化改正项 tc zeros(size(dem)); for i 1:size(dem,1) for j 1:size(dem,2) h0 dem(i,j); % 当前测点高程 sum_grav 0; % 遍历周围网格点简化示例实际需要矢量化和边界处理 for ii 1:size(dem,1) for jj 1:size(dem,2) if iii jjj, continue; end dr sqrt(((ii-i)*dx)^2 ((jj-j)*dy)^2); if dr r_max*1000 dh dem(ii,jj) - h0; % 使用平面棱柱体近似公式 sum_grav sum_grav G*rho*dh*dx*dy/(dr^2 dh^2)^1.5; end end end tc(i,j) sum_grav * 1e5; % 转换到毫伽 end end % 改正后的重力异常 g_anomaly g_obs - tc;这段代码的逻辑非常直白把每个网格点作为测点遍历周围所有点计算高差引起的引力变化累加后得到该点的地形改正值。注意我这只是示意实际工程会用矩阵运算或者fft加速两层循环在数据量大时跑到天荒地老。rho和r_max是两个关键参数前者代表地下岩层的平均密度后者代表你考虑多远的范围。密度取错了改正结果会系统性偏移半径取小了远区地形影响漏掉就成系统性误差。2.3 理论重力与改正项的计算流程摘要里提到需要先计算理论重力值。这一步通常用正常重力公式比如 1967 年的大地测量参考系统公式输入纬度和海拔就能得到该点没有地形影响的重力值。dixing.m里很可能把这部分简化了或者直接让你在data_r.txt里提供参考值。我的习惯是先算出正常重力再用实测值减去正常重力得到布格重力异常最后扣除地形改正项。整个流程可以分四步走准备 DEM 和重力观测数据 → 计算正常重力 → 计算地形改正项 → 合并得到余差。每一步都要检查中间数据尤其是地形改正项的数值范围如果出现几百毫伽的大值先别慌看看是不是网格间距填错了。3. 从文件到结果dixing.m、data_a.txt、data_h.txt、data_r.txt到底怎么配合3.1 文件清单与格式约定这资源里四个文件都不是摆设。data_a.txt存的是高程网格我假设它是纯文本矩阵行列分别对应东西和南北方向单位是米。data_h.txt存的是重力观测值单位一般是毫伽需要你确认自己的数据是不是同单位。data_r.txt比较微妙它能是测点的半径参考值也能是理论重力值得看dixing.m里怎么读。我建议你打开文件先瞄一眼如果data_r.txt只有一列且数字都差不多那可能是理论重力如果是二维矩阵那可能是参考高程或者密度模型。dixing.m是主脚本承担了读取、计算、输出全部工作。脚本里没有写文件名硬编码而是通过load函数直接加载那三个文本文件所以你的工作目录里必须同时放这四个文件。文件名不能改动否则脚本会找不到数据报错。3.2 脚本核心逻辑拆解先不急着跑用 MATLAB 的edit命令打开dixing.m把每个变量名的含义理清。我最关心的是三个变量dem、g_obs、tc。dem是高程矩阵g_obs是观测重力值tc是地形改正项。data_r.txt如果作为输入一般会在脚本里被用来做一个关键判断——比如地形修正半径的范围或是作为正常重力参与计算。脚本里大概率有一个关键常数rho表示密度默认值可能是 2.67这是 crustal rock 的平均密度。如果你的测区是玄武岩或者沉积岩区要按实际情况改到 2.8 或 2.4 左右。还有r_max这个半径参数决定参与计算的 DEM 网格范围。3.3 理论重力与改正项的计算流程我拆过不少这类脚本常见的实现顺序如下读取data_a.txt得到高程读取data_h.txt得到实测重力读取data_r.txt得到测点坐标或参考半径。然后计算每个测点的理论重力值正常重力用实测值减去理论值获得布格异常接着调用子函数计算地形改正项tc最后输出改正后的异常值。这里有个关键的物理点地形改正的符号。理论上地形质量对测点的引力向上为主所以改正项通常为负。但脚本里可能直接把它作为减法项比如g_anomaly g_obs - tc那你得检查符号方向对了没有。我见过有些人直接把原始公式抄漏一个负号导致整个异常都反转过来。4. 跑通流程参数设置、执行顺序与输出验证4.1 运行前置条件与参数调整在 MATLAB 里跑dixing.m之前确保当前文件夹里有那四个文件并且工作路径没有中文。直接用run(dixing.m)或点运行按钮。脚本运行结束后工作区会出现tc和g_anomaly变量这就是你需要的改正结果。参数调整主要在脚本开头那几行。密度rho影响改正项的量级一定要按测区岩性改。半径r_max控制参与计算的网格范围地形起伏大的区域建议设到 30 公里以上平原区可以缩小到 10 公里。还有网格间距dx和dy一定要和你实际 DEM 分辨率保持一致否则计算结果完全不可信。4.2 逐步执行与中间结果检查我不建议直接一口气跑到底而是分步执行。第一步只加载数据看看dem、g_obs的尺寸是否一致。data_a.txt和data_h.txt如果都是 100×100 的矩阵那就没问题如果一个是行向量一个是矩阵说明格式没对齐。用size()检查再用min()、max()看数值范围是否合理。第二步是单独算正常重力可以用 MATLAB 自带的重力函数或者自己写个公式。这一步不依赖地形数据纯粹用纬度算。如果算出来的理论重力值比实测值小很多别急着怀疑代码先确认单位换算了没有——毫伽和伽差 1000 倍。第三步算地形改正项tc检查它的数值范围。一个 500 米高差的山测区改正项通常应该在几毫伽到几十毫伽之间。如果出现上千的值十有八九是半径参数写成了公里但代码当米用。输出成图看surf(tc)或者imagesc(tc)能直观看到改正项和地形走势是否一致。4.3 改正效果的评价方法改正后的重力异常g_anomaly要怎么看第一画等值线图跟原g_obs对比看山脊和山谷位置的高值是不是被削平了。第二算一下改正前后异常的方差应该显著减小。如果方差反而变大说明改正方向反了。第三如果有已知地质体位置比如矿体或断层看异常图上是不是有对应的高值或低值带出现这能验证改正结果的地质合理性。我习惯用plot(g_obs, r)和plot(g_anomaly, b)叠加看总体形态变化。如果改正后曲线更平滑说明地形干扰被有效压制了。这个判断不是量化指标但非常直觉适合用来快速排除明显错误。5. 避坑与常见问题我踩过的五个地形改正大坑5.1 数据格式不匹配导致脚本崩溃现象运行dixing.m时 MATLAB 报错提示load的文件维度不一致。原因data_a.txt是 100×100 矩阵data_h.txt却是 1000 行单列脚本里两者直接相减维度不匹配。解决先读文件头部确认所有文件网格行列一致。如果是单列数据写成和 DEM 相同的矩阵形状。比如用reshape(data_h, size(dem))强制转成矩阵。5.2 地形改正项符号反了现象改正前后的重力异常差值巨大而且改正后异常比改正前更混乱等值线图一片花。原因公式里正负号搞反了。大多数情况下在山区重力仪受地形向上牵引读数偏小所以地形改正项应该加回去而不是减掉。我在脚本里看到g_anomaly g_obs - tc时先检查tc的实际正负——如果tc是正数减去它会让山顶的重力值更低这不对。解决改成g_anomaly g_obs tc或者把tc的符号取反具体看代码里公式是哪种约定。5.3 密度参数取值错误导致系统性偏差现象改正结果看起来光滑了但这和区域地质背景不一致比如已知有铁矿异常但图上看不出来。原因密度rho默认 2.67 均匀密度但实际测区地下若有大范围高密度体这显然不适用。解决查看测区地质图换用更合理的平均密度值。没有地质资料时可以用地形改正迭代反演法先解出一个最佳密度再用于正式改正。5.4 边缘效应把边界搞出巨大假异常现象改正后的异常图四周出现一圈极端的正负相间值而且离边界越近越夸张。原因DEM 和测区边界不重合边缘处的网格点只能计算区域内部分地形影响远区数据缺失形成人为截断。解决扩展 DEM 范围到测区外至少 5 公里用外推或填充背景高程。另外把测点区域适当缩在 DEM 内部别让测点靠近边界。5.5 坐标系单位混用让改正值离谱现象改正项数值大得离谱比如几十万毫伽一眼看出不科学。原因经纬度坐标和投影坐标混用或者把度当成米来计算距离。解决统一用投影坐标单位比如 UTM 坐标距离自然就是米。如果非要经纬度用deg2km函数做换算。我在dixing.m里会特意检查dx和dy确保它们是以米为单位而不是度。6. 进阶用法把脚本改造成自己的重力处理工具dixing.m虽说是一份现成资源但如果只学会点运行那有点浪费。我给你的进阶建议一是加入大气改正和零点漂移改正把它们作为独立子函数挂到主脚本后面这样你得到的最终异常值就完整了。二是批量化处理多测区、多条测线数据把脚本改造成函数形式输入文件名和参数输出改正后的异常后续混合地质解译会省力很多。一个非常实用的技巧是给脚本加一个输出参数控制比如增加一个save_result开关运行时选择是否把g_anomaly保存为 .txt 或 .mat。这样就不用来回在工作区复制数据。我一般是这么改function g_res terrain_correction(dem_file, obs_file, ref_file, rho, r_max) % 地形改正函数封装 dem load(dem_file); g_obs load(obs_file); % ... 中间计算过程略 ... g_res g_obs - tc; if nargin 5 strcmp(save_path, none) % 不保存 else save(g_result.txt, g_res, -ascii); end end封装的好处不言自明你能重复利用不会把数据一路写死在脚本里更不会因为换了个测区就要改源码。要验证结果最靠谱的方法是把改正后的异常与区域内已知地质剖面叠加看灰梯度隆起区和重力高是否对应。如果你手里没有这类地质资料退而求其次的做法是把改正前后绘成三图对比地形图、原始异常图、改正异常图交给有经验的地球物理解释人员判读。最后讲一个我的个人习惯从那次把密度参数写错导致整条剖面解释返工之后我每次跑地形改正都强制走一遍“先小范围、后全测区”的验证流程。先用 5×5 个点的小窗口跑一遍打印出中间变量确认符号、量级、分布都正常再放开全量数据。这么做虽然看起來多花五分钟但能省掉后面改参数重跑重图的几个小时。希望帮到你也祝你的重力异常图一次就能“削平”地形干扰。本文还有配套的精品资源点击获取