企业官网建设流程全解析

简介：这份资源聚焦重力测量中的地形改正处理，面向地球科学、遥感与GIS方向的学习者和研究人员，帮助解决地形起伏导致重力异常难以准确解析的问题。内容围绕高程数据准备、理论重力计算、地形改正模型选择、改正项计算与结果验证等环节展开，涉及积分法、差分法以及Kolosov、Harmon、Laplace等近似模型的思路，适合需要理解地形改正原理并动手处理重力数据的中级学习者。压缩包共4个文件，以3个txt数据文件和1个m脚本为主，整体约1KB，体积轻量，便于快速查看数据组织方式与算法实现框架。目前已有403人学习下载，可作为重力地形改正入门练习与代码参考的起点，帮助读者建立从DEM到重力异常改正的完整处理认知。

1. 地形改正到底在改什么：从一组重力数据说起

如果你手头有一批重力测点数据，高程落差几百米，布格异常图上总是出现和地形高度强相关的条带，那大概率不是仪器漂移，而是地形改正没做干净。地形改正（Terrain Correction）在重力勘探里属于中间层改正之后、布格改正之前绕不开的一步，它修正的是测点周围起伏地形对重力观测值的引力贡献。很多人第一次接触会觉得它和中间层改正重复，其实中间层改正只处理测点到基准面之间的规则板，而地形改正处理的是板外那些山包和沟谷——山体多出来的质量往上拉，沟谷缺失的质量往上“松”，两者符号相反，量级在山区能到几个毫伽，足以掩盖一个局部构造异常。这份地形改正.zip给的就是一套最小可复现的脚本加数据：dixing.m是主计算程序，data_a.txt、data_h.txt、data_r.txt分别对应角度、高程、半径或密度相关的输入，适合刚接手重力数据处理、需要快速跑通地形改正流程的人。

2. 拆开压缩包：dixing.m 与三个数据文件的角色分工

2.1 从文件名反推数据组织方式

拿到一个没有 README 的压缩包，第一步不是急着运行，而是看文件名和变量命名。dixing.m是 MATLAB 脚本，dixing就是“地形”的拼音，说明作者大概率是中文语境下做重力的从业者。三个 txt 文件里，data_h.txt最直观，h 对应 height，存的是高程；data_a.txt的 a 可能是 angle，也可能是 area，需要打开看列数；data_r.txt的 r 可能是 radius，也可能是 density 的某种缩写，但密度一般用 rho，所以半径或水平距离的概率更大。常见的地形改正实现里，会把测点周围按方位角和距离划分成扇形环块，每个环块的平均高程从 DEM 采样，然后对每个块做积分求和。如果data_a.txt是方位角、data_r.txt是径向距离，那这套数据就是按极坐标网格组织的，和dixing.m里的双重循环能对上。

打开文件先看维度，用 bash 下几条命令就能摸清结构：

# 查看每个文件的行数和前几行，判断是单列还是多列 wc -l data_a.txt data_h.txt data_r.txt head -5 data_a.txt head -5 data_h.txt head -5 data_r.txt

如果data_a.txt和data_r.txt都是单列，行数等于方位角划分数和径向划分数，那它们定义的是一个扇形网格的坐标轴；data_h.txt的行数应该等于两者行数的乘积，或者是一个二维矩阵按行展开。这一步不做，后面读数据必翻车。

2.2 地形改正的积分逻辑与代码骨架

地形改正的经典公式是對每个质量元求引力在垂直方向的分量，再对整个地形体积积分。实际编程不会真做三重积分，而是把地形离散成若干个棱柱或扇形块，用近似公式累加。dixing.m里大概率是类似这样的结构：

