MATLAB通量足迹计算与源区分析：从湍流物理到FFP模型实现

简介：涡动通量来源区分析模型是一套基于MATLAB开发的通量源区贡献与印痕分析工具，面向生态学、气候学、农业气象等领域的研究人员，帮助量化大气与地表之间能量、物质（如水汽、二氧化碳）交换的来源区域。压缩包共12个文件，主要包含5个.m脚本、3个PDF文档、2个txt说明及图片等辅助材料，整体仅5.51MB，轻量易部署。其中calc_footprint_FFP.m用于根据风向、风速等气象参数逐时刻计算通量印痕，calc_footprint_FFP_climatology.m处理长时间序列数据，生成通量源区的平均、最大与最小贡献区域等气候学特征；PDF文档详细说明了模型原理、输入参数与运行配置，许可文件则明确了使用规范，另有辅助文件夹提供扩展支持。目前已有216人学习下载，适合需要开展地表通量空间来源分析的研究生、科研人员或环境评估从业者。借助这套模型，可深入揭示农田蒸散模式、森林碳吸收季节变化及城市热岛效应等场景下通量的空间来源与动态过程，为碳平衡评估、气候变化影响预测和环保策略制定提供基础数据支撑。

1. 通量源区：你测到的信号，到底是哪一片地表的“贡献”？

涡动相关系统架在观测塔上，很多人默认它测的是塔基附近那一小片土地的碳通量或水汽通量。真实的答案是：传感器读数是上游一整片地表信号的加权平均，这个权重分布就是通量足迹（footprint），也叫源区贡献或印痕。它可能近到几十米，也可能拉长到几公里，随风速、大气稳定度和边界层高度剧烈变化。如果只拿通量数值说话而不看足迹，等于看一张没有比例尺的地图。更危险的是，农田通量可能被算进隔壁林地的贡献，生态系统的碳收支结论会从一开始就站歪。这就是为什么涡动通量来源区分析模型值得拆开来看：它把“数据从哪来”变成一个可量化的二维印痕，而 MATLAB 实现恰好把底层物理从粒子追踪降维成了可复现的函数调用。接下来的每一章，都围绕这套印痕分析的物理机制、解析近似和源码结构展开。

2. 通量足迹的物理根基：为什么传感器收到的是上游整条“记忆”

2.1 湍流输运：地表信号到达传感器的搬运链

涡动相关系统测的是垂直风速 w 与标量浓度 c 的协方差，得到一段时间内的平均垂直通量。这件事成立的前提是通量源区假设：空气块在抵达传感器之前贴着地表运动，沿途不断接收或释放标量，湍流涡旋则把这些携带信号的气团逐级混合、向上输送。于是传感器所在的固定点，实际上是一系列上游位置信号的“统计汇”。

我一般会这么理解：空气流过下垫面时，像一支接力队。离塔很远的涡旋，把接近地表的标量交给稍下游的涡旋，再下游的涡旋又混合一次。每传一次，原始信号被稀释、被重新混合，等到了传感器位置，早先位置的贡献已经衰减到几乎为零。这种湍流传递不是教科书里那种局地的分子扩散，而是准随机、大尺度、强间歇的，所以“来源”天然是一片区域，而不是一个点。这也是足迹模型跟简单的“上风向扇形范围”类经验公式最本质的差别——它描述的是连续概率密度，不是一个粗糙的半径。

这个过程的数学形式是面积分：实测通量 F 等于足迹函数 f(x, y) 与地表源强 Q(x, y) 的卷积。f(x, y) 是权重，Q(x, y) 是每块地表释放或吸收标量的速率。如果测到 F 之后想反推 Q，就必须先知道 f。这个反演问题本身是病态的，但在大气边界层相对均匀的前提下，f 可以通过流场统计量很好地近似出来。正是这种“用统计量替代逐点粒子追踪”的思路，让 FFP 这类解析模型成为可能。

再往细里说，足迹的形成与湍流的拉格朗日时间尺度有关。粒子在湍流场中保持速度记忆的时间越久，它在水平方向飘得越远，足迹就拉得越长。当边界层高度有限时，粒子碰到边界层顶会被反射回来，这相当于在垂直方向给湍流加了顶盖，改变了粒子落地的空间分布。这些物理过程全部需要显式出现在模型里，而不是简单套一个高斯烟羽公式就能蒙混过去。

