1. 项目概述与核心思路

搞风资源评估这行的都知道一句话：场址选对了，项目就成了一半。而场址判断的地基，就是气象塔测风数据的质量。不管是用再好的风速计、再精密的采集器，数据进了电脑之后这一关要是处理不对，后面算出的发电量、等效满发小时数可能全都会跟着跑偏。我这次要讲的，就是我们在一处拟建风电场地，用Matlab把一整年气象塔的历史测风数据从导入、清洗到分项评估的完整过程，代码自己也跟着整理了，可以直接拿去改改参数复用。

先说说这套数据本身。用的是某测风塔70米/50米/30米三个高度层的10分钟平均序列，采集频率是10Hz原始采样后自动归算成10分钟均值，字段包括平均风速、最大风速、最小风速、标准差、平均风向、环境温度、气压、湿度，还有数据采集器的有效记录状态位。这种数据格式在国内的测风项目中非常主流，很多单位用某常规气象仪器厂商的采集器，输出就是这个样子。但很遗憾，绝大多数同事拿到手之后，第一件事就是把文件拖进Excel里直接做透视表。这个路子不是不行，只是做一年数据的资源评估基本上会被Excel卡到怀疑人生，而且后续要做威布尔拟合、湍流强度分扇区统计的时候，Excel的效率是真的拉胯。

我这次用Matlab做，倒不是因为Matlab在数据处理上比Python有多少绝对优势，主要是团队的老代码和环境都在Matlab上，而且Matlab的统计工具箱调用威布尔拟合、绘图的美观度都比较顺手。如果你更熟悉Python，其实思路一样，数据清洗和统计逻辑完全可以平移过去。下面我就把从零到结果的全过程拆开讲清楚，包括每一步为什么要这么做、参数怎么定、哪些地方一旦处理不当整个评估就会失真。尽量让刚入行的人看完就能照葫芦画瓢，也让有经验的人能对比一下自己的处理流程有没有遗漏哪一环节。

2. 数据导入与格式预处理

2.1 原始数据的常见形态与检查清单

不管你用什么工具做分析，第一步永远不是急着读数据，而是先搞清楚你手里这份数据到底长什么样。我处理的这份数据是每10分钟一行的ASCII文本格式，固定宽度字段，文件按月份切割成12个txt，每个文件大约4万行左右。如果你拿到的是CSV或者Excel格式，也问题不大，后面读入方式换一下而已。

在真正动手写导入代码之前，有四个检查项我必须先做，每次都是如此：

• 检查文件头和数据说明，确认每个字段的物理单位，风速是m/s还是节，温度是摄氏度还是开尔文，这些一旦搞错，后面的功率密度计算会差出几倍。
• 检查时间戳覆盖范围和数据点数量是否对得上。一年365天按10分钟一条记录算，标准应该是365×24×6=52560条。如果差得很多，说明采集器可能停过机，这个要心里有数。
• 检查有没有重复的时间戳。有些采集器在重启之后会把前面的缓存重新写入，导致同一时刻出现两条记录。
• 确认数据状态位的含义。很多采集器会用一个整型字段表示记录是否有效，比如0表示正常、1表示传感器异常、8表示校准中。这个字段是你做数据清洗的重要依据。

我把这些检查项做成了Matlab脚本里的前置断言，不满足就直接报错并退出。这个习惯帮我在处理几百份测风数据时避免过很多低级错误。比如每次采集器的时钟偏差就是个容易被忽略的问题，10分钟数据如果每天快30秒，一个月就能积累出15分钟以上的偏差。这种系统性误差会在逐日平均风速序列上露出锯齿状痕迹。

2.2 Matlab快速导入：从文本到表格

Matlab读这种固定宽度文本，我通常用readtable配合自定义的Format参数，或者更省事的textscan。这里我倾向于用textscan，因为测风数据经常有缺失值表现成-999或者9999这类填充值，用textscan指定空值策略会更直接。

