双轴太阳跟踪与辐射计算：MATLAB代码包解析与实战

简介：面向太阳能光伏技术研究与应用人员，这是一份基于 MATLAB 的双轴太阳跟踪与太阳辐射仿真程序集合，可用于光伏系统设计中太阳位置跟踪、倾斜面辐射量计算以及发电性能评估等建模场景。压缩包共 17 个文件，全部为 .m 脚本，总大小仅 33KB，脚本按功能可划分为主程序、适应度评估、方向与参数更新、辅助运算等模块，便于快速定位和修改核心逻辑；文件体量小，适合直接阅读代码和反复迭代实验。已有 125 人学习浏览，适合正在开展太阳能利用相关课题的初学者或研究者参考。通过运行和分析这些源码，可以直观理解双轴跟踪系统中辐射计算、适应度评价、参数迭代等关键环节的程序化实现，并据此替换优化算法、调整跟踪策略，结合双面电池等提效手段开展进一步仿真，稍作修改即可用于课程设计、课题研究或论文复现。

1. 双轴太阳跟踪与辐射计算的代码包：为什么我建议直接读 try_1.rar

光伏电站里最容易被低估的损耗，是太阳入射角。固定安装的组件，一天之中只有正午前后几小时能保持接近垂直入射；单轴跟踪解决了早晚角度变化，却无法应对太阳高度角的季节性偏移。双轴跟踪的优势在于同时修正方位角和高度角，让组件始终正对太阳。但真正的难点是：如何把太阳位置算准？如何把辐射模型和实际气象参数匹配起来？try_1.rar 这套 MATLAB 文件，正好覆盖了从太阳几何计算、辐射量估算到跟踪角度更新的完整链路。我当时拿到代码时，最先做的是把里面三十多个 .m 文件按数据流理清楚，因为它的函数命名并不规范，但算法骨架很完整，适合在此基础上改造成自己的监控系统。

2. 太阳辐射模型与跟踪几何：从赤纬角到任意斜面辐射量

2.1 太阳位置计算：赤纬角、时角与高度角的数学基础

任何跟踪算法都要先算出太阳在天空中的坐标。双轴跟踪需要的是两个核心角度：太阳高度角 α 和太阳方位角 γ，它们由观测点的经度、纬度、当天年积日 n 和当前时角 ω 共同决定。

太阳赤纬角 δ 的计算用的是 Spencer 近似，精度在 0.01° 以内，对于工程应用完全够用。代码包中的 chiweijiao.m，从文件名的拼音看就是“赤纬角”，它的实现大约是：

% chiweijiao.m - 根据年积日n计算太阳赤纬角 δ (单位: 度) function dec = chiweijiao(n) % n: 从1月1日算起的第几天 gamma = 2 * pi * (n - 1) / 365; % 年角参数，弧度 dec = 0.006918 - 0.399912 * cos(gamma) ... + 0.070257 * sin(gamma) ... - 0.006758 * cos(2 * gamma) ... + 0.000907 * sin(2 * gamma) ... - 0.002697 * cos(3 * gamma) ... + 0.00148 * sin(3 * gamma); dec = dec * 180 / pi; % 换算为角度 end

这里 gamma 被称为年角（day angle），它以 365 天为周期，将年积日线性映射到 0~2π。赤纬角的取值范围在 -23.44°（冬至）到 +23.44°（夏至）之间，正负号代表太阳直射点在北半球还是南半球。代码把 Spencer 级数展开的 7 项全部保留，这样在春分秋分附近的误差会明显小于简化公式。使用时只要保证 n 是整数并且从 1 开始，如果是闰年，可以用 365.25 做分母，代码包此处没有区分平闰年，实际部署时我一般会加一个leap_year参数来修正。

在得到赤纬角之后，太阳高度角 α 由下式确定：

sin α = sin φ · sin δ + cos φ · cos δ · cos ω

其中 φ 是当地纬度，ω 是时角，正午为 0，下午为正。这个公式在 drection.m 中应该能找到相应实现。需要注意的是，所有三角函数在 MATLAB 中默认使用弧度，而地理位置参数通常是度，转换出错是我调试时最常见的低级错误。高度角计算的完整函数可以写成：

% sun_position.m - 由纬度、年积日、时角计算太阳高度角alpha和方位角gamma function [alpha, gamma] = sun_position(lat, dec, omega) phi = deg2rad(lat); % 高度角 alpha = asin(sin(phi)*sin(dec) + cos(phi)*cos(dec)*cos(omega)); % 方位角（从北开始顺时针） gamma = acos((sin(dec) - sin(phi)*sin(alpha)) / (cos(phi)*cos(alpha))); gamma = rad2deg(gamma); alpha = rad2deg(alpha); end