% 读取三个数据文件 a = load('data_a.txt'); % 方位角划分，单位可能是度或弧度 h = load('data_h.txt'); % 各网格点高程，单位米 r = load('data_r.txt'); % 径向距离划分，单位米 % 预分配改正值数组 tc = zeros(size(h)); % 双重循环遍历每个扇形块 for i = 1:length(a)-1 for j = 1:length(r)-1 % 取环块平均方位角与内外半径 az = (a(i) + a(i+1)) / 2; r1 = r(j); r2 = r(j+1); % 取该块平均高程，h 的索引方式取决于数据排列 h_avg = mean(h( ... )); % 计算该块对测点的垂直引力贡献 % 常用公式：G * rho * dA * h / (r^2 + h^2)^(3/2) 的积分近似 dA = 0.5 * (r2^2 - r1^2) * (a(i+1) - a(i)); % 扇形环面积 tc = tc + 6.674e-11 * 2670 * dA * h_avg / ... ( ( (r1+r2)/2 )^2 + h_avg^2 )^(3/2); end end % 单位换算：从 m/s^2 转到 mGal tc = tc * 1e5;

上面这段是骨架，不是dixing.m的原文，但逻辑八九不离十。关键参数有三个：密度rho一般取 2.67 g/cm³（地壳平均密度），面积元dA由方位角增量和半径增量决定，高程h_avg是该块的平均海拔。代码里用6.674e-11是万有引力常数，单位是 m³/(kg·s²)，密度要换算成 kg/m³，所以 2.67 g/cm³ 要写成 2670。最后乘1e5是把 m/s² 转成 mGal，因为重力异常通常用毫伽表示。如果你的dixing.m里没有这个换算，出来的数会小五个数量级，看着像零，其实是单位没转。

2.3 数据文件与代码的对接方式

data_h.txt的排列方式决定了索引怎么写。常见有两种：一种是按方位角优先，即先遍历所有方位角，再换半径，这样h的索引是(i-1)*nr + j；另一种是按半径优先，索引反过来。打开文件看前几行，如果数值随行号增加呈现周期性起伏，那可能是按角度排的。更稳妥的办法是在 MATLAB 里reshape一下：

% 假设方位角有 na 个，半径有 nr 个 na = length(a); nr = length(r); % 如果 h 是 na*nr 个元素，尝试按列优先重塑 h_grid = reshape(h, na, nr); % 然后取块平均时用 h_grid(i,j) 和 h_grid(i+1,j+1) 做平均

如果reshape报错说元素个数不匹配，说明h的维度不是na*nr，可能是(na-1)*(nr-1)，因为环块数比节点数少一个。这时候要改成reshape(h, na-1, nr-1)。这一步没有文档就只能试，试的时候先拿一小段数据跑，看输出量级对不对——山区地形改正一般在 0.1 到 5 mGal 之间，如果算出几百 mGal，肯定是半径单位用错了，把米当成千米或者反过来。

3. 跑通第一遍：从数据检查到结果验证的完整操作

3.1 环境准备与最小运行步骤

MATLAB 版本没有特殊要求，R2016a 以上都能跑，因为脚本里没用到太新的函数。如果你没有 MATLAB，Octave 也能凑合，但load和reshape的行为基本一致，注意 Octave 里mean对矩阵默认按列，和 MATLAB 相同。把压缩包解压到一个纯英文路径下，避免中文路径导致load读不到文件。然后按顺序执行：

% 第一步：清空环境，避免旧变量干扰 clear; clc; close all; % 第二步：检查文件是否存在 assert(exist('data_a.txt','file')==2, '缺少 data_a.txt'); assert(exist('data_h.txt','file')==2, '缺少 data_h.txt'); assert(exist('data_r.txt','file')==2, '缺少 data_r.txt'); % 第三步：运行主脚本 run('dixing.m'); % 第四步：查看输出变量 whos % 如果输出里有 tc 或 similar 变量，检查其数值范围 if exist('tc','var') fprintf('地形改正值范围：%.4f 到 %.4f mGal\n', min(tc(:)), max(tc(:))); end

assert那两行是血泪经验，很多人把脚本和数据的相对路径搞错，MATLAB 当前文件夹不是解压目录，load直接报错。加个断言能省去翻来覆去查路径的时间。运行完whos看变量列表，找到改正值数组，用min和max扫一眼量级。如果全是零，检查密度参数是不是写成了 0；如果出现NaN，检查data_h.txt里有没有空行或非数字字符。