2.2 主控变量：稳定度、粗糙度与边界层高度如何重塑足迹

对大多数实测条件，足迹的形状由四个量控制：测量高度 zm、摩擦速度 u*、奥布霍夫长度 L、边界层高度 h，再加上横向湍流强度 σv 和平均风速 U。下图（文字描述）是我刚接触这套模型时建立的对应关系，省了不少弯路：

这里最反直觉的是稳定层结。夜间晴空、风速很小时，湍流几乎被抑制，地面释放的标量只能在浅层缓慢输送，传感器测到的可能是几公里外的源。白天强对流时，湍流把信号快速混合，足迹往往压缩到百米量级。所以“足迹到底有多大”不是模型设定的固定值，而是由实时湍流状态决定的动态变量。这也是为什么 calc_footprint_FFP 要求输入的是半小时平均的湍流统计量，而不是随便给一个粗糙度参数就能出结果。

粗糙度长度 z0 和零平面位移 d 在这里也起关键作用。z0 决定地表对风的拖曳，间接影响 u*；d 则是高植被（森林、农田）下有效测量高度的修正量。城市和森林场景里如果忽略 d，会把有效测量高度算高好几米，峰值距离可能整体偏移几十到上百米。很多用户把 z0 当成“查表填一下”的参数，但实际上它对足迹的伸缩作用在近地层稳定条件下非常敏感。

2.3 足迹的数学语义：一个归一化的贡献权重场

足迹函数 f(x, y) 的单位是 1/m²，含义很直接：空间中每一块地表微元对观测通量的贡献权重。对整个平面做面积分，结果等于 1，因为传感器测到的标量总量必须全部来自上游那些位置。f 的峰值通常出现在迎风方向某个距离处，这个距离大致正比于测量高度与平均水平风速的比值，再除以一个湍流速度尺度。

想通了这个归一化约束，再去看程序输出就会容易很多。f_2D 矩阵就是离散化的权重场，矩阵里每个元素乘以网格面积，才是该网格的贡献占比；如果直接把矩阵元素相加，在不同网格密度下会得到完全不同的“贡献面积”。这也是后面避坑章里要展开的常见误用之一。理解 f(x, y) 是密度函数而不是通量本身，是读对结果的第一道门槛。

3. 从粒子漫游到解析公式：2015 年新足迹方案的“降维”思路

3.1 拉格朗日随机模型的算力账单

拉格朗日随机模型（LPDM 一类）的思路非常直接：从传感器位置反向释放大量粒子，让粒子在湍流风场中做随机游走，最后统计粒子落回地面的位置，落地密度分布就是足迹。这个方法物理上最可信，能处理非均匀地形、复杂边界层，但代价高得吓人。

一个典型的模拟要释放几万个粒子，每个粒子在积分时间步里都要根据局地湍流统计量随机更新速度，还要处理地面反射、边界层顶反射。风速、稳定度一变，整套模拟就得重跑一遍。这意味着你无法用它对几年、几十个站点的通量数据进行批处理，更不能作为常规数据处理流程的一部分。2015 年的新足迹方案（也就是 FFP 的底层模型）正是在这个背景下出现的：先承认拉格朗日模型是“正确答案”，再用大量模拟结果去拟合一个解析近似，把粒子追踪的开销压缩成一次函数求值。

3.2 解析近似的灵魂：无量纲化与统计矩回归

这个方案的核心可以拆成两步。第一步，在大量不同的边界层条件下跑拉格朗日模型，把得到的足迹函数全部放在同一个“无量纲坐标系”里比较。坐标轴的缩放因子不是随便选的，而是跟湍流统计量绑定：水平距离用测量高度和 σw 的组合来归一化，垂直方向的混合强度用边界层高度和奥布霍夫长度来刻画。经过无量纲化之后，不同气象条件下形状迥异的足迹曲线，会塌缩到一组相似的“标准曲线”上。

