简介:这是一份基于MATLAB实现的视日轨迹跟踪算法仿真资源,主要面向新能源、自动化等专业的本科与硕士生,以及从事光伏支架控制、太阳能利用等相关研究的工程技术人员。资源围绕太阳位置计算与视日运动轨迹模拟展开,能够为课程设计、毕业设计或科研课题提供可直接复现的算法参考与实验起点。压缩包内共4个文件,核心是一个MATLAB程序文件,可直接运行仿真流程;另有两个png格式的效果展示图,便于直观核对输出结果;一个txt格式说明文档,对版本与运行环境等作了简要交代。整包大小约459KB,轻量、紧凑,下载后可快速部署到本机环境。目前已有158人浏览学习此资源,使用时可结合源码、说明和截图快速理清视日轨迹跟踪的建模思路与实现步骤,并在此基础上进行参数修改、功能扩展或算法改进,对于深入理解光伏视日运动规律和开发相关控制策略都很有帮助。
1. 视日轨迹跟踪算法在Matlab里的工程化落地
视日轨迹跟踪算法解决的问题很直接:给定所在地经纬度、当前日期与时刻,算出太阳此刻在哪。与常见的光敏传感器闭环追踪方案相比,这种开环前馈计算不依赖天气和硬件反馈,在光伏双轴支架、塔式光热定日镜场、建筑遮阳与采光分析中都有直接用途。这个Matlab工程包内置了完整的main.m运行脚本,基于2014a及以上版本均可直接运行,输出太阳高度角与方位角的全天曲线。我拿到源码后的第一反应是核对它的时间基准与坐标系定义,因为这两处是后续所有计算正确与否的分水岭,也是论文里最容易被评委追问的细节。
2. 赤纬角、时角与真太阳时:时间基准先立住
视日轨迹算法的本质,是把太阳在黄道上的运动投影到观测者所在的水平坐标系。这一步要先后经过赤道坐标系和地平坐标系,中间涉及三个关键输入:赤纬角(太阳直射点纬度)、时角(太阳相对当地子午线的角位移)、以及一个干净的时间基准。任何一项有偏差,最终的高度角和方位角都会整体偏移,所以这一章先把天文基础打牢。
2.1 赤纬角的两个模型:Cooper与Spencer的取舍
赤纬角的定义是地心天球坐标系中太阳中心相对天赤道的角距离,取值范围在-23.45°到+23.45°之间。工程中最常用的简化模型是Cooper方程:
δ = 23.45 × sin(2π/365 × (284 + N))
其中N为年积日,1月1日为1。这个模型的推导基础是把地球公转近似为椭圆与自转轴的固定夹角,误差来源是忽略黄经与真近点角的差异,全年最大偏差约0.35°,出现在春秋分附近。
如果用来驱动双轴跟踪支架,Cooper模型的偏差换算成跟踪角度误差可能让光斑偏移一个塔式定日镜的镜面宽度,此时需要用Spencer模型做更高精度的逼近:
% 年积日计算 N = datenum(year, month, day) - datenum(year, 1, 1) + 1; % Spencer赤纬角(弧度),B为回归年比例系数 B = 2*pi*(N-1)/365; delta = 0.006918 - 0.399912*cos(B) + 0.070257*sin(B) ... - 0.006758*cos(2*B) + 0.000907*sin(2*B) ... - 0.002697*cos(3*B) + 0.001480*sin(3*B);这段代码把年积日换算成弧度相位B再代入7项傅里叶级数。Spencer模型的精度在0.01°量级,相比Cooper模型提升了一个数量级,代价只是多5次三角函数运算,对现代处理器来说可以忽略。工程包里选择Spencer还是Cooper,可以从仿真时间跨度来推断:单日曲线用哪个都行,全年逐时扫描就必须用Spencer,否则累加误差会在年发电量估算中放大。两个模型的适用边界可以这样对照:
| 模型 | 最大赤纬误差 | 适用场景 | 计算成本 | | Cooper | ±0.35° | 快速估算、教学演示 | 1次sin | | Spencer | ±0.02° | 跟踪控制、全年仿真 | 5次sin/cos |
提示:赤纬角公式里的365与366(闰年)差异很小,实际影响低于0.01°,多数工程实现不区分平闰年,直接用365即可。
2.2 时差方程与经度修正:从钟表时间到真太阳时
视日轨迹算法要求的时间是"真太阳时",日常使用的钟表时间是"平太阳时"。两者差值由两个因素叠加:一是地球公转速率不均匀导致的时差方程(Equation of Time),二是观测点经度与所在时区中心经度的偏移。工程里常见的错误是把钟表时间直接当太阳时使用,这样在北京(东经116.4°)时,正午会出现约14分钟的偏移,换算成时角就是3.6°,足以让高度角误差达到0.5°以上。
% 时差方程(分钟),公式与Spencer同源 EoT = 229.18 * (0.000075 + 0.001868*cos(B) - 0.032077*sin(B) ... - 0.014615*cos(2*B) - 0.040849*sin(2*B)); % 真太阳时(小数小时) solar_time = hour + minute/60 + ((lon - 15*tz)*4 + EoT) / 60; % 时角(度),每1小时对应15度 H_deg = 15 * (solar_time - 12);代码里lon代入的是东经正值、西经负值,tz代入的是UTC偏移小时数。经度修正项(lon - 15×tz)×4的含义是:经度每偏离时区中央经线1度,地方时间差4分钟;东侧时间领先,所以为正值。EoT的单位是分钟,除以60换成小时后再叠加。这个修正完成后,H_deg在日出时约为负值,正午时为0,日落时转为正值,方向约定与后续方位角公式严格配套。
2.3 高度角与方位角:球面三角形的两种解法
有了赤纬角δ和时角H,高度角α由球面余弦定理直接给出。这里需要把经纬度和角度全部统一为弧度再计算,否则Matlab的sin/cos函数会因为单位混用给出完全错误的结果。
lat_r = lat * pi/180; % 纬度转弧度 dec_r = delta; % delta已是弧度 H_r = H_deg * pi/180; % 时角转弧度 % 高度角(度) sin_alt = sin(lat_r)*sin(dec_r) + cos(lat_r)*cos(dec_r)*cos(H_r); alt = asin(sin_alt) * 180/pi; % 方位角(度,从正北顺时针) az_r = atan2(sin(H_r), cos(H_r)*sin(lat_r) - tan(dec_r)*cos(lat_r)); az = mod(az_r*180/pi + 180, 360);方位角的求解有两个常见路线。第一条路是通过余弦定理求出相对南方的角度,再用上午/下午判断东西侧,需要额外的逻辑分支;上面代码用的是atan2单一个表达式,分母的综合相位天然对应了太阳在天空的象限,省去判断分支。atan2的优势是值域覆盖-pi到pi,不会出现acos那种在0°、180°附近对微小扰动过于敏感的情况。mod(az+180,360)把参考零点从南方移到正北,得到的是北偏东为正的标准气象方位角。
这套函数封装完成之后,剩下的就是main.m里的逐时刻扫描与可视化。我在第3章把它展开到完整可运行的脚本。
3. main.m 实现拆解:从函数封装到全天曲线
工程包的main.m脚本承担三个任务:设定站点与日期参数、逐时刻调用太阳位置函数、绘制高度角与方位角随时间的曲线。下面给出的实现保留了原工程的核心流程,并补齐了日出日落界面的处理,是实际工作中最常用的写法。
3.1 输入参数区与函数封装
站点参数放在脚本顶部集中管理,改一个城市只需要改四个数字。这里用datenum统一处理日期,跨月、跨年循环时不用手工计算每月天数,也不容易漏闰年。
clc; clear; close all; % 站点:北京 lat = 39.9042; % 纬度(°),北纬为正 lon = 116.4074; % 经度(°),东经为正 tz = 8; % 时区(小时),UTC+8 % 日期与扫描步长 year = 2024; month = 6; day = 21; % 夏至日 step_min = 10; % 扫描步长(分钟) % 生成时间序列 t_start = datenum(year, month, day, 0, 0, 0); t_stop = datenum(year, month, day, 23, 59, 59); t_vec = t_start : step_min/(24*60) : t_stop;这里的datenum把时间编码为连续浮点数,1代表1天。step_min/(24×60)把分钟换算成天数增量,t_vec就得到从当天0点到23:59的全部时间戳。用浮点时间的好处是调用函数时直接解析时分秒,不需要单独维护每个小时刻的索引。
3.2 逐时刻扫描与数组预分配
太阳位置函数需要独立的year、month、day、hour、minute参数,所以循环里要从t_vec反解这些量。Matlab的datevec可以一次完成:
n = length(t_vec); alt = zeros(n, 1); az = zeros(n, 1); for i = 1:n [y, mo, d, h, mi, ~] = datevec(t_vec(i)); [alt(i), az(i)] = solar_position(lat, lon, y, mo, d, h, mi, tz); enddatevec返回的第6个元素是秒,用波浪线丢弃。数组alt、az预先用zeros分配,避免循环里动态扩容导致的性能下降——这在扫描全年、步长1分钟的场景下(约52万次循环)能省下10倍以上的时间。solar_position函数内部按第2章的公式逐步求值,其输入输出约定与绘图脚本完全解耦,单独维护函数文件的好处是后续切换Cooper/Spencer模型时只动一处。
function [alt_deg, az_deg] = solar_position(lat_deg, lon_deg, ... year, month, day, hour, minute, tz) % 输入均为原始单位,输出高度角、方位角(°,从北顺时针) N = datenum(year, month, day) - datenum(year, 1, 1) + 1; B = 2*pi*(N-1)/365; delta = 0.006918 - 0.399912*cos(B) + 0.070257*sin(B) ... - 0.006758*cos(2*B) + 0.000907*sin(2*B) ... - 0.002697*cos(3*B) + 0.001480*sin(3*B); EoT = 229.18*(0.000075 + 0.001868*cos(B) - 0.032077*sin(B) ... - 0.014615*cos(2*B) - 0.040849*sin(2*B)); solar_time = hour + minute/60 + ((lon_deg - 15*tz)*4 + EoT)/60; H_deg = 15 * (solar_time - 12); lat_r = lat_deg * pi/180; dec_r = delta; H_r = H_deg * pi/180; sin_alt = sin(lat_r)*sin(dec_r) + cos(lat_r)*cos(dec_r)*cos(H_r); alt_deg = asin(sin_alt) * 180/pi; az_r = atan2(sin(H_r), cos(H_r)*sin(lat_r) - tan(dec_r)*cos(lat_r)); az_deg = mod(az_r*180/pi + 180, 360); end这个函数里,纬度和经度、时区作为普通数值参数传入,便于后续用数组或表格驱动批量仿真。一个需要注意的边界是:当高度角为负(太阳在地平线下)时,方位角数值仍有意义,但水平坐标系下不含大气层的几何计算会给出一个"虚拟太阳"方向,画图时通常用逻辑索引过滤掉。
3.3 绘图与时间轴格式化
主程序的绘图部分需要同时呈现高度角与方位角两条曲线,并把横轴格式化为"时刻"而非浮点日期数:
figure('Color','w','Position',[100 100 900 400]); subplot(1,2,1); plot(t_vec, alt, 'b-', 'LineWidth', 1.2); hold on; plot(t_vec, zeros(size(t_vec)), 'k--'); xlabel('时刻'); ylabel('太阳高度角(°)'); title('北京 2024-06-21 高度角'); datetick('x', 'HH:MM'); xlim([t_start t_stop]); ylim([-10 90]); grid on; subplot(1,2,2); plot(t_vec, az, 'r-', 'LineWidth', 1.2); xlabel('时刻'); ylabel('太阳方位角(°)'); title('北京 2024-06-21 方位角'); datetick('x', 'HH:MM'); xlim([t_start t_stop]); ylim([0 360]); grid on;datetick('x','HH:MM')是Matlab里处理时间轴的关键函数,它把datenum浮点数横轴重新标记为人类可读的时分格式。xlim与t_start、t_stop绑定可以防止datetick在数据范围外多画出多余刻度。黑色虚线标记地平线,高度角曲线与虚线的两个交点就是日出与日落时刻,从图上可以直接读出来。如果实测中发现交点位置与天文年历不同,问题基本都出在第2章的时间修正项上。
4. 仿真结果判读与精度校验:误差藏在这些地方
模型跑出曲线只是第一步,判断曲线对不对、误差来自哪里,是工程落地前必须做的功课。这一章给出三条校验路径:正午峰值对照、全天形态比对、边界时刻验证。
4.1 正午高度角的特征检验
太阳高度角在真太阳时正午达到峰值。对北半球中纬度的观测者,正午高度角的理论值由90°-纬度+赤纬角直接给出。以北京(39.9042°N)2024年夏至日为例:
| 节气 | 日期(近似) | 赤纬角(°) | 正午高度角理论值(°) | | 春分 | 3月20日 | 0 | 50.1 | | 夏至 | 6月21日 | +23.45 | 73.6 | | 秋分 | 9月23日 | 0 | 50.1 | | 冬至 | 12月22日 | -23.45 | 26.6 |
跑完main.m后,把夏至日曲线峰值与73.6°对比,如果偏差超过0.5°,优先检查时区tz与经度lon是否输错。这里有一类常见误用:有人把经度代入正值但忘了东八区的tz=8,或者把时区当成0(UTC)而经度还是东经116°,真太阳时直接偏晚8小时,高度角曲线整体从正午向右平移,峰值不再出现在12:00附近而出现在19:00附近,这是最容易肉眼识别的故障。
4.2 误差来源的量化拆解
即使参数全部正确,计算值与理想几何模型之间仍然存在系统性偏差。这套算法产出的是"几何太阳位置",不是"视太阳位置",两者的主要差异来自大气折射。高度角大于10°时,折射造成的抬高量小于0.1°,可以忽略;但在日出日落附近,折射会造成约0.5°以上的表观抬升,工程上判断昼夜分界时一般把高度角修正为-0.833°再判零,也就是把折射补偿近似固定为0.833°。
时差方程与赤纬角的模型误差也会叠加入最终结果。Spencer模型里EoT的误差约在±0.5分钟以内,对应时角误差约0.125°;赤纬角0.02°的误差乘以cos(φ)的影响系数后,最终高度角误差在0.02°量级。三者叠加后的综合误差不超过0.2°,满足绝大多数光伏跟踪支架0.5°以内的转角控制需求。如果发现峰值偏差超过这个范围,我一般先用行星历表抽查某个特定时刻,再回过来核对公式里的经纬度符号。
4.3 日出日落时刻与负高度角过滤
日出日落时间的解析解可以直接从高度角方程反推,更简单的做法是在已有曲线上做阈值检测。对于10分钟步长,线性插值就能把边界时刻误差控制在几分钟内;如果步长是60分钟,日出日落附近曲线曲率高,线性插值会产生较大偏差,这时就要用fzero做精确求根。
% 精确计算日出时刻:在高度角定义中嵌入折射补偿-0.833° refraction_offset = -0.833 * pi/180; % 定义高度角函数,输入为儒略日,输出为高度角-补偿值 f = @(t) sin_alt_function(t, lat, lon, tz) - sin(refraction_offset); % 日出在0~12点之间搜索 sunrise = fzero(f, [t_start, t_start + 0.5]); % 日落在12~24点之间搜索 sunset = fzero(f, [t_start + 0.5, t_stop]);fzero需要函数在区间两端异号,所以日出区间必须选在高度角从负到正的范围内,日落区间选在正到负。这里的sin_alt_function复用solar_position函数,只需把datenum分解出时分秒再调用。对比几何法(按0°判零)与折射修正法(按-0.833°判零),后者得到的白昼时长在春秋季平均长4-6分钟,这在实际光热电站的运行策略里会直接影响早晨吸热器启动时刻的选择。
5. 进阶:从视日轨迹到双轴跟踪与全年辐照评估
视日轨迹计算的最终价值在于驱动机构和评估产能。这一章给出两个最直接的进阶用法。
5.1 双轴支架转角解算
对双轴跟踪系统,最常见的结构是高度角-方位角型(AZ-EL)和俯仰-滚转型(Tilt-Roll)。前者直接使用本算法输出的alt和az作为两个转轴的指令角;后者需要先转换到天顶角θz = 90° - alt,再把天顶角分解到支架的倾斜轴和旋转轴上。需要注意,视日轨迹跟踪是纯几何位置跟踪,没有考虑阴天散射辐照的优势方向,因此云量较大地区应考虑低成本时控策略代替全时跟踪,以节省驱动能耗。
% AZ-EL双轴转角指令 theta_z = 90 - alt; % 天顶角 azimuth_cmd = az; % 旋转轴指令 % 检查角度变化率,防止跳变 delta_az = [0; diff(az)]; delta_az = mod(delta_az + 180, 360) - 180; % 折回处理delta_az的处理是为了应对方位角从355°跨到5°时差分出现的大幅跳变。直接用diff会得到一个约-350°的错误变化率,折回处理后变成+10°,这个值才是伺服电机的真实角速度需求。省掉这一步,在仿真报告里会出现瞬时上千度每秒的转速尖刺,实际电机选型时会误判。
5.2 全年逐时仿真与年发电量粗估
把main.m的单日计算扩展到全年,可以得到365×24的位置矩阵。结合一个简单的晴空辐照模型(例如Hottel模型),就能在没装传感器的情况下估算双轴跟踪相对固定角度安装的发电增益。
% 全年逐时太阳位置矩阵 days = 1:365; hours_of_day = 0:23; alt_matrix = zeros(length(days), 24); az_matrix = zeros(length(days), 24); for d = 1:length(days) [y, mo, da] = datevec(datenum(year, 1, 1) + d - 1); for h = 1:24 [alt_matrix(d,h), az_matrix(d,h)] = solar_position(... lat, lon, y, mo, da, h, 0, tz); end end % 过滤掉高度角小于5°的低光照时段 valid = alt_matrix > 5;这个矩阵的生产成本在Matlab里约在0.5秒以下。后续处理时,把晴空直射辐照IDNI乘以cos(入射角),再按valid掩码累加,就能得到月均发电量相对值。相比用PVsyst做全年仿真,这套自研流程胜在可控和透明——每个系数、每个滤波条件和辐照模型都能在论文或者技术报告里写清楚,对于硕士课题和预研项目已经足够。矩阵里的高度角变化率也可以顺带统计出来,用来校核支架电机的最大角速度需求,这一步在招标技术参数表里经常作为硬性指标出现。
本文还有配套的精品资源,点击获取