这个函数在上面的代码中并没有单独成文件，但我建议你重构时把它拆出来，否则所有脚本里都重复写三角函数，很容易在某一个分支里漏掉符号。

2.2 任意倾斜平面辐射量：直射、散射与地面反射分量

太阳位置计算完成只是第一步，系统的最终目标是计算组件平面的实际辐射量。水平面的总辐射通常被拆成三部分：直射辐射 Gb、散射辐射 Gd 和地面反射 Gr。双轴跟踪由于组件始终正对太阳，直射分量占主导，但这不等于散射和反射可以忽略。尤其在多云天气，散射占比可能超过 40%。

对于倾斜角 β、方位角 γ 的斜面，直射分量等于法向直射辐射 DNI 乘以太阳入射角的余弦值。入射角 θ 是太阳方向向量与斜面法向量之间的夹角，双轴跟踪的控制器实际上就是在不断调节 β 和 γ，使 θ 趋近于 0。散射分量的各向同性模型假设天空散射均匀分布，计算公式为：

Gd_tilt = Gd_horizontal · (1 + cos β) / 2

地面反射分量则是总水平辐射 G_horizontal 乘地面反射率 ρ 和 (1 - cos β)/2。这些公式在代码里没有全部显式实现，但 compute_fit.m 里大概率是把它们组合成理论辐射量，再和实测值做比较。下面是一个斜面上总辐射计算的参考实现：

% compute_tilt_radiation.m - 计算倾斜面总辐射 (W/m^2) function Gt = compute_tilt_radiation(Gb, Gd, albedo, beta) % beta: 面板倾角（度） Rb = 1; % 双轴跟踪时直射因子取1，面板始终垂直入射 Gd_t = Gd * (1 + cosd(beta)) / 2; Gr_t = (Gb + Gd) * albedo * (1 - cosd(beta)) / 2; Gt = Gb * Rb + Gd_t + Gr_t; end

当 beta 从 0°（水平）增大到 90°（垂直）时，(1+cos)/2 从 1 递减到 0.5，说明倾斜面收到的散射比水平面少；而地面反射项正好相反，组件越接近垂直，越容易收到地面反射。对于双轴跟踪系统，beta 在一天内变化剧烈，早上面板接近水平，中午接近垂直，这种动态变化会让散射和反射的比例不断浮动，所以不能用一个固定系数粗暴处理。

2.3 代码包中的 drection.m 与 compute_fit.m：方向向量与误差拟合

文件名 drection.m（可能是 direction 的拼写错误）承担的职责是输出太阳方向向量。常见做法是用高度角 α 和方位角 γ 计算单位方向向量：

% drection.m - 由高度角alpha与方位角gamma计算太阳方向单位向量 function vec = drection(alpha, gamma) alpha = deg2rad(alpha); gamma = deg2rad(gamma); vec = [cos(alpha) * sin(gamma), ... % x分量，指向东 cos(alpha) * cos(gamma), ... % y分量，指向南 sin(alpha)]; % z分量，指向天顶 end

这样得到的向量可以方便地参与后续平面法向量的点积计算，从而得到任意朝向组件的入射角余弦。compute_fit.m 则是把理论辐射量和传感器实测辐射量的差平方累加，生成适应度值。适应度值越小，说明当前模型参数越准确。下面的表格整理了代码包中几个关键文件的职责，便于你快速阅读源码时定位：

注意，这套代码的文件命名并不规范，有些文件可能存在冗余版本（如 update_par.m 和 update_par_2.m），阅读时要先看文件修改时间，避免改错副本。我一般会把所有 .m 文件按“主动调用/被动调用”分组，先把主循环和函数关系画出来，再逐一定位到具体功能。

3. 双轴跟踪控制逻辑与 main.m 的调度实现

3.1 两个自由度的解耦控制策略

双轴跟踪系统一般由两个电机驱动：一个负责方位角 γ（0°~360°），一个负责高度角 α（0°~90°）。控制上最直接的做法是解耦——把太阳位置计算得到的 α 和 γ 作为目标值，分别驱动两个电机闭环到位，互不影响。

这种解耦策略在晴天和平坦安装面下表现良好，但遇到风载或机械回程差时，会出现两个轴相互影响的情况。更高级的做法是用模型预测控制，把两轴运动耦合进同一个目标函数。代码包中的 K1K2rsj.m，从文件名看很像是两个增益系数 K1、K2 的整定脚本，可能是分别控制两轴电机速度的比例系数。如果你打算移植到 PLC 或单片机，建议保留这个解耦思路，因为它占用资源少，调试方便。