第二步，对这组标准曲线做统计矩分析。每条足迹曲线都有峰值位置（一阶矩）和水平展宽（二阶矩），这两个量随稳定度参数、边界层高度、测量高度变化，可以被拟合成连续函数。拟合完成后，模型就不再需要粒子追踪了。你只要输入 u*、L、h、σw 这些实测统计量，就能立刻得到足迹的峰值位置和展宽，再由解析公式重建整个权重场。

calc_footprint_FFP.m 的内部结构大致就是这条链路的代码化，常见的实现如下：

function [x_2D, y_2D, f_2D] = calc_footprint_FFP(zm, z0, d, umean, h, ... sig_v, ustar, L, wind_dir, peak, sigma_xy) % 1. 由测量高度、零平面位移计算无量纲稳定度参数 zeta = (zm - d) / L; % 2. 根据 ustar、h、L 估计垂向湍流速度标准差 sig_w sig_w = estimate_sig_w(ustar, L, h, zm - d); % 3. 用回归公式计算无量纲峰值距离 x_star 和展宽 sigma_star [x_star, sigma_star] = footprint_moments(zeta, h, zm - d, sig_w); % 4. 按实际风速尺度和测量高度还原为有量纲距离，构造二维权重场 [x_2D, y_2D, f_2D] = build_field(x_star, sigma_star, ... umean, sig_w, wind_dir); end

这个结构里的关键是第 3 步：x_star 和 sigma_star 不是靠模拟粒子算出来的，而是从回归公式直接得到的。如果拿放大镜看 footprint_moments 函数，会发现里面全是针对不同稳定度区间分段定义的拟合系数，这正是“用离线模拟换在线速度”的代价——精度上限完全取决于拟合样本覆盖的参数范围。超出拟合范围的极端稳定层结，结果会明显失真。

3.3 单峰足迹与交叉风积分：二维矩阵的成形顺序

得到峰值位置和展宽之后，二维足迹矩阵不是直接算出来的，而是分两步“组装”的。先沿主风方向构造一维足迹剖面，这个剖面是典型的单峰曲线，峰值出现在迎风向某个距离处，峰值前后不对称，上风向衰减较慢、下风向衰减较快。然后再叠加上垂直于风向的扩散分布，通常假设为高斯型，扩散宽度由横向湍流强度 σv 和风速共同决定。

这两个方向的信息组合成一个完整的二维贡献场后，程序会再把坐标旋转到实际地理方向。这里的 wind_dir 参数必须使用气象学约定的“风吹来的方向”，而这个约定在代码注释里往往被一笔带过，实际用的时候特别容易翻车。二维矩阵的网格范围和分辨率也不是随意给定的：网格得足够细，否则峰值位置会被抹掉；网格范围得足够大，否则 90% 贡献区域会跑到边界外面。下一章会对输出语义做更细的拆解。

4. MATLAB 实现链路：气象输入到足迹矩阵的完整管线

4.1 输入参数清单与气象数据预处理

calc_footprint_FFP 的输入接口设计得很克制，只吃湍流统计量，不吃原始高频数据。也就是说，在调用它之前，数据已经完成了从 10Hz/20Hz 原始采样到半小时平均统计量的压缩。这也隐含着一个常见要求：先对原始通量数据做质量控制和稳态检验，再把合格时段的气象统计量喂给足迹计算。如果原始信号本身不稳，算出的足迹再精确也没有物理意义。

一个典型的调用是这些参数：

这些参数里最容易被忽略的是 h。很多人随手给一个固定值（比如 1000 米），在小尺度不稳定条件下影响不大，但遇到稳定边界层，h 直接决定足迹尾部能延伸多远。h 给大了，足迹末尾多出一大段虚假的远距离贡献；h 给小了，峰值可能被错误地推到塔下。

4.2 输出结构里到底装了什么：x_2D、y_2D、f_2D 的语义

函数返回的三个输出 x_2D、y_2D、f_2D，前两个是网格坐标矩阵，第三个是对应的足迹权重密度。这三个矩阵同一维度，描述了一个以观测点为中心、固定旋转角度的二维网格。风的方向决定 x 轴正方向，y 轴与之垂直。单独看这三个变量，无法直接判断“源区有多大”，必须结合网格间距一起理解。

