企业官网建设流程全解析

简介：本资源为风能资源评估方向的Matlab代码工具包，面向计算机、电子信息工程、数学等专业的大学生及科研入门者，可用于课程设计、期末大作业与毕业设计。包内共25个文件，以19个m脚本为核心，辅以mat、txt、accdb与kml等数据与地理文件，压缩包约5.99MB，覆盖风速分布分析、风功率密度计算、风能玫瑰图绘制、数据质量控制与长期风数据抓取等模块。代码采用参数化编程，关键参数可灵活调整，注释清晰，并附赠案例数据可直接运行，便于快速搭建评估模型并验证结果。已有66人学习，适合需要完整实现思路与可复用脚本的读者参考借鉴。

1. 风能资源评估的 MATLAB 实现：从测风数据到可投资结论

拿到一份测风塔的原始数据，直接算个平均风速就敢写资源评估报告，这事我干过，后来被现实狠狠教育。风能资源评估远不止“平均风速”四个字，它要回答的是：这个点位一年能发多少电、风机的选型是否匹配、湍流和切变会不会让机组提前报废、50 年一遇的极端风速能不能扛住。MATLAB 在这套流程里的价值，是把 Weibull 拟合、风功率密度、风频玫瑰图、湍流强度、极大风速推算这些环节串成一条可复现的流水线，而不是每次换一个测风塔就重写一遍脚本。这篇文章面向的是手里有测风数据、需要出评估结论的风资源工程师和新能源方向的研究生，我会把每个环节的代码、参数含义和踩过的坑都摊开讲，让你拿到数据就能跑通。

2. 测风数据的预处理与质量控制：脏数据不清理，后面全是玄学

风资源评估的第一步永远不是建模，而是把测风塔导出的原始数据洗干净。测风塔通常按 10 分钟或 1 小时记录，字段包括各高度层的风速、风向、温度、气压，还有一堆状态码。原始数据里常见的脏东西有：传感器结冰导致的持续零风速、通讯中断造成的整段缺失、风向在 0/360 度附近的跳变、以及超出物理量程的野值。这些不处理，后面 Weibull 拟合出来的 k 和 c 参数会偏得离谱，年发电量估算能差出 20% 以上。

2.1 读入与字段对齐

我一般先把数据统一成 MATLAB 的 timetable，时间列做严格校验，因为不同测风塔导出的时间戳格式五花八门，有的是本地时间，有的是 UTC，混用会让风频玫瑰图整体旋转。

% 读入测风塔导出的 CSV，假设列顺序为：时间, 10m风速, 30m风速, 50m风速, 风向, 温度 opts = detectImportOptions('tower_raw.csv'); opts.VariableNamingRule = 'preserve'; T = readtimetable('tower_raw.csv', opts); % 统一时间列为 datetime 并设为行时间 T.Properties.DimensionNames{1} = 'Time'; T = sortrows(T, 'Time'); % 按时间排序，防止乱序 T = rmmissing(T, 'DataVariables', ... {'WindSpeed_10m','WindSpeed_30m','WindSpeed_50m'}); % 先删关键列缺失 % 物理量程过滤：风速合理范围 0~40 m/s，超出视为野值 validIdx = T.WindSpeed_50m >= 0 & T.WindSpeed_50m <= 40; T = T(validIdx, :);

这段代码的关键点有三个。readtimetable保证时间列被识别为 datetime 而不是字符串，后续做时间重采样才不会出错。sortrows是血泪经验，有些数据采集软件在断电恢复后会把缓存数据追加到文件末尾，时间顺序是乱的。物理量程过滤用 40 m/s 作为上限，是因为绝大多数陆上测风塔的杯式风速仪量程就在这个范围，超过的基本是电磁干扰产生的野值。

2.2 缺失值插补与结冰标记

缺失值不能简单删掉，否则会破坏时间序列的连续性，影响后续湍流强度的计算。常见做法是：连续缺失小于 3 个点的用线性插值，大于 3 个点的整段标记为无效，不参与统计。结冰的识别靠温度和风速的联合判断——温度低于 2 摄氏度且风速持续为 0，大概率是传感器结冰而不是真的无风。

% 线性插值填补短缺失（最多 3 个连续点） T.WindSpeed_50m = fillmissing(T.WindSpeed_50m, 'linear', 'MaxGap', 3); % 结冰标记：温度 < 2 且风速 = 0，标记为无效 iceFlag = T.Temperature < 2 & T.WindSpeed_50m == 0; T.WindSpeed_50m(iceFlag) = NaN; % 统计有效数据完整率 validRate = sum(~isnan(T.WindSpeed_50m)) / height(T); fprintf('50m 风速有效数据完整率：%.2f%%\n', validRate * 100);

