做风资源评估的人，手里最重要的是什么？很多刚从学校出来的朋友会觉得是WAsP或者WindPRO这类高级软件，但真正落到“我手里只有一座70米测风塔一年的10分钟平均数据，要怎么判断这片区域值不值得继续往下做”的时候，最踏实的第一步反而是用Matlab把数据老老实实读进来，清洗干净，算出一批基础参数，画几张能说服自己的图。这篇文章就用一个典型的气象塔历史风力数据案例，从导入、处理到计算与出图，完整走一遍Matlab流程。无论你是刚入行的评估工程师、准备做毕业设计的研究生，还是临时接手测风数据项目的开发者，都可以照着这套思路把数据摸透。
为了避免光贴代码，我会把每个步骤背后的原理也讲清楚。你不需要事先掌握很多统计学知识，跟着操作就能理解为什么平均风向不能用普通mean，为什么风功率密度要用风速的三次方，为什么Weibull拟合对估算年发电量这么重要。
1. 风能资源评估的核心指标与气象塔数据基础
1.1 气象塔测风数据里到底有什么
测风塔不是简单的立一根杆，不同高度会安装不同传感器。常见测风高度是10m、30m、50m、70m，再高的塔还会测到90m、100m或更高，具体看项目需求。每个高度通常有一台风速计和一个风向标，塔身上还会配置温度、气压和湿度传感器。更专业一点的测风塔会加装超声风速仪，用来记录高频湍流数据。
对常规风资源评估来说，气象塔输出的数据一般是10分钟平均记录，其中包含该时段平均风速、最大风速、最小风速以及标准差。如果是1Hz采样，再聚合成10min均值，原始文件里可能还有每秒原始数据，但这类文件体积很大，通常只保留统计值。你看一个典型数据表，列名大概长这样：TIMESTAMP、WD10、WS10、WD30、WS30、TEMP70、PRESS70。时间列是采样时刻，WD和WS分别代表风向和风速，后面的数字是测量高度。
但真正打开这种文件之后，坑特别多。时间戳格式不统一、个别行列错位、风向标半夜冻住、风速传感器被扬尘堵住，这些问题我在项目里都遇到过。所以数据处理的第一步不是急着算数，而是先搞清楚传感器类型、采样原理、数据约定，再动手清洗。否则后面所有计算结果都建立在一个不牢靠的基础上。
1.2 风能资源评估的关键指标
风资源评估看起来指标很多，真正核心的可以压缩成下面这六项。
第一是平均风速。通常指某个高度上10min平均风速的算术平均，是最常用的资源描述。但平均风速不能单独使用，因为它没有反映风速分布。
第二是风功率密度。这是风能资源丰歉的核心指标，公式是0.5乘以空气密度再乘以风速三次方的平均。因为风机捕获的功率与风速三次方成正比，所以风功率密度比平均风速更能说明问题。很多初学者会用平均风速直接立方，这是错的，正确做法是对每个样本的风速三次方取平均，再乘以0.5和空气密度。
第三是Weibull分布参数。风速频率分布通常用双参数Weibull分布描述，形状参数k和尺度参数c。k决定了分布的形态和风速的集中程度，c则反映整体风速量级。机组发电量估算和极端风计算都离不开这两个参数。
第四是湍流强度，等于风速标准差除以平均风速，公式是TI=σ/v。湍流强度影响机组疲劳载荷和实际出力，在复杂地形尤其重要。
第五是主导风向。包括风向频率分布和风能密度方向分布，这决定了机位怎么排布、机组之间的尾流影响有多大。
第六是风切变指数。风速随高度变化，通常用幂律公式v2/v1=(h2/h1)^α来推算不同高度的风速，α就是风切变指数。它影响从测风高度向轮毂高度折算风速的可靠性。
1.3 数据的时间分辨率与质量控制前提
评估一个场址，测风数据至少要覆盖一个完整年，因为风资源的季节性非常强。夏天和冬天的平均风速可能相差一倍以上。如果项目只给了你三个月的数据，顶多能做初步判断，不能直接拍板年发电量。
时间分辨率上，10min平均是行业最常用的基线。1h平均会把阵风信息平滑掉，算出来的湍流强度会偏低；1s原始数据则用于湍流和极端风分析，文件量巨大。常规风资源评估，我建议至少保留10min平均数据，再做小时和日聚合分析。
质量控制必须放在所有统计计算之前，不能跳过。数据有效率低于90%不代表不能使用，关键要看缺失分布。如果一年里只有某一个月集中缺测，哪怕总有效率95%，也不能简单插补后当作完整年计算。野点剔除也需要规则，比如风速超出0~60m/s、风向超出0~360°、一个小时内风速突变超过阈值的点，都要打标签。气象塔数据永远不可能完美，但只要你把清洗规则记录下来，评估结论就可追溯，可复现。
2. Matlab数据导入与预处理实操
2.1 不同格式数据的导入方法比较
说回Matlab。我常用的导入方式是根据文件大小和格式来选。小规模Excel文件，readtable最省事；几十万行的CSV，readmatrix或textscan更快；如果文件自带可解析的时间列，readtimetable能一步到位。
先用readtable作为通用入门方案，代码是这样的：
% 读取CSV，保留原始列名 raw = readtable('met_tower.csv', 'PreserveVariableNames', true); % 把列名统一成后续方便使用的名称 % 假定原始列顺序：时间、风向、风速、温度、气压 raw.Properties.VariableNames = {'Time','WindDir','WindSpeed','Temp','Press'}; % 解析时间戳，常见格式 yyyy-MM-dd HH:mm:ss raw.Time = datetime(raw.Time, 'InputFormat', 'yyyy-MM-dd HH:mm:ss'); % 转成timetable，后面很多函数都围绕timetable设计 tt = table2timetable(raw); % 如果你有多个高度，比如WS10和WS50，用cell数组统一列名也可以 % raw.Properties.VariableNames = {'Time','WindDir','WindSpeed','Temp','Press','WS10','WS50'};如果文件时间列本身能被Matlab自动识别，直接用readtimetable更简单：
tt = readtimetable('met_tower.csv');但readtimetable对时间格式比较敏感，时间戳如果带有时区后缀“+08:00”，或者日期直接用“20230101”这种整数表示，它不一定能解析。此时还是要走readtable加datetime转换的老路。
Excel文件也不难处理，多个sheet分别对应不同高度时，最好用readtable分别读取，再按时间戳做outerjoin合并。我见过一个项目把70m风速放在一个sheet，温度放在另一个sheet，日期还不对齐，最后只能先分别读入，用时间戳join，把所有列拼在一起。这种时候千万不要凭行号对齐，时间戳才是唯一基准。
2.2 缺失值处理：插补还是剔除
导入数据后，第一件事是看缺失情况。用summary可以快速知道每列有多少missing值，也可以自己统计：
% 统计每列缺失数量 missCount = sum(ismissing(tt)); disp(missCount);缺失处理的原则是“零散缺失可插补，连续成片要剔除”。温度和气压存在明显的日周期，零散缺失用线性插补影响不大。但某几个连续小时的10min风速缺口，如果直接用线性插补填上，相当于人为制造了一段并不存在的持续风，会污染湍流强度计算。这种情况应该保留NaN，后续计算时用omitnan跳过。
我自己的经验阈值是：连续缺失超过6个10min点，也就是1小时，就不插补了。代码可以这么做：
% 先对零散缺失做线性插补 tt.WindSpeed = fillmissing(tt.WindSpeed, 'linear'); % 找出连续缺失超过6个点的缺口，重新置为NaN missVec = double(ismissing(tt.WindSpeed)); d = diff([0; missVec(:); 0]); starts = find(d == 1); ends = find(d == -1) - 1; runLengths = ends - starts + 1; % 只把长缺失段设为NaN longRuns = runLengths >= 6; for i = find(longRuns) tt.WindSpeed(starts(i):ends(i)) = NaN; end注意fillmissing默认是线性插补，它会完全用前后相邻点连一条直线。如果缺口只有两三个点，这种处理几乎不影响统计量。但如果你做湍流和极端阵风分析，连1个点的插补都要谨慎，最好保留原始标记。
2.3 异常值与野点剔除
异常值不能只靠“看着不对”来剔除，要有可解释的规则。第一道规则是物理合理性判断：
% 风速合理范围，单位m/s validWS = tt.WindSpeed >= 0 & tt.WindSpeed <= 60; % 风向合理范围，单位° validWD = tt.WindDir >= 0 & tt.WindDir <= 360; % 温度合理范围，单位°C validTemp = tt.Temp > -50 & tt.Temp < 60; % 组合成综合有效标记 valid = validWS & validWD & validTemp; tt.WindSpeed(~valid) = NaN;第二道规则是针对突变点。常用的两种方法：滑动中位数和变化率阈值。滑动中位数适合捉孤立的尖峰，变化率阈值适合抓传感器跳变。
% 用1小时滑动中位数判断离群点 out = isoutlier(tt.WindSpeed, 'movmedian', hours(1)); tt.WindSpeed(out) = NaN; % 10min平均风速前后变化超过10m/s，属于异常跳变 dV = abs([NaN; diff(tt.WindSpeed)]); tt.WindSpeed(dV > 10) = NaN;这里有个很容易犯的错：不要直接把大于“mean+3*std”的值全删掉。风速本身是右偏分布，强风过程完全属于正常现象，3σ会把正常的极端大风误删。物理范围和变化率才是更安全的过滤器。
风向传感器还有一个典型故障：冻结。温度低于0°C时，风向标可能被冰卡住，输出长时间保持同一个角度。处理方法是检查每个风向值连续重复的次数，超过一定数量就标记为异常。风速和气温过低可以作为辅助判断条件。这类问题只靠自动算法不一定抓得准，最好再叠加人工抽查时间序列。
3. 风能资源参数计算与Matlab代码实现
3.1 平均风速与风向的统计
预处理完成后，先做最基础的统计。Matlab里的retime函数可以把10min数据聚合成小时平均或日平均，在分析季节性变化时非常方便：
% 小时平均风速 ws_hourly = retime(tt(:, 'WindSpeed'), 'hourly', 'mean'); % 日平均风速 ws_daily = retime(tt(:, 'WindSpeed'), 'daily', 'mean'); % 年平均风速 ws_yearly = retime(tt(:, 'WindSpeed'), 'yearly', 'mean');风向就不能这么算了。风向是角度，0°和360°本身是同一个方向。如果你直接对350和10取平均，得到的是180°，完全错误。正确做法是把风向当成圆上的向量，对sin和cos分别取平均，再通过atan2还原角度：
% 风向圆平均 u = mean(sin(deg2rad(tt.WindDir)), 'omitnan'); v = mean(cos(deg2rad(tt.WindDir)), 'omitnan'); dirMean = rad2deg(atan2(u, v)); if dirMean < 0 dirMean = dirMean + 360; end这样做350和10的平均就会接近0°，而不是180°。除了平均方向，还可以用合成向量长度R来量化风向稳定性，R接近1说明风向非常集中，R接近0说明风向散乱。这个R值在报告里作为风向稳定性指标也很有用。
3.2 风功率密度与空气密度修正
风功率密度是整个评估里最不能算错的一个数。空气密度ρ并不是常数，它随气温、气压和湿度变化。标准海平面条件下大约是1.225kg/m³，但在海拔3000m的高原，密度可能只有0.9kg/m³，直接影响风功率密度和发电量。
常规做法是用温度和气压反算密度：
% 假设温度单位是°C，气压单位是kPa Tk = tt.Temp + 273.15; Ppa = tt.Press * 1000; rho = Ppa ./ (287.05 .* Tk); % 风功率密度，逐个样本计算后取平均，单位W/m² wpd = mean(0.5 .* rho .* tt.WindSpeed.^3, 'omitnan');注意这里用的是逐样本的rho乘以v³再取平均，不是用平均rho乘以平均风速的三次方。这两个结果差异很大，因为风速立方放大了大风时段的贡献。我在项目里见过有人把平均风速代入v³，算出来的风功率密度低估了20%以上。风功率密度对高风速段极其敏感，数据里只要有几个强风过程，就能明显抬升整体值，如果用平均风速替代，恰恰把最关键的贡献丢掉了。
3.3 Weibull分布参数估计与拟合质量判断
风速频率分布通常用双参数Weibull分布描述。Matlab里的fitdist可以直接做极大似然估计：
% 取出有效风速，忽略NaN ws_valid = rmmissing(tt.WindSpeed); % Weibull分布拟合 pd = fitdist(ws_valid, 'Weibull'); c = pd.A; % 尺度参数 k = pd.B; % 形状参数如果用的Matlab版本比较老，可以用wblfit，返回参数顺序也是尺度参数在前，形状参数在后，别搞反了。另一个我常用的快速估计方法是矩估计：
mu = mean(ws_valid); sigma = std(ws_valid); k0 = (sigma/mu)^(-1.086); c0 = mu / gamma(1 + 1/k0);这个公式在绝大多数情况下给出的初值都很稳，也可以用来核对fitdist的结果。k通常在1.5到3.5之间，c大约等于平均风速的1.05到1.15倍。如果你拟合出来的k只有0.8或者高达6，就要警惕数据分布是不是被异常值污染了。
参数算出来不是终点。要把频率直方图和拟合曲线叠加，重点看大风速段是否拟合合理。如果拟合曲线在高风速端明显低于柱状图，说明尾段被低估，年发电量可能会被低估；反之高估则要检查野点。做这种判断时，最好配合年发电量计算一起看，单纯看曲线贴合程度是不够全面的。
3.4 风玫瑰图与主导风向识别
风玫瑰图是判断主导风向最直观的工具。Matlab里的polarhistogram可以快速画出风向频率图：
dirRad = deg2rad(tt.WindDir); figure; polaraxes('ThetaZeroLocation', 'top', 'ThetaDir', 'clockwise'); polarhistogram(dirRad, 16, 'Normalization', 'probability'); title('风向频率玫瑰图');这段代码把北放在正上方，角度按顺时针方向递增，符合国内气象报告的习惯。16个扇区每个跨22.5°，是风资源报告里最常见的扇区划分方式。
风能玫瑰比频率玫瑰更实用，因为它反映了哪个方向传来的风真正携带更多的能量。polarhistogram不好直接支持加权概率，我通常在Matlab里算扇区统计表，再决定画图方式：
edges = 0:22.5:360; dirBin = discretize(tt.WindDir, edges); ok = ~isnan(dirBin) & ~isnan(tt.WindSpeed); % 每个扇区风向频率 freq = accumarray(dirBin(ok), 1, [16 1]) / sum(ok) * 100; % 每个扇区平均风速 meanWS = accumarray(dirBin(ok), tt.WindSpeed(ok), [16 1], @nanmean, NaN); % 每个扇区风能密度 E = 0.5 .* rho .* tt.WindSpeed.^3; meanE = accumarray(dirBin(ok), E(ok), [16 1], @nanmean, NaN); pctE = 100 .* meanE ./ nansum(meanE); T = table(edges(1:end-1)' + 11.25, freq, meanWS, pctE, ... 'VariableNames', {'扇区中心方向', '频率百分比', '平均风速', '风能占比'}); disp(T);这张统计表本身就可以写进报告。要注意的是扇区边界不要重复，0°和360°其实同一个方向，discretize如果直接用0:22.5:360作为边界，360°会被分到最后一个扇区，容易多算。稳妥的做法是把所有等于360°的风向先改成0，再参与分扇区。
3.5 湍流强度与风切变
湍流强度是风速标准差除以平均风速。如果你的气象塔数据文件已经包含每10min的风速标准差列，直接相除即可：
% 假设有标准偏差列 WS_std TI = tt.WS_std ./ tt.WindSpeed; % 风速越低，TI越不稳定，一般只统计风速大于4m/s时的湍流强度 sel = tt.WindSpeed > 4; TI_mean = mean(TI(sel), 'omitnan');风速接近0时，TI会极端偏大，这种情况下统计出来的平均值没有工程意义，所以要加一个风速下限。如果你的数据没有标准差列，只能用原始高频数据算，或者用小时标准差代替，但结果会偏小，报告里必须注明。
风切变指数用两高度风速计算。假设有10m和50m风速，按下式求α：
v10 = mean(tt.WS10, 'omitnan'); v50 = mean(tt.WS50, 'omitnan'); alpha = log(v50 / v10) / log(50 / 10);切变指数通常在0.1到0.5之间。白天太阳加热、大气混合充分，α接近0.1；夜间地面冷却、大气层结稳定，α可能超过0.3。如果项目关心轮毂高度风速，最好分季节和昼夜统计α，用同一个年均值推算可能产生明显偏差。切变指数还和地形、地表粗糙度密切相关，复杂地形区域更要用实测数据统计，不能随便套公式。
4. 结果可视化和评估报告输出
4.1 时间序列图与高度对比
评估报告里，图表比表格更能直接说明问题。最基础的时间序列图要解决两件事：整体趋势和细节突变。拿一整年10min数据直接画一条曲线会密成一片黑，什么也看不清。所以我通常先选前30天看细节，再用月平均或日平均图看长期趋势。
% 前30天，大约4320个10min点 figure; plot(tt.Time(1:4320), tt.WindSpeed(1:4320), 'LineWidth', 0.8); hold on; plot(tt.Time(1:4320), tt.WS50(1:4320), 'LineWidth', 0.8); legend('70m风速', '50m风速', 'Location', 'best'); xlabel('时间'); ylabel('风速（m/s）'); grid on;全年数据改成月平均柱状图更清楚，一眼能看出哪几个月风大、哪几个月风小。箱线图也能看每个月的风速分布范围：
% 提取月份 tt.Month = month(tt.Time); % 绘制月风速箱线图 figure; boxchart(tt.Month, tt.WindSpeed); xlabel('月份'); ylabel('风速（m/s）');如果你的Matlab版本没有boxchart，可以用传统的boxplot代替，效果类似。重点在于报告里除了“年平均风速”，还要有季节变化的描述。
4.2 风速频率直方图与Weibull叠加
频率直方图和Weibull拟合曲线叠加，是风资源报告里最常出现的图，也是最能看出拟合问题的图。前面已经算出了c和k，画图这一步很直接：
ws_valid = rmmissing(tt.WindSpeed); pd = fitdist(ws_valid, 'Weibull'); c = pd.A; k = pd.B; figure; histogram(ws_valid, 'Normalization', 'pdf', 'NumBins', 40, ... 'FaceColor', [0.7 0.7 0.7]); hold on; vline = linspace(0.5, max(ws_valid), 100); plot(vline, wblpdf(vline, c, k), 'r-', 'LineWidth', 2); legend('实测风速概率密度', 'Weibull拟合'); xlabel('10min平均风速（m/s）'); ylabel('概率密度'); grid on;观察曲线时，重点看高风速尾段。如果实测柱状最右边还有不少超过20m/s或30m/s的点，而拟合曲线已经衰减到接近0，说明Weibull分布没有完全抓住强风过程。这时候不要急着手动改参数，先检查这些大风是不是异常值；确认无误后，可以在报告中单独指出强风过程对评估的影响，必要时通过极值分布模型补充描述。
4.3 风玫瑰图的美化与导出
风玫瑰图画完后，导出图片时要注意分辨率和尺寸。Matlab里exportgraphics是一个很稳定的选择：
figure; polaraxes('ThetaZeroLocation', 'top', 'ThetaDir', 'clockwise'); polarhistogram(deg2rad(tt.WindDir), 16, 'Normalization', 'probability'); title('风向频率玫瑰图'); % 导出300dpi的PNG exportgraphics(gcf, 'wind_rose.png', 'Resolution', 300);风能玫瑰的绘制稍微麻烦一点。polarhistogram不支持权重，我通常把扇区统计表算好后，在Matlab里用polaraxes加自定义扇形填充，或者直接导出表后用其他绘图工具画。如果你愿意折腾，polaraxes底层的cartesian坐标转换能实现扇形填充，但代码量不小。更建议的做法是：先输出上文中的T表格，再根据表格画图，这样统计与绘图解耦，数据核对也更容易。
4.4 一键生成评估表格
最后把关键参数汇总成一个表，方便写报告，也方便项目交接。不要在分析结束后临时拼凑，从导入那天起就建立汇总变量，每算完一个指标就写入一行，后面引用会非常方便。
result = table; result.meanWS = mean(tt.WindSpeed, 'omitnan'); result.maxWS = max(tt.WindSpeed, [], 'omitnan'); result.wpd = wpd; result.TI_mean = TI_mean; result.k = k; result.c = c; result.dirMean = dirMean; writetable(result, 'wind_resource_summary.csv');如果项目有多个测风高度，把高度信息也加到表里，比如多一行“测风高度=70m”。这样汇总表本身就是一个小型数据库。导出图片时，统一用exportgraphics设置300dpi，保存为PNG，插入Word或PDF报告都很清晰。
5. 实操中的常见问题与避坑经验
5.1 塔影效应和传感器故障识别
塔影是测风塔支架和传感器之间相互作用产生的测量偏差。风从塔支架下游方向吹来时，风速会被扭曲，典型特征是某个固定扇区的平均风速明显低于相邻扇区。这个问题不画风玫瑰很难被发现。
处理方法是识别受扰扇区并剔除该扇区数据。常见的受扰扇区是塔体正对来流的下风向区域，范围大约30°到60°。具体做法是绘制不同扇区的平均风速图，找持续低值扇区，再对照现场照片或塔结构图确认受扰方向，最后在统计计算时把该扇区数据标记为无效。有些评估软件内部也有扇区有效率概念，但如果你在用Matlab处理，必须自己完成这个过程。
5.2 Weibull拟合的收敛性问题
风速数据如果大量集中在0附近，fitdist可能出现Warning，甚至拟合出非常不合理的k值。这是因为Weibull分布的对数似然函数在风速等于0时会出现极端值，数据里的静风点越多，影响越大。
一个稳的办法是自定义负对数似然函数，用mle优化并设置参数边界和初值：
ws_valid = rmmissing(tt.WindSpeed); ws_valid = ws_valid(ws_valid > 0.1); % 过滤静风点，但记录数量变化 nlogLL = @(params) -sum(log(wblpdf(ws_valid, params(1), params(2)))); param0 = [c0, k0]; opts = statset('MaxIter', 5000); phat = mle(ws_valid, 'pdf', @wblpdf, 'start', param0, ... 'Options', opts, 'LowerBound', [0.01, 0.2], ... 'UpperBound', [100, 10]); c = phat(1); k = phat(2);过滤掉静风点会影响发电量估算，所以要在报告中注明过滤比例，特别是静风占比高的地区。更稳妥的做法是把静风时段单独统计，不参与Weibull拟合，但后续发电量计算时再把静风小时数加回来。
5.3 风向圆统计的边界陷阱
前面提过风向平均不能用普通mean。这里再强调一个容易踩的坑：当定义扇区边界时，0°和360°是同一个方向，但程序会把它们当成两个数。如果你用0:22.5:360作为边界，discretize会把0分到第一个扇区，360分到最后一个扇区，相当于最后多算了一次北极向。
解决办法很简单，先把所有风向修正到[0, 360)区间：
wd = mod(tt.WindDir, 360); wd(wd == 0) = 360; % 如果希望以360代表正北也可以，但要保证边界统一实际处理时要固定一种约定，并在脚本注释里写清楚，否则下一个人接手时很难看懂。这种细节最容易在多人协作时引入错误，所以代码注释不要省。
5.4 数据代表性不足怎么办
如果测风塔数据只有八个月，不能直接把评估结果当作长期资源结论。常规做法是找附近长年代气象站或再分析数据的同期风速序列，和测风塔数据做线性回归，再把测风塔数据订正到长周期平均风速。这个流程叫MCP，本质上是把短期测量“挂靠”到长期参考序列上。
MCP听起来高级，实际做起来有几个前提。参考站距离不能太远，地形特征不能差太多，参考序列质量也要检查。再分析数据用来补空间覆盖可以，但直接替代现场测站风险很大，因为再分析数据的分辨率不一定能反映局地地形影响。
所以拿到一份历史测风数据时，先不要急着写代码，先问清楚数据的时间范围、原始采样频率、维护记录和故障标记。这些背景信息比任何算法都重要。数据文件本身不会告诉你传感器坏过，只有维护日志会。
5.5 报告里的“平均风速”陷阱
很多评估报告只给年平均风速，这远远不够。平均风速相同的两个场址，风功率密度可能差30%，因为风功率密度对高风速段极其敏感。正确做法是同时给出年平均风速、风功率密度、Weibull参数、湍流强度、主风向扇区、静风小时占比。
另外，季节变化一定要展示。单一个年均值会掩盖“冬季大风、夏季小风”或“春秋稳定、夏冬波动”这些重要特征。报告评审时，一张月平均风速表和一张风向频率表往往比一段文字更有说服力。
最后再分享一个我自己早期的教训：有一次我在算湍流强度时，直接用全套数据求标准差和均值，结果12月晚上几个夜间强风过程把TI平均值拉得很高。后来把数据按月、按昼夜分组，发现真正异常的是那几天传感器结冰。从那之后我养成了“先分组观察、再整体统计”的习惯，也建议你们都保留一个清洗日志。气象塔数据不可能完美，但只要你把每一步清洗规则和剔除量记录下来，评估结论就是可追溯、可复现的。这比任何高级算法都重要。