一个常见的示例代码如下：

fid = fopen('met_tower_2021_07.txt', 'r'); % 按空格分隔，每行13列 C = textscan(fid, '%s %f %f %f %f %f %f %f %f %f %f %f %f', ... 'HeaderLines', 1, ... 'MultipleDelimsAsOne', true, ... 'TreatAsEmpty', {'-999', '9999', '-9999'}); fclose(fid); % 组合成timetable t = datetime(C{1}, 'InputFormat', 'yyyy-MM-dd HH:mm:ss'); T = timetable(t, C{2}, C{3}, C{4}, C{5}, C{6}, C{7}, ... C{8}, C{9}, C{10}, C{11}, C{12}, C{13});

这里有两个细节值得说。第一，TreatAsEmpty可以把填充值直接转成标准缺失值，避免后期到处找-9999这种魔法数。第二，用timetable而不用普通表格，是因为后面做重采样、按月份聚合、滚动统计时，timetable有非常顺滑的原生支持，尤其在Matlab新版本里对时间序列的处理函数比老版本强了太多。

在处理完所有单月文件后，我会把12个月的timetable直接纵向拼接起来，用vertcat，然后做一次去重。去重的逻辑是：先按时间戳排序，再检查相邻时间差是否小于10分钟，如果小于，就说明存在重复记录，保留状态位更可靠的那一条。

2.3 一次完整的读入效果检查

数据读入之后不能急着算指标，我先做一次快速的可视化检查。画三个东西：全年的平均风速时间序列、风向出现频次、以及温度气压的逐日变化趋势。不要小看这一步骤，它能在几分钟内暴露大规模的数据异常。

我在这份数据上就发现了一个问题：6月中旬的逐日平均风速整体明显偏低，低到比相邻月份的均值低50%以上。后来查了采集器的日志，才知道那段时间风速计被鸟粪严重遮挡，传感器输出的风速值大量低于真实值。这种程度的数据必须做标定替换或者整段剔除，绝不能留在样本里参与后面威布尔拟合。如果不做这步检查，而直接进入计算，测风塔数据的均值会被人为拉低，资源评估结果就会偏保守，导致项目决策失误。

3. 数据质量控制与清洗

3.1 三个层级的质控策略

原始测风数据不能直接用，这是行业共识。我在项目里遵循的是行业通行做法：一级质控查合理性，二级质控查相关性，三级质控查变化趋势。简单说，第一层检查物理量是否在合理范围，比如风速不能为负，极大风速不能超过某个阈值（比如60m/s，一般陆地测风塔不会出现这个量级），风向在0-360度之间，气压在合理海拔范围内。第二层利用同一塔上不同高度间的相关性，比如70米和50米风速的相关性通常非常高，如果某个时段两层的相关性突然跌破0.9，就要怀疑其中某一层的传感器出问题了。第三层是检查变化趋势，因为真实风速有强自相关特征，相邻10分钟的风速差通常不会超过一个经验阈值，如果连续多条记录的瞬时跳变过大，大概率是传感器结冰或雷击后的假数据。

这个三层体系听起来简单，实际操作中却经常被跳步。很多人只做第一层，把负值和超范围值剔掉就交差，结果数据里的结冰问题、传感器卡涩问题、塔影效应问题完全没有被识别出来，到了算湍流强度的时候，得出的数值高得离谱。

3.2 缺失值和无效值的处理办法

数据清洗的大头在于缺失值和无效值的填补取舍。首先是填充值的统一识别。我在代码里定义了一个统一的规则：

speedValidIdx = (speed > 0 & speed < 50) & (statusFlag == 0) ... & ~isnan(speed);

这个逻辑并不复杂，但要把用户自定义的边界值暴露成参数。不同场址有不同的极端气候特征，比如沿海台风区的风速上限就应该放宽到60甚至70，而内陆山地可能40以上就极可疑了。