MaxGap参数设成 3 是经验值，对应 10 分钟数据里的半小时缺口。完整率低于 90% 的测风塔，我一般会建议延长观测周期或者换塔，因为 IEC 61400-12 对数据完整率有明确要求，低于这个阈值做出来的年发电量估算置信度不够。结冰标记这一步很多人会漏，结果冬季月份的平均风速被拉低，Weibull 分布的尺度参数 c 偏小，最后算出来的发电量偏保守，投资决策就会误判。

2.3 风向数据的环形处理

风向是环形数据，0 度和 360 度是同一个方向。做风频玫瑰图时如果直接按数值分箱，0 度附近的数据会被拆到首尾两个箱子里，玫瑰图会出现一个假的缺口。正确做法是用circ_系列函数或者手动做角度折叠。

% 风向折叠到 [0, 360) T.WindDirection = mod(T.WindDirection, 360); % 按 16 方位分箱，每个方位 22.5 度 edges = 0:22.5:360; binIdx = discretize(T.WindDirection, edges, 'IncludedEdge', 'left'); binIdx(isnan(binIdx)) = 16; % 360 度归入第 16 个方位 % 统计各方位风频 freq = histcounts(binIdx, 1:17) / sum(~isnan(binIdx));

discretize的IncludedEdge参数设成'left'是为了让每个区间的左边界闭合，避免 0 度被重复计数。最后那个binIdx(isnan(binIdx)) = 16是处理 360 度整点的情况，因为discretize对等于最后一个边界的值会返回 NaN。这些细节看起来琐碎，但风频玫瑰图是给业主看的第一张图，画错了会直接影响对主导风向的判断，进而影响风机排布方案。

3. Weibull 分布拟合与风功率密度计算：参数选错，发电量差一个档

风资源评估的核心是把风速的统计分布拟合出来，最常用的是两参数 Weibull 分布。形状参数 k 反映风速的集中程度，尺度参数 c 反映风速的整体大小。这两个参数决定了风功率密度，而风功率密度直接决定了这个点位值不值得开发。很多人直接用wblfit一把梭，但测风数据里有效样本量、分箱方式、拟合方法的选择都会影响结果，下面把几种常见做法和适用场景讲清楚。

3.1 三种拟合方法的对比与选择

MATLAB 里拟合 Weibull 分布至少有三种路径：最大似然估计wblfit、矩估计、以及最小二乘法拟合累积分布。最大似然估计在样本量大于 1000 时最稳定，矩估计计算快但对野值敏感，最小二乘法适合数据已经分箱成风频直方图的情况。

% 方法一：最大似然估计（推荐，样本量充足时） pd = fitdist(T.WindSpeed_50m(~isnan(T.WindSpeed_50m)), 'Weibull'); k_mle = pd.A; % 形状参数 c_mle = pd.B; % 尺度参数 % 方法二：矩估计（快速估算，用于初步筛查） v_mean = mean(T.WindSpeed_50m, 'omitnan'); v_std = std(T.WindSpeed_50m, 'omitnan'); k_mom = (v_std / v_mean)^(-1.086); c_mom = v_mean / gamma(1 + 1/k_mom); % 方法三：最小二乘拟合累积分布（数据已分箱时） [sortedV, idx] = sort(T.WindSpeed_50m(~isnan(T.WindSpeed_50m))); empCDF = (1:length(sortedV))' / length(sortedV); X = log(sortedV); Y = log(-log(1 - empCDF)); p = polyfit(X, Y, 1); k_lsq = p(1); c_lsq = exp(-p(2) / p(1));

三种方法的结果通常不会完全一致，k 值差 0.1、c 值差 0.3 是常有的事。我的习惯是：先用矩估计快速看一眼，再用最大似然估计出正式结果，最后用最小二乘法交叉验证。如果三种方法算出来的风功率密度差异超过 5%，说明数据本身有问题，得回头查预处理环节。fitdist返回的pd.A是形状参数，pd.B是尺度参数，这个命名和 Weibull 分布的标准写法是反的，第一次用很容易搞混。

3.2 风功率密度与可开发等级判定