3.2 参数敏感性：密度、半径划分与角度分辨率

地形改正对三个参数最敏感：地壳密度、径向划分的起始半径、方位角划分数。密度取 2.67 g/cm³ 是常规做法，但在沉积盆地地区，近地表密度可能只有 2.0 到 2.3，这时候用 2.67 会高估改正量。dixing.m里如果密度是硬编码的，找到那一行改成你工区的实测密度。径向划分的起始半径不能太小，因为近测点的地形块对重力贡献最大，但公式在 r 趋近 0 时会出现奇点，常见做法是从 1 米或 5 米开始，而不是从 0 开始。方位角划分数一般取 8 的倍数，比如 16、32、64，太少会导致方位方向上的积分误差，太多则计算量暴涨。下面这张表是不同参数组合下的经验值：

如果dixing.m里这些参数是写死的，建议改成脚本开头的变量，方便批量测试。改完跑一遍，对比改正前后布格异常图的等值线是否还和地形等高线平行——如果平行度明显下降，说明改正起了作用。

3.3 结果验证：改正前后对比与量级检查

跑出改正值后，不能直接加到布格异常上就完事，得先验证。最简单的办法是画两张图：一张原始布格异常，一张加上地形改正后的布格异常，看地形相关条带是否减弱。MATLAB 里几行就能画：

% 假设 bg 是原始布格异常，tc 是地形改正值 bg_corrected = bg + tc; % 注意符号：地形改正通常加到观测值上 figure; subplot(1,2,1); imagesc(bg); colorbar; title('原始布格异常'); subplot(1,2,2); imagesc(bg_corrected); colorbar; title('地形改正后'); colormap(jet);

符号问题容易搞反。地形改正的定义是：把测点周围多余地形质量去掉，所以改正值要加到观测重力上，使异常更接近真实地下密度分布。如果你加完发现异常和地形相关性更强了，那就是符号反了，改成减。量级检查用统计量：计算改正值和测点高程的相关系数，如果相关系数绝对值大于 0.8，说明改正量主要受高程控制，这是正常的；但如果改正后布格异常和高程的相关系数没降下来，说明改正没做对，回去查半径单位或密度。

4. 避坑与排查：地形改正脚本最容易翻车的五个地方

4.1 现象：改正值全是零或极小 → 原因：单位换算缺失 → 解决：检查 1e5 换算

这是最高频的翻车点。dixing.m里如果直接输出tc而没有乘1e5，结果会是 10^-6 量级，看着像零。打开脚本搜索1e5或100000，如果没有，在最后输出前加上。另外检查密度单位，如果代码里写的是2.67而不是2670，同样会小三个数量级。两个错误叠加，结果会小八个数量级，彻底没法看。

4.2 现象：运行报错“索引超出矩阵维度” → 原因：h 的排列与循环不匹配 → 解决：先 reshape 再取块

data_h.txt如果是按节点排列，而循环按环块索引，就会越界。解决办法是在循环前把h重塑成二维矩阵，用size确认维度。如果h的长度是(na-1)*(nr-1)，就reshape(h, na-1, nr-1)，然后循环里用h_grid(i,j)。不要凭感觉猜，用numel(h)和(na-1)*(nr-1)比一下，相等就重塑，不相等就看是不是na*nr。

4.3 现象：改正后异常出现同心圆状假异常 → 原因：径向划分太粗或起始半径太大 → 解决：加密近场划分

如果径向环块宽度是固定的，比如每 100 米一环，近测点 100 米内的地形变化就被平均掉了，改正后会在测点周围留下同心圆状残差。常见做法是近场用对数间隔，远场用等间距。如果dixing.m不支持变间距，就把起始半径改小，同时增加环数，比如从 1 米到 1000 米分 50 环，每环宽度按几何级数增长。

4.4 现象：结果对密度参数极其敏感 → 原因：近场块高程过大导致非线性放大 → 解决：检查近场高程是否异常