处理缺失的策略我比较保守：短于3小时的连续缺失，用相邻时段的线性插值补上；超过3小时的连续缺失，整段标为无效，后续计算有效数据完整率的时候如实统计，不做任何填补。为什么用3小时这个阈值？因为3小时内的插值误差对均值的影响通常小于1%，而更长的连续缺失，受天气系统演变影响，插值结果很容易失真。如果你做的是瞬时极端风分析，比如计算极大风速重现期，那插值更是大忌，直接跳过缺失时段，用极值序列法处理。

3.3 塔影效应与扇区剔除的具体操作

测风塔本身的塔体、横臂和支架会对风场产生干扰，导致某一方向范围内的风速测量值出现系统性偏差。这是气象塔测风数据在做资源评估时绕不开的一个问题。常见的影响区域就是主导风向上的塔影区，通常根据风向传感器和风速计的相对方位，剔除以塔影方向为中心的60度左右扇区。

我这边的做法是利用对比测风设备的数据来判断塔影扇区。70米层有两个风速计，一个装在塔身的一侧，另一个装在相对的横臂上，它们各自会暴露在不同方位的影区中。统计两个风速计在每个风向扇区的平均风速差值，差值持续大于某个阈值的方位就是要被剔除的扇区。剔除之后，在这些扇区内不能直接留空处理，而是用另一支传感器在相同时间段的数据做替换，只有两者同时缺失时才标记为无效。

这种扇区剔除看起来是在损失数据，实际上是在保住整组数据的可信度。塔影区的数据如果不处理，各个扇区的平均风速就会出现一个虚假的“凹陷”或“凸起”，直接影响风向玫瑰图中的能量分布，后续的风机布局设计就会偏向一个并不存在的优势方位。

3.4 结冰和雷击损坏的数据识别

对于冷地区或者雷暴多发地区的测风塔，传感器结冰和雷击损坏是最让人头疼的两类异常。结冰时风速计叶片被冻住，转速几乎归零，或者随风摆动产生一系列不规则的低速峰值。这些数据在合理性检查里完全不会暴露，因为它们的数值都在合法范围内，但它们的持续时间和出现时段往往是模式化的，比如多出现在夜间湿度大温度低且风速在3-8m/s的区间。

我的识别方法是用“标准差异常”配合“持续低速”两个条件。具体来说，当温度低于2摄氏度、湿度高于85%、风速低于3m/s且持续超过2小时，就标记为疑似结冰区间。此外我也会对比同一塔上不同高度的风速序列，如果低层开始结冰而高层风速正常，就能判断是结冰导致的局部传感器失效。对于这部分时段，我会整体剔除，而不是用插值填补，因为结冰期间的记录已经失真，插值等于在人为制造虚假的风速变化曲线。

4. 核心评估指标与Matlab实现

4.1 平均风速与风功率密度

平均风速是最基础的指标，但算起来有几个细节需要注意。比如“按算术平均”和“按时间加权平均”的区别。10分钟平均风速序列本身就是等间隔的，直接算术平均没问题。但如果数据因为缺失导致时间间隔不均匀，就必须用时间权重，否则缺失多的时段权重会被低估。我一般用retime方法统一到标准的10分钟网格，再做平均。

风功率密度表达式是：

[ WPD = \frac{1}{2} \rho V^3 ]

这个式子里空气密度ρ不能拍脑袋取常数。在同样的风速下，ρ差了10%，功率密度就能差10%，这比很多测量误差还要大。我会用理想气体方程根据实测气压和温度计算：

[ \rho = \frac{P}{R_d \cdot T} ]

其中(P)是环境气压（Pa），(R_d)是干空气比气体常数，约等于287.05 J/(kg·K)，(T)是开尔文温度。这样每天、每十分钟都有对应的空气密度值，计算风功率密度时用的是逐条的密度而不是全年的平均值。我在这份数据的计算中发现，冬季冷高气压期间的空气密度比夏季高出约12%左右，这意味着同一台风机的冬季出力天然会高于夏季，全年发电量评估如果忽略了这种季节密度差异，会带来数个百分点的不确定性。