风功率密度是评估点位好坏的核心指标，计算公式里用到了 gamma 函数。根据国家标准，风功率密度等级从 1 级到 7 级，3 级以上才具备较好的开发价值。

% 空气密度修正（按实测温度和气压） rho = mean(T.Pressure ./ (287.05 * (T.Temperature + 273.15)), 'omitnan'); % 风功率密度计算 WPD = 0.5 * rho * c_mle^3 * gamma(1 + 3/k_mle); % 年发电量初步估算（假设风机扫风面积 A，功率系数 Cp） A = 12000; % 扫风面积，单位 m^2，对应约 4MW 机组 Cp = 0.45; % 功率系数，现代机组典型值 hours = 8760; E_annual = WPD * A * Cp * hours / 1000; % 单位 kWh fprintf('空气密度：%.3f kg/m^3\n', rho); fprintf('风功率密度：%.1f W/m^2\n', WPD); fprintf('年发电量估算：%.0f kWh\n', E_annual);

空气密度这一项很多人直接取 1.225 kg/m³，但在高原或者高温地区，实际空气密度可能只有 1.0 甚至更低，风功率密度会差出 20%。gamma(1 + 3/k_mle)这一项是 Weibull 分布三阶矩的解析解，不能省。年发电量估算里的Cp取 0.45 是当前主流机组的典型值，但实际应该查具体机型的功率曲线，这里只是初步筛查用。如果算出来的风功率密度低于 150 W/m²，我一般会直接建议放弃这个点位，除非有特殊的上网电价政策。

3.3 风频玫瑰图与主导风向分析

风频玫瑰图不只是画着好看，它决定了风机排布的主轴方向。主导风向对应的风能方向，机组应该尽量垂直于这个方向排列，减少尾流影响。

% 计算各风向方位的风频和平均风速 dirEdges = 0:22.5:360; dirIdx = discretize(T.WindDirection, dirEdges, 'IncludedEdge', 'left'); dirIdx(isnan(dirIdx)) = 16; freqByDir = zeros(1, 16); meanSpdByDir = zeros(1, 16); for i = 1:16 mask = dirIdx == i; freqByDir(i) = sum(mask) / sum(~isnan(dirIdx)); meanSpdByDir(i) = mean(T.WindSpeed_50m(mask), 'omitnan'); end % 绘制风频玫瑰图 theta = deg2rad(dirEdges(1:16) + 11.25); figure; polarplot([theta, theta(1)], [freqByDir, freqByDir(1)], 'b-o', 'LineWidth', 1.5); title('50m 高度风频玫瑰图');

polarplot画玫瑰图时，角度需要取每个方位区间的中心值，所以加了 11.25 度的偏移。风频和平均风速要分开统计，因为主导风向不一定等于最大风能方向——有时候某个方向风频不高但风速很大，风能贡献反而更大。这个区别在排布方案里很关键，尾流计算要用风能玫瑰图而不是风频玫瑰图。

4. 湍流强度与极端风速推算：决定机组选型和塔筒寿命

风资源评估报告里如果只有平均风速和风功率密度，那是不完整的。湍流强度和极端风速直接关系到机组的安全等级和塔筒的疲劳寿命。湍流强度超标，机组会频繁触发振动保护，发电量损失不说，关键部件寿命也会大幅缩短。极端风速算错，塔筒设计载荷不够，后果更严重。

4.1 湍流强度的计算与分级

湍流强度定义为风速标准差与平均风速的比值，通常按 10 分钟窗口计算，再取 90% 分位数作为代表值。IEC 61400-1 把湍流强度分为 A、B、C 三个等级，对应不同的设计标准。

% 按 10 分钟窗口计算湍流强度 T10 = retime(T, 'regular', 'mean', 'TimeStep', minutes(10)); T10.TI = T10.WindSpeed_50m_std ./ T10.WindSpeed_50m; % 取 90% 分位数作为代表湍流强度 TI_rep = prctile(T10.TI(~isnan(T10.TI)), 90); % 按风速分箱统计湍流强度（IEC 要求 15 m/s 处的值） binEdges = 0:2:30; spdBin = discretize(T10.WindSpeed_50m, binEdges); TI_by_bin = zeros(1, length(binEdges)-1); for i = 1:length(binEdges)-1 mask = spdBin == i; if sum(mask) > 10 TI_by_bin(i) = prctile(T10.TI(mask), 90); end end