取 90% 贡献面积是常用的下游操作，典型的实现是这样：

% 从输出矩阵提取 90% 足迹等值线内的源区面积 dx = abs(x_2D(1,2) - x_2D(1,1)); dy = abs(y_2D(2,1) - y_2D(1,1)); cell_area = dx * dy; % 按贡献权重降序排列，累计求和到 0.9 f_sorted = sort(f_2D(:), 'descend'); csum = cumsum(f_sorted) * cell_area; idx90 = find(csum >= 0.9, 1, 'first'); threshold = f_sorted(idx90); % 所有权重大于阈值的网格面积之和，即 90% 源区 area90 = sum(f_2D(f_2D >= threshold)) * cell_area;

这段代码里最容易出错的是降序累计这一步。因为足迹权重分布高度偏斜，峰值附近的网格贡献极其集中，如果不用降序累计，而是直接找“f > 0.1 倍峰值”的网格，算出的面积会远大于真实贡献区域。降序累计保证了 90% 是真正贡献的 90%，而不是网格数的 90%。同样的逻辑可以换成 50%、80%，取决于你的分析目的是判断主源区还是完整影响半径。

另一个输出细节：f_2D 的元素值通常很小，尤其在网格间距较大的时候。这是因为足迹函数是密度函数，约束条件是面积分等于 1，而不是元素求和等于 1。如果想要的是“各个方向各占几成”，建议先对 f_2D 做按方向的积分，而不是直接对矩阵做行列求和。

4.3 气候学版本：长时间序列怎么聚合成一张代表性源区图

calc_footprint_FFP_climatology 是这套脚本里第二个核心入口，对应长期数据的需求。单次足迹只能回答“这半小时信号来自哪里”，但做碳收支评估时，需要回答的是“这个站点一整年、或者某个季节，信号主要来自哪里”。这是一种统计性的源区描述，不能靠把每半小时的足迹矩阵简单相加来得到。

常见的实现思路是：先按季节或月份分组，把组内的风向、风速、稳定度条件做分箱统计，选出具有代表性的气象条件组合，再对每个组合分别计算足迹，最后按各组合在时间序列中的出现频率做加权叠加。这样得到的平均足迹矩阵，代表的不是一个特定时刻，而是一个气候平均意义下的源区概率。它的形状通常比任何单次足迹都平滑，多个风向的贡献叠在一起会呈现扇形或双峰形态。

从代码结构上看，climatology 版本底层会反复调用瞬时版函数，但外层多了一层统计聚合逻辑。如果只读主函数，很容易以为它只是循环调用的“批处理壳子”，但实际上分箱、加权、平滑这些步骤决定了气候学结果跟瞬时结果在统计上是否自洽。后面避坑章里有一条专门讲这两者为什么对不上。

5. 通量足迹计算避坑：最容易翻车的四种现场

5.1 夜间稳定层结下的大范围“假尾巴”

现象：白天结果正常，夜间单独跑某一时段，足迹矩阵显示贡献范围延伸到几十公里外，明显超出站址周围下垫面类型，完全没法跟土地利用图对照。

原因：稳定层结下湍流极弱，足迹公式外推到拟合数据覆盖范围以外的区域。2015 年方案的解析近似主要基于对流边界层和中等稳定条件下的拉格朗日模拟，强稳定、低风速场景属于外推区，模型输出会夸大远端贡献。

解决：先按 L 或无量纲稳定度参数做筛选，只在模型有效范围内计算足迹。常见做法是设定稳定度阈值——当 (zm-d)/L 大于某个正值时，直接放弃该时段足迹分析，改用其他方式描述源区。阈值写在数据处理脚本开头，比跑完再检查矩阵节省大量时间。

5.2 足迹峰值出现在塔的下风方向

现象：某次计算结果的峰值位置坐标和预期完全相反，看图像是“信号从塔后吹来”，气象上说不通。

原因：wind_dir 参数的风向约定搞反了。气象风向北往报告风向是“风的来向”，但有些观测系统记录的是“风的去向”，或者代码内部旋转坐标系时使用了不同的角度基准。峰值跑到塔后，通常意味着旋转方向反了 180 度。