计算代码核心部分如下：

rho = pressure ./ (287.05 * (tempC + 273.15)); WPD = 0.5 * rho .* (speed .^ 3); annualAverageWPD = mean(WPD, 'omitnan');

在实际项目中，我还会把逐月的功率密度算出来，观察季节变化，同时也方便和周边测站的历史数据做趋势对比。一年数据太短，不足以完全代表长期气候态，但趋势本身最能说明这个场址的气候特征是否与附近长期气象站一致。

4.2 威布尔分布拟合：极大似然估计与最小二乘对比

风频分布是风力发电评估里的核心要素，而威布尔分布是最通用的描述风速概率分布的模型。两参数威布尔分布的概率密度函数为：

[ f(v) = \frac{k}{c}\left(\frac{v}{c}\right)^{k-1}\exp\left(-\left(\frac{v}{c}\right)^k\right) ]

其中(k)是形状参数，(c)是尺度参数。(k)值决定了风频曲线的形态：(k)小，风速变化分散；(k)大，风速集中在特定区间。一般来说(k)在1.5到3之间都算正常，我这个场址的70米高度拟合出的(k)大约为2.1，(c)大约为7.8m/s。

我推荐用极大似然估计来做拟合，而不是把风速分成区间后做最小二乘。最大似然的优点是充分利用了原始样条数据的信息，不会因为分箱宽度的选择导致结果变化。分箱最小二乘如果箱子宽度取不好，拟合出的威布尔参数波动很大，尤其是10m/s以上的高风速区间，箱子里样本太少，拟合就非常不稳定。

Matlab里一行代码可以直接调：

pd = fitdist(speedValid, 'wbl'); k = pd.a; % 尺度参数c c = pd.b; % 形状参数k，注意Matlab这里的命名顺序和常用符号不一致

这里有一个特别容易踩的坑，Matlab的fitdist威布尔分布的参数顺序和统计学教材里的符号有出入。教材上一般写作(c)(尺度)和(k)(形状)，但Matlab里返回的pd.a其实是尺度参数，pd.b才是形状参数。以前我写报告时就搞反过一次，输出的k和c对调了，导致估算的年发电量偏了将近5%。后来我在所有代码里明确规定：cScale = pd.a; kShape = pd.b;并加注释，彻底杜绝此类问题。

4.3 风玫瑰图与主导风向分析

风玫瑰图是风资源报告里必不可少的一张图。它既呈现风向频次的分布，也叠加各扇区的平均风速，以此识别主导风向和功率密度集中的方向。作图本身在Matlab里并不困难，难点在于如何选择风向扇区的数量。我比较常用的是16扇区，每个扇区22.5度。用更少的扇区会掩盖比较窄的盛行风通道，而使用更多扇区又会因为每个扇区内样本量太少而出现随机起伏。

这里我给出一段可以直接跑的绘制代码框架：

dirEdges = 0:22.5:360; dirCenter = dirEdges(1:end-1) + 11.25; [counts, ~] = histcounts(dir, dirEdges); freq = counts / sum(counts) * 100; % 计算每个扇区的平均风速 meanSpeedBySector = zeros(16, 1); for i = 1:16 idx = dir >= dirEdges(i) & dir < dirEdges(i+1); meanSpeedBySector(i) = mean(speed(idx), 'omitnan'); end

在绘图时，我习惯把风向频率的平均风速同时画在一个极坐标下，用条形高度表示频率，用颜色或第二层折线表示风速。这样一张图就能看出来哪个方向风出现得最多，哪个方向风来得最猛烈。有时候频率最高的方向和功率密度最高的方向并不是同一个，这个时候风机排布就要优先考虑功率密度高的方向，而不只是风向频率。