密度是线性因子，正常情况改密度 10%，改正值也变 10%。如果改密度导致改正值变化超过 50%，说明近场某个块的高程异常大，比如把 5000 米的高程误读成 50000 米。打开data_h.txt用max和min扫一眼，正常高程范围在 -500 到 9000 米之间，超出这个范围就是数据错误。

4.5 现象：Octave 下运行结果和 MATLAB 不一致 → 原因：mean 函数对矩阵维度处理不同 → 解决：显式指定维度

MATLAB 的mean对矩阵默认按列求平均，Octave 也是，但如果代码里用了mean(h)而h是向量，两者一致；如果h是矩阵且代码依赖隐式扩展，Octave 可能报错或给出不同结果。稳妥写法是mean(h, 'all')或mean(h(:))，强制按所有元素求平均。另外 Octave 对assert的支持没问题，但run函数在旧版 Octave 里可能不存在，直接用脚本名运行即可。

5. 进阶技巧：把地形改正嵌进批量处理流程

5.1 用函数封装替代脚本，支持多测点循环

dixing.m如果是脚本，每次只能处理一个测点的数据。实际工区有几百个测点，每个点周围的地形块都不同，需要把改正逻辑封装成函数，输入测点坐标、DEM 网格、密度，输出改正值。下面是一个函数骨架：

function tc = terrain_correction(x0, y0, dem, dx, rho) % x0, y0: 测点坐标 % dem: 数字高程模型矩阵 % dx: DEM 网格间距，单位米 % rho: 地壳密度，单位 kg/m^3 % 返回 tc: 地形改正值，单位 mGal G = 6.674e-11; [nr, nc] = size(dem); tc = 0; % 遍历 DEM 每个网格，计算对测点的引力贡献 for i = 1:nr for j = 1:nc % 网格中心相对测点的水平距离 x = (j - 0.5) * dx - x0; y = (i - 0.5) * dx - y0; r = sqrt(x^2 + y^2); if r < dx/2 continue; % 跳过测点所在网格，避免奇点 end h = dem(i,j); % 矩形棱柱引力近似公式 dA = dx^2; tc = tc + G * rho * dA * h / (r^2 + h^2)^(3/2); end end tc = tc * 1e5; % 转 mGal end

这个函数比扇形划分更直接，适合 DEM 是规则网格的情况。调用时把测点坐标传进去，循环所有测点即可。注意r < dx/2那个判断，跳过测点所在网格，否则r=0会导致除零。如果测点不在网格中心，这个近似会有小误差，但地形改正本身精度要求不高，通常可接受。

5.2 用向量化替代双重循环，提速十倍以上

上面的双重循环在 DEM 很大时慢得让人想砸键盘。MATLAB 擅长矩阵运算，把距离和引力贡献全部向量化：

function tc = terrain_correction_fast(x0, y0, dem, dx, rho) G = 6.674e-11; [nr, nc] = size(dem); [X, Y] = meshgrid((1:nc)*dx - x0, (1:nr)*dx - y0); R = sqrt(X.^2 + Y.^2); R(R < dx/2) = NaN; % 标记测点所在网格 H = dem; dA = dx^2; integrand = G * rho * dA * H ./ (R.^2 + H.^2).^(3/2); tc = nansum(integrand(:)) * 1e5; end

meshgrid生成所有网格中心相对测点的坐标，R是距离矩阵，NaN标记要跳过的点，nansum求和时忽略NaN。这段代码在 1000×1000 的 DEM 上比双重循环快几十倍。注意H和R的维度要一致，如果 DEM 不是方阵，meshgrid的输入顺序要对应行列。

5.3 验证批量结果：用剖面线检查改正连续性

批量跑完后，抽一条穿过地形起伏的剖面线，画改正前后的布格异常。如果改正后剖面曲线变得平滑，和地形剖面的相关性明显降低，说明批量处理没问题。我一般会固定抽三条线：一条沿山脊、一条沿沟谷、一条跨断层，三条都通过才算过。从那以后我每次拿到新的重力数据，都强制先跑一遍单点测试，确认量级和符号，再上批量，这个习惯帮我省了至少两次返工。希望帮到你。

本文还有配套的精品资源，点击获取