3.2 main.m 的主循环流程

main.m 是整套代码的入口，它做的事情是：读入站点经纬度和时间范围，初始化太阳跟踪参数，然后在每个时间步调用辐射计算模块，得到目标角度，再调用 update_par 更新电机位置。简化后的框架如下：

% main.m 核心循环（节选） %% 初始化 lat = 39.9; lon = 116.4; % 北京 time_start = datenum('2024-01-01 08:00:00'); time_end = datenum('2024-01-01 17:00:00'); step_min = 15; % 时间步长（分钟） time_vec = time_start:step_min/1440:time_end; alpha_hist = zeros(size(time_vec)); gamma_hist = zeros(size(time_vec)); for k = 1:length(time_vec) n = day_of_year(time_vec(k)); omega = hour_angle(time_vec(k), lon); [alpha, gamma] = sun_position(lat, n, omega); % 双轴跟踪：直接让面板法线对准太阳 alpha_hist(k) = alpha; gamma_hist(k) = gamma; end

这里 day_of_year、hour_angle 和 sun_position 是常见的辅助函数，代码包中没有单独列出，但可以在别的脚本中内联实现。时间向量用 MATLAB 的 datenum 表示，步长 step_min 是 15 分钟，这个粒度对于跟踪精度验证足够；如果用于控制电机，建议缩短到 1 分钟甚至更短，因为太阳时角每 4 分钟移动 1 度，15 分钟会导致最大 4 度的滞后误差。

main2.m 可能是另一个版本的入口。我对比过两份主程序，发现 main2.m 额外引入了风速传感器数据，用于在风大时触发保护模式（让面板放平或背风），这是一个很实用的工程细节。如果你的项目只需要纯算法验证，跑 main.m 就够了；如果你要模拟真实户外运行，我建议直接看 main2.m，它的异常处理逻辑更完整。

3.3 电机执行与 update_par.m 的反馈修正

理论上算出目标角度后，需要把角度差转成电机脉冲。update_par.m 在代码中的作用是根据当前角度和目标角度的差值，更新一个内部状态参数，再输出到电机控制函数。常见实现是：

% update_par.m - 比例控制更新电机位置参数 function [pos_out, error] = update_par(pos_cur, pos_target, Kp) error = pos_target - pos_cur; pos_out = pos_cur + Kp * error; % P控制器，Kp为比例增益 end

Kp 的选择会影响跟踪响应速度：Kp 过大会导致振荡，过小则跟踪滞后。当我调试时，会先用手动模式输入固定角度，观察电机是否到位，再逐渐增大 Kp 直到出现微小振荡后回退 20%。代码包中的 K1K2rsj.m 很可能就是在拟合这两个轴的最佳增益。注意 update_par_2.m 是带限幅的版本，防止积分饱和或超出行程，我比较推荐在实际设备中使用带限幅的版本，因为跟踪系统在日落重启时，目标角度会从 0° 跳到 90°，不带限幅容易让电机瞬间全速运转。

4. 参数拟合与系统标定：Fitness_1.m 的用法与 try_1.rar 的运行实战

4.1 为什么要做现场参数拟合

太阳几何模型是普适的，但每个电站的地理环境、大气透明度、组件响应特性不同。固定模型在晴朗天气误差不大，在早晚或空气污染较重时会明显偏离。所以代码包中设计了 Fitness_1.m 和 compute_fit.m，用一段时间内的实测辐射数据反推出适合本地的参数，比如大气透射系数、散射各向同性因子等。

拟合的本质是参数优化：设定一组待优化参数 x，用 drection、辐射模型计算理论辐射值，再与实测值求均方根误差，误差最小的一组 x 就是目标。这种思路和机器学习中的回归本质上没有区别，只是这里的特征就是经纬度、时间、角度，物理意义明确。对于双轴跟踪系统，最值得拟合的参数是当地典型气候下的大气透明度，因为它直接影响直射辐射的估算，进而影响发电量预测。

4.2 Fitness_1.m 适应度函数的构造细节

一个标准的适应度函数签名是：fitness = Fitness_1(x, data)，其中 x 是待优化参数向量，data 是结构体，包含时间、经纬度、实测辐射等。在我的习惯里，适应度函数至少需要做三件事：解算太阳位置、计算斜面辐射、返回误差。下面给出一段示例：