4.4 湍流强度与风切变指数

湍流强度定义为风速标准差与平均风速的比值，是衡量风况稳定性的重要参数，直接关系到风机的疲劳载荷等级和湍流等级选择。行业里有明确的湍流等级划分标准，如果测得的湍流强度偏高，选型时就得选用高湍流等级的风机或者适当降低轮毂高度的预期发电量。我的计算方法是先按10分钟数据逐条计算(I = \sigma / V)，然后按风速段和风向扇区分别做平均。这里论证一下，直接对所有数据求一个总体平均湍流强度，会背离湍流强度随平均风速的递减特性，因为低风速时段湍流强度天然偏高。真正有意义的是在风机运转风速区间（比如3-25m/s）内，按风速段统计湍流强度的平均值和90分位值。

风切变指数则反映了风速随高度的变化规律。经典的幂指数公式：

[ \frac{V(z_2)}{V(z_1)} = \left(\frac{z_2}{z_1}\right)^{\alpha} ]

由此可以推算出(\alpha)：

[ \alpha = \frac{\ln(V_2/V_1)}{\ln(z_2/z_1)} ]

我使用70米和50米两个高度的数据计算每小时的平均风切变指数。这个参数对风机轮毂高度的选择有直接意义。如果风切变指数偏大，说明风资源随高度的增加很明显，那么加高塔筒的收益就非常可观。我这边的数据算下来，全年平均(\alpha)约为0.18，这个值比教科书上常用的0.2略低，但在当地复杂地形背景下很合理。而且我还注意到冬季夜间(\alpha)明显升高，这与夜间地表逆温层的形成有关，对风机叶片结冰的风险评估也有参考价值。

5. 实操过程与关键环节代码解析

5.1 完整流程逻辑与模块划分

为了便于复用，我把整个处理流程编写成了模块化的Matlab脚本集。总共有五个模块：

1. 读入与初查模块：负责读入txt、拼接、时间校对、重复检查。
2. 质控清洗模块：负责任何标记异常、扇区剔除、插值和无效时段标注。
3. 统计计算模块：计算逐条WPD、威布尔参数、风切变指数、湍流强度、分扇区统计表。
4. 绘图输出模块：生成平均风速时间序列、风玫瑰图、威布尔拟合对比图、风速日变化/月变化图。
5. 报告表格模块：以表格形式输出逐月统计表、分扇区统计表、有效数据完整率统计表。

模块化带来的最大好处是：换一个场址的数据，只需要修改路径、塔高、扇区参数和数据状态位的配置项，不需要改动计算逻辑。我在好几个项目里都是这样干的，省掉了大量重复劳动，也不用担心因为改代码导致某个环节悄悄出错。

5.2 逐月与整体统计的参数选择

在做逐月统计的时候，要特别注意有效数据完整率在不同月份的差异。测风塔全年完整率在95%以上是比较理想的。如果某个月完整率掉到90%以下，那这一个月的统计结果就不应该参与整体的年平均值计算，而要做相应的合理性说明。我的做法是，先输出一张逐月完整率表，凡是完整率低于90%的月份，在最终年均风速计算中按无效处理，并在报告里注明。

另一个关键参数是风速区间的分段。我按1m/s为一档做风速频率分布表，对低风速区（0-3m/s）单独归拢，对高风速区25m/s以上合并为一档。这样做一方面是为了威布尔图的展示清晰，另一方面是因为高风速端的样本量天然较少，过分细化没有统计意义。

5.3 绘图细节与输出格式

Matlab绘图有几个小技巧我觉得很值得分享。比如做时间序列长图时，如果画出全年的逐10分钟曲线，整张图就是一片密密麻麻的黑色，什么都看不清楚。我一般会先画逐时平均的曲线，然后叠加逐日平均的粗线，这样既能看出波动细节，又不会失去宏观趋势。

