处理二维地震剖面或者三维勘探数据体的时候常常会遇到一个问题我们看到的原始记录里重要的结构并不是水平的而是带着各种倾角、呈现出某种局部走向。如果只知道“平均值”“方差”这类全局统计量根本描述不了这种空间结构。真正有用的是知道每一个局部位置的边坡朝向——也就是局部斜率local slope或倾角然后顺着这个方向去增强信号、压制噪声。这种“先估计方向、再沿方向滤波”的思路在MATLAB里做起来并不复杂但涉及不少容易踩的坑。这篇文章从一个实际数据处理者的视角把2D/3D局部边坡估计与结构滤波的完整思路捋一遍。从结构的数学定义、结构张量的计算原理到2D和3D的MATLAB落地实现再到结构导向平滑的具体做法都会给到可复现的思路和代码。适合正在做地震数据处理、图像方向性分析或者需要对体数据做结构增强的MATLAB使用者和研究人员。1. 先从“干什么”说起局部边坡到底在描述什么1.1 一个现象数据里的“同相轴”是有走向的早年我第一次处理实际采集的地震剖面时拿到的一整块二维数据最直观的感觉就是数据里面有非常明显的倾斜纹理——某些能量团像波浪一样斜着排布。当时我想要做平滑降噪却总是不小心把这种有方向的纹理也抹掉了。后来才理解这类数据里的核心信息不是每个像素单独的灰度大小而是空间的延展方向。从数学上看对于二维数据 (d(x,t))如果局部存在一个以斜率 (p \Delta t / \Delta x) 变化的平面波分量那么沿着某一个方向数据的数值几乎不变而跨过这个方向数值迅速变化。这里要估计的“局部边坡”就是这个方向或斜率。三维数据也是一样只是从单一斜率扩展为首要倾角和方位角两个角度。对于地震勘探这对应反射层位的视倾角对于图像这很像边缘纹理的局部方向场。1.2 为什么这个估计“局部”而不是全局全局性的倾角扫描比如把整条剖面放到频率-波数域统计主能量方向有个致命问题实际数据几乎不会是单一平面波。浅层、深层可能有不同的倾角甚至断层两侧的结构方向突然变化。你算出一个全局面上的方向可能在大部分区域都不适用。所以这里的核心要义是在每一个空间位置附近的小窗口上做方向分析。常见的做法是用结构张量structure tensor作为局部方向的描述子——它把局部梯度信息的二阶统计量压缩成一个2×22D或3×33D的对称矩阵。然后通过矩阵的特征分解拿到两个关键信息一个是主方向一个是在该方向上信号强度的对比度。这个过程天然是“局部”的因为梯度本身是每个像素的邻域操作而窗口平均又进一步把信息约束在局部范围。1.3 一个直觉检验白噪声没有方向如果数据是随机噪声局部梯度会在任何方向上都杂乱无章对应的两个特征值大小会非常接近方向不稳定。如果数据具有清晰纹理沿着纹理梯度的方向特征值会明显大于垂直方向方向信息也很稳定。因此特征值的相对大小可以当作“方向可信度”来使用这就为后面结构滤波打下了基础——知道哪些地方的方向可用哪些地方只能当噪声处理。2. 结构张量的数学拆解从梯度到主方向2.1 张量构造先算梯度外积再做空间平滑要计算任意点 ((x_0,t_0)) 的结构张量第一步是取该点邻域内的梯度向量 (\mathbf{g} [g_x, g_t]^T)。如果直接使用单个点的梯度外积噪声会非常大必须在一个窗口内做平均。因此结构张量 (J) 的标准定义为[ J \begin{bmatrix} \overline{g_x^2} \overline{g_x g_t} \ \overline{g_x g_t} \overline{g_t^2} \end{bmatrix} ]其中上划线表示空间局部平均。实际实现时通常先分别计算 (g_x^2)、(g_x g_t)、(g_t^2) 三个分量图再对这三个图做高斯平滑而不是先平均梯度。对于3D多一个维度即可得到3×3的矩阵[ J \begin{bmatrix} \overline{g_x^2} \overline{g_x g_y} \overline{g_x g_z} \ \overline{g_x g_y} \overline{g_y^2} \overline{g_y g_z} \ \overline{g_x g_z} \overline{g_y g_z} \overline{g_z^2} \end{bmatrix} ]这种构造方式把局部梯度的二阶信息完整保留了下来。平滑核的大小直接控制“局部”的程度核窗口越大方向估计越稳健但空间分辨率越差。2.2 特征分解为什么特征向量一定是主方向对称矩阵必然可以对角化而且特征向量正交这是一个非常好用的性质。对2×2矩阵特征值 (\lambda_1 \ge \lambda_2) 分别代表梯度能量在特征向量方向上的投影大小。最大特征值对应的特征向量就是局部梯度最集中的方向它旋转90度就是局部结构延伸的方向。在二维地震剖面里若把横轴看成 (x)纵轴看成 (t)最大特征向量 ([e_{1x}, e_{1t}]^T) 偏向哪个轴就直接对应倾角大小。利用特征向量算斜率时需要把方向转成物理角度。常用公式[ \theta 0.5 \cdot \arctan2\left(2 J_{12}, J_{11} - J_{22}\right) ]这是通过解析解直接给出的不需要真的去调用 ( \mathrm{eig} )计算效率非常高。后面第3节给出完整实现。2.3 一致性和各向异性判断方向是否可信特征值的另一个重要用途是计算一致性coherence或者叫各向异性系数[ c \frac{\lambda_1 - \lambda_2}{\lambda_1 \lambda_2} ]取值在0到1之间。理想平面波处 (c) 接近1随机噪声处接近0。这个值在结构滤波里特别关键——低一致性区域说明没有稳定方向此时如果强行按照估计的方向去滤波反而会引入伪结构。所以后来做滤波时我通常会把一致性作为混合系数在结构导向平滑和各向同性平滑之间做线性过渡。3D下情况类似特征值有3个常见的分析重点变为((λ_1-λ_2)/(λ_1λ_2))平面性((λ_2-λ_3)/(λ_2λ_3))线状性/边界性对于地质层面这种二维片状结构通常最大和第二大特征值都比较明显第三个小得多对于断层一侧的边缘第二和第三特征值的关系也会有相应差异。3. 二维MATLAB完整实现从函数到可视化3.1 核心函数用解析解避开循环二维的2×2对称矩阵特征分解完全可以用解析式算出来避免对每个像素循环调用eig。下面这段代码是我常用的一种矢量化实现适用于单精度或双精度的二维数据function [theta, coh, lambda1, lambda2] structure_tensor_2d(img, sigma_grad, sigma_win) % img: 二维灰度数据可能是地震剖面、图像等 % sigma_grad: 梯度计算前的可选预平滑尺度常取 1~3 % sigma_win: 结构张量分量的空间平滑尺度常取 3~10 if sigma_grad 0 img imgaussfilt(img, sigma_grad); end [gx, gt] gradient(img); J11 imgaussfilt(gx .* gx, sigma_win); J12 imgaussfilt(gx .* gt, sigma_win); J22 imgaussfilt(gt .* gt, sigma_win); % 解析特征分解2x2对称矩阵 tmp sqrt((J11 - J22).^2 4 * J12.^2); lambda1 0.5 * (J11 J22 tmp); lambda2 0.5 * (J11 J22 - tmp); theta 0.5 * atan2(2 * J12, J11 - J22); % 梯度主方向 % 注意结构的走向 theta pi/2如果要用结构延伸方向 coh (lambda1 - lambda2) ./ (lambda1 lambda2 eps); end这段代码有几个细节值得说明。首先gradient默认使用中心差分边界附近效果不理想往往需要先做边缘扩展例如padarray(img, [pad, pad], symmetric)算完再裁掉。其次梯度预平滑很重要如果数据本身噪声较强直接对原始数据求梯度估计出的方向会抖动得非常厉害小尺度的高斯预平滑能明显改善稳定性。最后coh分母加eps是为了防止全零区域出现除零。3.2 参数到底怎么取三把尺子实际用下来参数选取基本可以按数据维度和信噪比来定sigma_grad预平滑尺度。信噪比低时我会取3到6信噪比高、希望看清细小结构时可以取0到1。如果想检测更小尺度结构甚至可以不预平滑但后续方向稳定性差。sigma_win张量分量平滑尺度这是最核心的调参项。数据中的结构延伸越长这个值应当越大但要检测断层这类短尺度突变时这个值不能太大否则会把断层“糊”过去。数据范围和物理单位gradient默认按网格步长1计算。如果横纵两个方向的物理采样间距不同比如纵向时间采样2ms、横向道距25m那么算出来的梯度分量在数值上会完全不同。这种情况建议先做各向同性化处理即按物理间距归一化或者直接把采样间距换算进梯度计算。我通常先用一个中等窗口比如5~13像素跑一遍看方向场是否平滑且符合地质常识再根据效果微调。这类算法对参数有一定的“宽容度”不需要精确到某个值。3.3 自检用合成数据验证方向估计写算法最忌讳直接上真实数据因为很难判断方向估计到底准不准。我的习惯是先生成若干个已知倾角的平面波合成数据再叠加高斯噪声用来评估方向估计的误差。% 生成斜率为 p 0.3 的平面波数据 nx 256; nt 256; [x, t] meshgrid(1:nx, 1:nt); p_true 0.3; data sin(2*pi*0.08*(t - p_true*x)); data data 0.3*randn(size(data)); % 估计结构张量 [theta, coh] structure_tensor_2d(data, 1.5, 8); % 将梯度主方向转换为斜率 slope_est tan(theta); % 注意需要根据坐标定义确定符号 mean_slope mean(slope_est(coh 0.8)); fprintf(真实斜率: %.3f, 估计斜率: %.3f\n, p_true, mean_slope);实测效果通常都不错在合成数据上只要窗口参数不是特别离谱估计斜率误差一般能控制在百分之几。真正会出问题的往往不在算法本身而在坐标符号约定。二维gradient(img)的第一个输出是沿第一维的梯度而第一维在MATLAB中对应的是行方向。如果你习惯把横轴作为空间轴需要特别当心方向矩阵的转置关系。建议在合成数据阶段就把符号问题消化掉不然真实数据上很容易出现角度反向的问题。4. 三维扩展从切片思考切换到体块思考4.1 特征分解从2×2变成3×3角度从一个变成两个二维情况只关心一个倾角三维情况下局部结构有了真正的方向层面同时在横向和纵向倾斜需要两个角度倾角和方位角来描述。所以三维结构张量是3×3对称矩阵特征分解不再有一句简单的解析式通常需要写循环调用eig或者用 MATLAB 的页式特征分解函数如是R2021b之后的版本可以用pagemtimes 循环思路。三维里的“一致性”同二维类似但多层特征值的组合可以用来区分不同结构类型。比如层面型结构最大特征值显著其它两个特征值小断层或河道型线性构造则两个特征值接近且大于第三个。许多开源的方向属性算法本质都是在这种三维特征分解上做文章。4.2 MATLAB三维实现多写一个维度三维结构张量的核心代码和二维几乎同构function [v1, coh, lambda] structure_tensor_3d(vol, sigma_grad, sigma_win) % vol: 三维体数据 [nz, ny, nx] 注意维度顺序在MATLAB中是 z, y, x % 返回值: v1 是对应最大特征值的单位特征向量 [vz, vy, vx] 三个分量图 % coh: 一致性系数 % lambda: lambda1, lambda2, lambda3 三个分量图 if sigma_grad 0 vol imgaussfilt3(vol, sigma_grad); end [gz, gy, gx] gradient(vol); J11 imgaussfilt3(gz .* gz, sigma_win); J12 imgaussfilt3(gz .* gy, sigma_win); J13 imgaussfilt3(gz .* gx, sigma_win); J22 imgaussfilt3(gy .* gy, sigma_win); J23 imgaussfilt3(gy .* gx, sigma_win); J33 imgaussfilt3(gx .* gx, sigma_win); [nz, ny, nx] size(vol); v1 zeros(nz, ny, nx, 3); lambda zeros(nz, ny, nx, 3); coh zeros(nz, ny, nx); for iz 1:nz for iy 1:ny for ix 1:nx J [J11(iz,iy,ix), J12(iz,iy,ix), J13(iz,iy,ix); J12(iz,iy,ix), J22(iz,iy,ix), J23(iz,iy,ix); J13(iz,iy,ix), J23(iz,iy,ix), J33(iz,iy,ix)]; [V, D] eig(J, vector); [~, idx] max(D); v1(iz,iy,ix,:) V(:,idx); lambda(iz,iy,ix,:) D; lmax D(idx); lmin min(D); coh(iz,iy,ix) (lmax - lmin) / (lmax lmin eps); end end end end这段代码保持了思路清晰但三维循环确实慢。数据体较大比如500×500×500时三层嵌套循环的代价非常明显。我有两个优化建议数据结构上把三维体转换为二维矩阵批量处理对每个体素构造3×3矩阵并按一定排列展开用向量化的方式得到所有特征向量和特征值。具体做法是把9个分量矩阵分别保存然后对每个体素提取特征分解——在MATLAB中可以使用arrayfun或自定义的MEX函数但最直接有效的是只对感兴趣区域做特征分解而不是遍历整个体。对特别大的数据先把体数据降采样或分块估计方向用粗网格即可。方向场本身是平滑变化没必要在每个体素上都做精细计算这样能节省大量时间。4.3 三维方向的解释和可视化三维方向场不像二维那样直接画个箭头就完了。常用做法切片分析沿某些关键方向切片在二维剖面上叠加该位置的方向矢量。方向着色把三维方向分解成RGB分量用颜色表示方向场。体素箭头每隔若干体素画一个三维箭头适合展示大体走向。quiver3可以直接在三维切片上画方向箭头。有个坑是方向向量的正负是不定的——特征向量和负特征向量对应同一个特征值所以画箭头时需要用邻域一致性来统一方向否则看起来忽正忽反很乱。统一符号最简单的办法固定第一行的符号为正比如让vz 0但这样在水平层面处也会产生跳变需要根据实际应用处理。5. 结构滤波落地顺着结构方向去平滑而不是跨过去5.1 思路转变从各向同性到各向异性拿到局部方向和一致性数据后“结构滤波”就变成了一件非常自然的事。普通高斯平滑是各向同性的在噪声压制的同时会把断层、尖灭等不连续结构磨平。结构导向滤波的思路恰恰相反在结构延伸方向上多平滑在跨结构方向上少平滑甚至不平滑。这样既能大幅压制随机噪声又不会破坏结构边界效果通常比普通中值滤波好看得多。二维实现上我比较喜欢“按角度分桶”的思路因为结构张量得到的倾角是逐点变化的不可能对每个像素单独构造卷积核。实际做法是把估计到的角度量化成 (M) 个离散角度比如每10度为一个桶。对每个离散角度生成一个旋转后的高斯核长轴沿该角度。对原始数据分别做 (M) 次滤波。输出结果时每个像素按它实际角度所在桶直接选用对应滤波结果或者基于一致性做混合。这个思路既避免了“逐像素自定义核”的复杂循环又利用了imfilter的快速卷积。下面是一个示例代码function out structure_filter_2d(img, theta, coh, sigma_along, sigma_cross) % 角度离散化 M 18; % 每10度一个桶 angle_bins linspace(-pi/2, pi/2, M1); filtered zeros([size(img), M]); [gy, gx] size(img); [xx, yy] meshgrid(1:gx, 1:gy); for k 1:M ang angle_bins(k) pi/2; % 结构延伸方向 kernel directional_gaussian_kernel(ang, sigma_along, sigma_cross); filtered(:,:,k) imfilter(img, kernel, symmetric); end % 对每个像素按其角度选择对应滤波结果 [~, idx] min(abs(theta - angle_bins(1:end-1)), [], 3); % 简化按最近角度选桶 out_all zeros(size(img)); for k 1:M mask_k (idx k); out_all(mask_k) filtered(mask_k, k); % 注意这里只是示意逻辑实际需逐像素赋值 end % 结合一致性低一致性区域退化为各向同性平滑或原图 out_iso imgaussfilt(img, sigma_along); out coh .* out_all (1 - coh) .* out_iso; enddirectional_gaussian_kernel可以这样写先用fspecial(gaussian, [len, len], [sigma_cross, sigma_along])得到一个沿行方向延伸的高斯核在MATLAB里第一维会被当作行方向所以长轴沿列方向也就是垂直方向再用imrotate转成任意角度。旋转之后角落会出现零值为了让滤波保持能量守恒通常还需要对该核做归一化。这一步一定要做不然会出现条带明暗变化的伪影。5.2 一致性滤波强度的保护开关一致性系数 (c) 在这里起的作用是保护非结构区。当时我在一个三维体上做结构导向平滑时没有用一致性做混合结果在低信噪比区域出现了明显的“纹理伪造”——本来就全是噪声的地方被方向滤波强行拉出了假结构。加了 (c) 做线性混合之后非常有效地解决了这个问题。实际操作中一致性阈值也可以按百分位选取例如取所有像素一致性最高的30%作为强结构区其余区域完全做各向同性平滑。这个办法在层位清晰但部分被噪声覆盖的地震数据上效果不错。5.3 结构导向中值滤波一种更抗野值的选择高斯型结构导向平滑对高斯噪声效果好但对极端离群值不如中值滤波稳健。改进方案也简单把“结构导向”的思想套在中值滤波上。做法是把沿结构方向一个小窗口内的数据收集起来取中值作为当前点输出窗口很窄但沿结构方向够长所以既保持了不跨越结构又能压制脉冲噪声。实现上如果追求简单可以在方向拉平域处理把局部窗口沿着结构方向旋转到水平在水平窗口内取中值再旋转回来。MATLAB里可以用imrotate或者用插值实现局部坐标变换。虽然性能不如频域方案但胜在容易理解也容易调试。6. 算完方向之后还有几件容易忽略的事6.1 三组方向定义经常把人绕晕结构张量最大特征向量指向梯度方向结构延伸方向与它垂直。在实际做滤波时真正需要的是延伸方向而不是梯度方向。所以上面结构滤波代码里做方向旋转时我给倾角加了 (\pi/2) 的偏移。如果不加你会发现滤波方向反过来了完全沿着结构垂直方向去平滑噪声压不掉结构却全被抹平了。建议在自己的代码里方向定义一目了然地注释清楚theta是梯度主方向还是结构走向最后做一次合成数据验证确认没有符号翻转。6.2 平滑核尺度与数据尺度的匹配问题结构滤波最容易被忽视的是核尺寸和实际结构尺度的匹配。如果地震剖面上的同相轴在纵向占5个采样点而你结构滤波核横向宽度给了50个网格那就算方向正确也会把可以分辨的薄层细节全部抹掉。我的经验是先统计一下目标结构在数据中的大概尺度设置滤波核时让核的长轴和该尺度相当短轴取长轴的1/5以下即可。6.3 边缘效应无论结构张量估计还是方向滤波都会受边界影响。建议统一采用symmetric对称扩展加一小段衰减窗口的策略。坏处是边缘结果相对平滑好处是边界附近不会出现剧烈的起伏伪影。实际项目里如果分析区域很大直接裁掉边界10~20个采样点就可以。6.4 计算负载的务实取舍三维体数据上完整跑一遍逐体素特征分解加自定义方向滤波时间往往难以接受。务实方案是两步走先在粗网格上每48个体素取一个估计结构方向再对方向场做插值加密最后滤波阶段使用插值后的方向场。这样方向估计的计算量能缩小几十倍滤波本身仍然可以充分利用矩阵运算来快速完成。我在处理512×512×256的体数据时用这个策略把单次处理的耗时从十几分钟压缩到两三分钟。个人经验是局部边坡估计和结构滤波这类方法真正卡住人的往往不是算法本身也不是MATLAB代码实现而是“方向定义是否统一”“参数尺度是否匹配数据”这类看似琐碎的地方。每次新拿一块数据先做合成测试、再跑真实数据能省下大量调试时间。把上面这些环节理顺后这套方法无论放在图像增强、地震属性分析还是其他方向性结构凸显任务上都是非常实用的工具箱。