retime做 10 分钟重采样时，标准差需要单独算，不能直接用'mean'。我一般会写一个自定义函数传给retime，同时输出均值和标准差。prctile取 90% 分位数是 IEC 规定的做法，不能用平均值代替，因为湍流强度的分布是右偏的，平均值会低估极端情况。如果 15 m/s 处的湍流强度超过 0.16，这个点位就需要选 IEC A 类机组，机型选择范围会窄很多，造价也更高。

4.2 极大风速的极值推算

极端风速推算用极值理论，常用 Gumbel 分布拟合年最大风速序列。测风塔通常只有 1 到 2 年的数据，直接拟合年最大风速样本量不够，常见做法是用月度最大值或者用风速的日最大值序列做拟合。

% 提取各月最大风速 T_monthly = retime(T, 'monthly', @(x) max(x, [], 'omitnan')); monthlyMax = T_monthly.WindSpeed_50m; % Gumbel 分布拟合 pdGumbel = fitdist(monthlyMax(~isnan(monthlyMax)), 'ExtremeValue'); mu_g = pdGumbel.mu; sigma_g = pdGumbel.sigma; % 推算 50 年一遇极大风速 T_return = 50; V_50 = mu_g - sigma_g * log(-log(1 - 1/T_return)); fprintf('50 年一遇极大风速（50m）：%.1f m/s\n', V_50);

Gumbel 分布的mu是位置参数，sigma是尺度参数。log(-log(1 - 1/T_return))是 Gumbel 分布的分位数函数，T_return是重现期。这里用月度最大值而不是年最大值，是因为测风塔观测周期通常不够长，月度最大值能把样本量扩大 12 倍。但要注意，月度最大值序列不满足极值理论的独立性假设，算出来的结果会偏保守，实际工程中会再乘一个 0.9 左右的折减系数。这个折减系数没有统一标准，各家咨询公司用的值不一样，报告里要写清楚。

4.3 风切变与轮毂高度风速外推

测风塔通常只测到 50m 或 70m，而现代机组的轮毂高度在 100m 以上，需要做风切变外推。常用的是幂律模型，风切变指数 alpha 由两个高度的风速比取对数得到。

% 计算风切变指数（用 30m 和 50m 数据） mask = ~isnan(T.WindSpeed_30m) & ~isnan(T.WindSpeed_50m) & ... T.WindSpeed_30m > 2 & T.WindSpeed_50m > 2; alpha = log(mean(T.WindSpeed_50m(mask)) / mean(T.WindSpeed_30m(mask))) / ... log(50 / 30); % 外推到 100m 轮毂高度 V_hub = mean(T.WindSpeed_50m, 'omitnan') * (100 / 50)^alpha; fprintf('风切变指数：%.3f\n', alpha); fprintf('100m 轮毂高度平均风速：%.2f m/s\n', V_hub);

风切变指数 alpha 的典型值在 0.1 到 0.3 之间，小于 0.1 说明大气层结不稳定，大于 0.3 说明可能有逆温层。mask里限制风速大于 2 m/s 是为了排除静风时段，因为低风速下的风向和风速测量误差都很大，会污染 alpha 的计算。外推到 100m 时，如果 alpha 是用 30m 和 50m 算的，外推距离比较远，误差会放大，最好用 50m 和 70m 的数据算。如果测风塔只有一层数据，那风切变只能取经验值，报告里要注明这个不确定性。

5. 避坑与排查：那些让评估报告翻车的细节

风资源评估这活，理论都不难，难的是数据里藏着的各种坑。下面这几条是我这些年踩过的、也见过别人踩的，每条都按现象、原因、解决来写，你对着排查能省不少时间。

5.1 现象：Weibull 拟合的 k 值大于 3，风功率密度异常偏高

原因通常是数据里混入了大量低风速的无效记录，比如传感器故障时的持续零风速，或者静风时段的重复计数。这些低风速值会把分布的形状参数拉大，c 值也跟着偏大。解决方法是回到预处理环节，检查零风速的占比，如果超过 20%，要逐条排查是真实静风还是传感器问题。另外，用fitdist之前先画一下风速直方图，看看分布形状是不是双峰，双峰说明数据里混了不同来源的记录。

5.2 现象：风频玫瑰图的主导风向和附近测风塔完全相反

原因多半是风向传感器的安装方向错了，或者数据导出时做了镜像翻转。有些测风塔的风向传感器在安装时没有对准正北，而是按塔筒的某个标记对齐，导致整体偏移一个固定角度。解决方法是找同一区域其他测风塔的数据做交叉验证，如果偏移是固定的，可以整体旋转修正。如果偏移不固定，那可能是传感器故障，这段数据不能用。