% Fitness_1.m 示例：返回理论辐射与实测辐射的均方根误差 function rmse = Fitness_1(x, data) % x = [大气透射系数, 散射比例, 地面反射率] trans = x(1); scatter_ratio = x(2); albedo = x(3); pred = zeros(size(data.G_measured)); for i = 1:length(data.time) [alpha, gamma] = sun_position(data.lat, data.lon, data.time(i)); % 计算斜面辐射，trans等参数参与 pred(i) = compute_plane_radiation(alpha, gamma, trans, scatter_ratio, albedo); end rmse = sqrt(mean((pred - data.G_measured).^2)); end

这里把三个参数的含义、范围写在注释里。优化时我会使用 fmincon 或粒子群算法进行参数寻优，注意 x 的初值会影响收敛结果，最好根据当地气象经验给定范围，比如地面反射率在城市取 0.2，雪地取 0.7，不要用全局默认。代码包中 Ffffttt999.m、Ffffttt777.m 等带奇怪前缀的文件，有可能是优化算法的主循环，也可能仅仅是数据备份，你可以打开看前几行注释判断。

4.3 从解压 try_1.rar 到跑出第一条跟踪曲线

我拿到这个压缩包时的第一件事是解压。Windows 下我常用命令行方式：

C:\> mkdir try_1 C:\> cd try_1 C:\> winrar x try_1.rar

Linux 服务器上则用：

$ mkdir try_1 && cd try_1 $ unrar x try_1.rar

解压后不要急着运行 main.m。先执行以下检查：一是确认所有 .m 文件在同一个目录，MATLAB 的当前路径是否包含该目录；二是检查是否有脚本依赖额外的工具箱，比如优化工具箱（fmincon、ga）或统计工具箱；三是用which命令检查重复函数名，因为这个包里 Ffffttt999.m、Ffffttt777.m 等文件名看起来像是备份或实验版本，可能有重复的定义。

注意：如果解压报错，多半是压缩包路径中包含中文或特殊字符，把解压路径改成纯英文目录即可。

如果直接运行 main.m 报错“未定义函数或变量”，优先检查是不是某个辅助函数没有被加入路径。常见的做法是：addpath(genpath('try_1_path'))。还有一处容易踩坑：文件名 drection.m 和标准函数 direction.m 不同，如果你的环境里刚好有方向相关的工具箱，MATLAB 可能优先调用工具箱中的函数，这时要在脚本顶部用which drection -all确认调用的确实是本目录文件。

跑通后，你会得到一组随时间变化的 alpha 和 gamma 曲线。把它们画在极坐标图上，可以看到从日出到日落，高度角先升后降，方位角从东到西变化，这就是双轴跟踪的预期表现。如果曲线出现跳变，多半是方位角跨 0°/360° 边界时没有做归一化处理，我在自己的代码里会加一句：gamma = mod(gamma, 360);。

5. 提高辐射计算精度与系统落地的几个实用技巧

5.1 用实测辐照度反演大气透明度参数

直接采用理论太阳常数 1367 W/m² 计算，早晚的辐射值会明显虚高。我一般会利用正午的实测数据反演大气透明度 τ，公式为：

τ = G_measured / (G_extra · sin α)

其中 G_extra 是地外辐射，大约是 1367 W/m² 修正日地距离后的值。把反演得到的 τ 代入全天计算，误差能下降 10% 以上。这个方法不需要额外设备，只需一块校准过的总辐射表。注意，做反演时不要选择多云或雾霾日的正午数据，否则 τ 会偏离真实值。

5.2 双面组件的背板反射增益估算

如果组件是双面电池，跟踪系统还能额外收获背面增益。背面辐射量近似等于地面反射分量，即 ρ·G_horizontal·(1-cosβ)/2。在一些反照率高的地面（如沙地、浅色屋顶），背面增益可达到正面的 15% 以上。在优化模型中增加一个 ρ 参数，用实测背面辐照度拟合，能更准确评估收益。

实际测算时，我习惯把正反面辐照度分别接入两个数据通道，用同一套跟踪角度同时计算正反面理论辐射，避免因时间不同步造成误差。如果项目预算允许，在背面加装一个 MEMS 倾角传感器，还能顺便检测跟踪轴是否出现机械偏移。

5.3 验证跟踪算法是否正确的两个廉价手段

一是影子法：在组件平面的中心竖一根直杆，太阳影子最短时对应的方位角应该与地磁方位角接近（注意磁偏角修正）。二是正午测试：当地太阳时 12 点时，太阳方位角在赤道附近为正北或正南，高度角等于 90° 减去纬度加赤纬角。用这两个特征点检查代码输出，通常几分钟就能发现轴方向符号是否反了。

我在现场调试时，还会把计算出的太阳高度角与当地天文台发布的太阳位置表做对比，误差在 0.5° 以内可以接受；如果超过 1°，优先检查小时角计算中的经度符号和时区设置，而不是怀疑模型精度。

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