解决：检查风向变量的单位制和基准方向。把输入 wind_dir 与同时段的风向计原始记录对照，确认是 0 度表示北风还是南风。统一后重新计算，峰值位置就会回到迎风方向。这个坑很难靠读代码发现，最好在测试阶段就用一组已知条件的样例数据验证输出。

5.3 气候学平均足迹和逐时足迹叠加结果对不上

现象：用 calc_footprint_FFP_climatology 算出的平均源区，跟手动把几百个小时的瞬时足迹做平均的结果差别很大，怀疑是程序有 bug。

原因：气候学版本不是简单平均，而是对气象条件作了分箱、筛选和加权处理，它可能排除了不满足稳态条件、或超出模型适用范围的时段。也就是说，两个结果的对象根本不是同一个样本集合，数值对不上是正常的。

解决：先确认你的“手动平均”是否做了同样的样本筛选。如果只是把有效时段足迹矩阵逐像元求平均，记得按每个时段的贡献权重归一化，而不是简单算术平均。做碳收支分析时，建议以气候学版本的结果为准，因为它天然规避了极端气象条件的干扰。

5.4 足迹矩阵的网格边界截断了部分贡献区

现象：90% 贡献面积算出来比预期小，或者某些风向下的足迹延伸区域在矩阵边缘处被切掉，等值线不完整。

原因：输出网格范围设置太小。默认网格是为常见气象条件设计的，高塔、强稳定或大风条件下，足迹特征距离会超出默认网格边界，程序不会自动扩展网格。

解决：看 f_2D 矩阵最外圈元素值，如果仍有显著非零权重，就要增大网格范围重新计算。另一个快速验证方法：检查 x_2D 的最大值是否大于该时段峰值距离的 5 到 10 倍，不够就扩展。网格越大、分辨率越高，计算时间线性增加，这一参数值得在批处理脚本里做成自适应。

5.5 输入了未经质量控制的湍流统计量

现象：足迹计算结果数值跳动剧烈，相邻十分钟算出的峰值位置相差几十倍，完全不像连续变化。

原因：涡动协方差输出的 u*、L 在小风速或降雨时段会有异常值，比如 u* 接近零、L 符号频繁翻转，这些异常值进入回归公式后被放大成极大的足迹尺度。

解决：在调用足迹函数之前，先加一轮通量数据质量控制。常见做法包括：u* 低于阈值（比如 0.1 m/s）的时段剔除，L 绝对值的合理范围做限制，缺失或超出物理范围的 L 用相邻时段插值。这段逻辑加在数据预处理里，不算模型本身的事，但能挡住九成以上的离群结果。

6. 解析近似的边界：什么时候该信，什么时候该撤

FFP 这类解析足迹模型的价值，本质上是拿精度换规模。拉格朗日随机模型准确，但面对几年、十几个站的通量数据，计算成本不可接受；解析近似牺牲掉复杂地形和非均匀层结下的精度，换来的是可以批量处理长序列的能力。理解这个取舍后，你就能给每个应用场景画一条边界：均匀平坦的下垫面、中等稳定度范围、对流边界层，是这套模型的主场；复杂山地、城市街区、强稳定夜间的数据，最好不要强行套用。

我一般会用一个很直接的办法验证模型在这个站点是否适用：把足迹矩阵跟土地利用栅格叠加，看看 80% 贡献面积内是否包含了明显不同的下垫面类型。如果 80% 贡献区里有超过两种地类，或者足迹边缘跨越了水体与森林的边界，就该警惕单点观测的代表性了。解析足迹算法本身声明了自己均匀地表的假设，这种跨越本身就是模型不适用的信号。

从那以后，我每次跑完足迹结果，都强制走一遍数据检查流程——先看峰值位置是否与风向一致，再看 90% 贡献面积是否落在下垫面栅格的合理范围内，最后抽查几个夜间稳定时段的输出有没有外推失真。这三步花不了几分钟，但能挡住绝大部分翻车现场。足迹模型不是黑匣子，它每一处输出的形状背后都有湍流物理的约束；读懂了这些约束，你就能放心地在论文里写“通量源区主要贡献来自……”，而不是心里打鼓。希望帮到你。

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