5.3 现象：年发电量估算比实际运行数据高 30% 以上

原因通常是尾流模型没考虑，或者空气密度用了标准值而实际点位海拔较高。尾流损失在大型风电场里能到 10% 到 20%，单台机组的估算不能直接乘以机组数量。解决方法是先用 Jensen 模型或者更精细的 CFD 做尾流计算，再乘一个综合折减系数，包括可利用率、电气损耗、叶片污染等。空气密度必须用实测温度和气压算，不能偷懒取 1.225。

5.4 现象：极端风速推算结果每次跑都不一样

原因是极值拟合对样本量很敏感，月度最大值序列里如果有一两个异常大的值，拟合结果就会跳。解决方法是先做异常值检验，用 Grubbs 检验或者箱线图法把离群点挑出来单独看。如果离群点是真实的极端天气事件，要保留；如果是传感器野值，要剔除。另外，重现期推算的置信区间也要算出来，只给一个点估计值是不严谨的。

5.5 现象：MATLAB 脚本在换了一台电脑后报错

原因通常是路径分隔符或者编码问题。Windows 和 Linux 的路径分隔符不一样，中文路径在某些 MATLAB 版本里也会出问题。解决方法是所有路径用fullfile拼接，不要手写斜杠。数据文件尽量用英文命名，CSV 的编码用 UTF-8 无 BOM。如果用了自定义函数，确保函数文件在 MATLAB 的搜索路径里，或者用addpath显式添加。

6. 从评估到报告：把 MATLAB 结果变成可交付的结论

跑完上面所有计算，你手里会有一堆数字：k、c、风功率密度、湍流强度、极大风速、年发电量。但这些数字本身不是报告，报告要回答的是“这个点位能不能投、投多少、风险在哪”。我一般会把 MATLAB 的输出整理成三张核心表：风资源参数汇总表、风频风能分布表、以及不同机型方案的年发电量对比表。下面这段代码把关键结果导出成结构化的表格，方便后续写报告时直接引用。

% 汇总关键参数 results = table( ... {'形状参数 k'; '尺度参数 c'; '平均风速 (m/s)'; '风功率密度 (W/m^2)'; ... '湍流强度 (90%)'; '50年一遇极大风速 (m/s)'; '风切变指数'; '年发电量估算 (kWh)'}, ... [k_mle; c_mle; mean(T.WindSpeed_50m,'omitnan'); WPD; TI_rep; V_50; alpha; E_annual], ... 'VariableNames', {'Parameter', 'Value'}); % 导出为 CSV，供报告引用 writetable(results, 'wind_resource_summary.csv'); % 不同机型方案对比（假设三种轮毂高度和功率曲线） hubHeights = [90, 100, 110]; ratedPower = [3.6, 4.0, 4.5]; % MW for i = 1:length(hubHeights) V_h = mean(T.WindSpeed_50m,'omitnan') * (hubHeights(i)/50)^alpha; % 简化估算，实际应查功率曲线插值 CF = 0.35; % 容量系数，需根据功率曲线和风速分布计算 AEP = ratedPower(i) * 1000 * CF * 8760; fprintf('轮毂高度 %dm，额定功率 %.1fMW，估算年发电量 %.0f kWh\n', ... hubHeights(i), ratedPower(i), AEP); end

results表里的每一项都要在报告里给出计算方法和不确定度。比如风功率密度的不确定度主要来自风速测量误差和空气密度估算误差，一般给 ±10%。年发电量的不确定度更大，能到 ±15%。这些不确定度不是摆设，投资决策里的内部收益率计算要用到。机型对比那段代码里的CF是简化处理，实际应该用功率曲线和 Weibull 分布做卷积积分，但那个计算量大，初步筛查用容量系数就够了。

最后说一个我自己的习惯：每次跑完评估，我都会把原始数据、预处理后的数据、以及所有中间结果存成不同的文件，命名里带上日期和版本号。风资源评估经常要反复迭代，业主会提各种新要求，比如换个机型、改个轮毂高度、加个尾流方案。如果中间结果没存，每次都要从头跑一遍，浪费时间不说，还容易在重复计算中引入不一致。MATLAB 的save命令配合-v7.3参数能存大文件，timetable 也能直接存，恢复的时候用load就行。这个习惯看起来笨，但关键时刻能救命。希望帮到你。

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