★ FEATURED ARTICLE 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×652560条。如果差得很多说明采集器可能停过机这个要心里有数。检查有没有重复的时间戳。有些采集器在重启之后会把前面的缓存重新写入导致同一时刻出现两条记录。确认数据状态位的含义。很多采集器会用一个整型字段表示记录是否有效比如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(i1); 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脚本集。总共有五个模块读入与初查模块负责读入txt、拼接、时间校对、重复检查。质控清洗模块负责任何标记异常、扇区剔除、插值和无效时段标注。统计计算模块计算逐条WPD、威布尔参数、风切变指数、湍流强度、分扇区统计表。绘图输出模块生成平均风速时间序列、风玫瑰图、威布尔拟合对比图、风速日变化/月变化图。报告表格模块以表格形式输出逐月统计表、分扇区统计表、有效数据完整率统计表。模块化带来的最大好处是换一个场址的数据只需要修改路径、塔高、扇区参数和数据状态位的配置项不需要改动计算逻辑。我在好几个项目里都是这样干的省掉了大量重复劳动也不用担心因为改代码导致某个环节悄悄出错。5.2 逐月与整体统计的参数选择在做逐月统计的时候要特别注意有效数据完整率在不同月份的差异。测风塔全年完整率在95%以上是比较理想的。如果某个月完整率掉到90%以下那这一个月的统计结果就不应该参与整体的年平均值计算而要做相应的合理性说明。我的做法是先输出一张逐月完整率表凡是完整率低于90%的月份在最终年均风速计算中按无效处理并在报告里注明。另一个关键参数是风速区间的分段。我按1m/s为一档做风速频率分布表对低风速区0-3m/s单独归拢对高风速区25m/s以上合并为一档。这样做一方面是为了威布尔图的展示清晰另一方面是因为高风速端的样本量天然较少过分细化没有统计意义。5.3 绘图细节与输出格式Matlab绘图有几个小技巧我觉得很值得分享。比如做时间序列长图时如果画出全年的逐10分钟曲线整张图就是一片密密麻麻的黑色什么都看不清楚。我一般会先画逐时平均的曲线然后叠加逐日平均的粗线这样既能看出波动细节又不会失去宏观趋势。威布尔拟合图我习惯画成风速累积概率和概率密度两个子图。左边是经验累积频率与拟合累积分布对比右边是概率密度直方图与拟合曲线对比。看到两条曲线贴合得好拟合质量才有说服力。除了图以外我还会用kstest做拟合优度检验不过要提醒一句当样本量达到几万条时KS检验几乎总会拒绝原假设这是因为真实风速分布和理论威布尔分布之间总有细微偏差而大样本把这种偏差放大到“显著”。所以检验结果仅供参考不宜作为唯一标准。6. 常见问题与排查技巧实录6.1 测风数据处理常见问题速查表能把下面这张表里的问题都事先想一遍处理测风数据时能少走很多弯路。问题现象可能原因解决方法某时间段风速持续为0风速计结冰、堵塞或信号断线检查状态位和温度湿度数据标记为无效连续记录出现固定值采集器传输故障或通道锁定检测序列中重复恒定值长度做剔除风向集中在某个角度风向标卡涩或冰冻对比风向频率分布检查累积频率奇异点风速昼夜差异过大可能正常也可能传感器受周边影响与相邻高度层做相关性分析确认气压突跳气压计受潮或雷击损坏检查维修记录剔除异常时段两个高度层风速几乎完全相同传感器高度配置写错或标定混乱用风切变公式反推高度差是否合理6.2 几个实操之后才明白的细节第一件事是关于时间格式的解析。Matlab的datetime在读取日期字符串时偶尔会因为系统时区设置而产生偏移。如果你用的是UTC时间而本机时区设置成了东八区那datetime缺省会直接按本地时区解析导致所有时间戳偏移8小时。这种错误在后续的日变化分析中会表现为风速峰值相位偏移。我现在的做法是在读入后立刻检查时间戳的首末记录是否和文件名中的月份范围一致不一致就停下来排查。第二件事是不要迷信“所有风向扇区都需要剔除塔影”。不同塔型、不同横臂布置塔影区范围并不相同。我见过有人直接把整个正北方向60度范围内的数据全部删掉理由是多数报告都这么写但实际上那支塔的风速计装在朝北一侧塔影区应该在正南。这种照搬模板的做法很危险直接削减了有效数据量还引入偏差。必须根据本塔的实际传感器几何布局来判断。第三件事是关于数据备份的版本管理。我习惯在每跑完一个清洗步骤后先把结果保存成中间文件而不是覆盖原文件。这样做的好处是如果有人质疑最终结果我可以回溯到每一步快速定位是哪一步的筛除范围偏大或者插值影响了统计。有一回我就是靠中间文件发现扇区剔除参数里有一处角度边界写错了导致尾流方向的几个扇区被错误剔除。要不是留有中间结果这个错误可能要等到报告写完才能发现。还有一个小而关键的细节计算空气密度的时候气压计海拔修正。测风塔上的气压传感器如果安装在塔基处读数对应的海拔与轮毂高度处的气压有差别。虽然气压随高度的衰减对密度的影响并不大不到明显误差级别但既然我们是做评估严谨一点总没坏处。我的处理是用测站海拔统一到一个参考面保证同一项目里不同测风塔的数据具有可比性。7. 这份评估的后续延展数据清洗完毕威布尔参数和风功率密度都算出来了并不代表风资源评估工作就结束了。气象塔数据只是点上的测量真正要给整个场址下一个结论还要结合中尺度模拟降尺度数据、地图地形分析、尾流模型测算以及周边已运行风电场的实际发电数据做交叉验证。如果是复杂山地地形一个测风塔的数据远不足以代表整个场址的各个机位就需要补充临时测风设备或者使用相关分析推算法进行空间推算。我个人实际操作下来的体会是数据处理过程的确定性其实很强工具用什么、代码怎么写都不是最难的关口。真正考验经验的是决定哪些数据可用、哪些不可用哪段插值合理、哪段必须整段作废。这些判断没有统一的算法公式靠的是对整个测量系统误差来源的深刻理解和对数据曲线特征的敏锐直觉。希望这次完整的处理流程拆解能让刚开始接触测风塔数据的朋友少走些弯路。最后再分享一个小技巧把处理好的timetable和最终的统计结果在项目结束后原样保留一份连同关键参数配置文件一起打包归档。以后如果项目要扩容或者做后评估这套数据能直接复用省下的时间远比你当时整理它所花的时间多得多。
阅读完成 · 觉得有帮助？