威布尔拟合图我习惯画成风速累积概率和概率密度两个子图。左边是经验累积频率与拟合累积分布对比，右边是概率密度直方图与拟合曲线对比。看到两条曲线贴合得好，拟合质量才有说服力。除了图以外，我还会用kstest做拟合优度检验，不过要提醒一句：当样本量达到几万条时，KS检验几乎总会拒绝原假设，这是因为真实风速分布和理论威布尔分布之间总有细微偏差，而大样本把这种偏差放大到“显著”。所以检验结果仅供参考，不宜作为唯一标准。

6. 常见问题与排查技巧实录

6.1 测风数据处理常见问题速查表

能把下面这张表里的问题都事先想一遍，处理测风数据时能少走很多弯路。

6.2 几个实操之后才明白的细节

第一件事是关于时间格式的解析。Matlab的datetime在读取日期字符串时偶尔会因为系统时区设置而产生偏移。如果你用的是UTC时间而本机时区设置成了东八区，那datetime缺省会直接按本地时区解析，导致所有时间戳偏移8小时。这种错误在后续的日变化分析中会表现为风速峰值相位偏移。我现在的做法是在读入后立刻检查时间戳的首末记录是否和文件名中的月份范围一致，不一致就停下来排查。

第二件事是不要迷信“所有风向扇区都需要剔除塔影”。不同塔型、不同横臂布置，塔影区范围并不相同。我见过有人直接把整个正北方向60度范围内的数据全部删掉，理由是多数报告都这么写，但实际上那支塔的风速计装在朝北一侧，塔影区应该在正南。这种照搬模板的做法很危险，直接削减了有效数据量，还引入偏差。必须根据本塔的实际传感器几何布局来判断。

第三件事是关于数据备份的版本管理。我习惯在每跑完一个清洗步骤后，先把结果保存成中间文件，而不是覆盖原文件。这样做的好处是，如果有人质疑最终结果，我可以回溯到每一步，快速定位是哪一步的筛除范围偏大或者插值影响了统计。有一回我就是靠中间文件发现扇区剔除参数里有一处角度边界写错了，导致尾流方向的几个扇区被错误剔除。要不是留有中间结果，这个错误可能要等到报告写完才能发现。

还有一个小而关键的细节：计算空气密度的时候，气压计海拔修正。测风塔上的气压传感器如果安装在塔基处，读数对应的海拔与轮毂高度处的气压有差别。虽然气压随高度的衰减对密度的影响并不大，不到明显误差级别，但既然我们是做评估，严谨一点总没坏处。我的处理是用测站海拔统一到一个参考面，保证同一项目里不同测风塔的数据具有可比性。

7. 这份评估的后续延展

数据清洗完毕，威布尔参数和风功率密度都算出来了，并不代表风资源评估工作就结束了。气象塔数据只是点上的测量，真正要给整个场址下一个结论，还要结合中尺度模拟降尺度数据、地图地形分析、尾流模型测算以及周边已运行风电场的实际发电数据做交叉验证。如果是复杂山地地形，一个测风塔的数据远不足以代表整个场址的各个机位，就需要补充临时测风设备或者使用相关分析推算法进行空间推算。

我个人实际操作下来的体会是，数据处理过程的确定性其实很强，工具用什么、代码怎么写，都不是最难的关口。真正考验经验的，是决定哪些数据可用、哪些不可用，哪段插值合理、哪段必须整段作废。这些判断没有统一的算法公式，靠的是对整个测量系统误差来源的深刻理解和对数据曲线特征的敏锐直觉。希望这次完整的处理流程拆解，能让刚开始接触测风塔数据的朋友少走些弯路。最后再分享一个小技巧：把处理好的timetable和最终的统计结果在项目结束后原样保留一份，连同关键参数配置文件一起打包归档。以后如果项目要扩容或者做后评估，这套数据能直接复用，省下的时间远比你当时整理它所花的时间多得多。