企业官网建设流程全解析

简介：面向导航、遥感与天文应用场景的MATLAB太阳位置计算程序包，解决给定经纬度下太阳方位角与高度角的精确求解问题，也适用于偏振导航研究中由太阳方向推算偏振角等参数。压缩包仅6KB，包含9个文件，以8个.m脚本为主，另附1个说明文档。脚本覆盖从经纬度与时间输入、天文公式运算到地平坐标转换、最终角度输出的完整流程，并涉及deg2rad、atan2等典型函数用法；说明文档可辅助理解各脚本调用关系。目前已有463人学习下载。借助这份代码，读者既能直接获得可运行的太阳方位角与高度角计算工具，也能通过SolarAngle、skew_symmetric、FaiToSouth等模块了解偏振角建模、向量叉乘等关键细节，适合具备一定MATLAB基础、希望深入太阳位置算法与偏振导航原理的开发者参考学习。

1. 这套太阳方位程序为什么值得拆

做偏振光导航或者天文定位的人，手里大部分太阳位置算法都是C写的，换到MATLAB环境经常要重写一遍坐标转换。而这个名为“太阳方位matlab程序.rar”的压缩包，恰好把一套完整的太阳方位角、太阳高度角计算流程用MATLAB函数封装好了。压缩包里既有SolarAngle.m这样的主函数，也有skew_symmetric.m、FaiToSouth.m这类辅助工具，一看就是从实际项目里拆出来的代码，不是教学示例那种只算一个公式的玩具。对于需要把太阳位置写进仿真链路、或者做偏振角解算的工程师来说，这套代码的价值在于它帮你省掉了查天文年历和调试坐标系的重复劳动。本文就从天文算法、程序结构、偏振角应用和工程排错四个角度，把这份资源里值得移植的部分拆开讲透。

2. 太阳方位角和高度角的底层天文算法

2.1 你需要先理解的三个坐标系

很多人拿到太阳位置程序第一反应是去找公式，但公式里的变量代表什么才是真正容易出错的地方。计算太阳方位角时涉及三套坐标系：赤道坐标系（赤经、赤纬）、时角坐标系（时角、赤纬）和地平坐标系（方位角、高度角）。程序内部做的事情本质上就是：根据时间和地点求出太阳在赤道坐标系的位置，再转到时角坐标系，最后投影到地平坐标系。

在三套坐标系里，时角坐标系是最容易被忽略的。时角HA的定义是太阳所在子午圈与当地子午圈之间的夹角，以正南为0，向西为正，每小时对应15度。而经纬度输入longitude在整个计算里最大的作用，就是把UTC时间换算成当地太阳时。太阳过当地子午圈的时刻不是12点整，而是12:00 - longitude/15，加上均时差修正后才是真正的太阳正午。这一项不校正，方位角误差会超过2度，在偏振导航里这个误差完全不能接受。

2.2 太阳赤纬的近似模型足够用

太阳赤纬delta随时间变化，程序包里的SolarAngle.m大概率用的是以下两种模型之一。第一种是Cooper近似公式：

% 输入：一年中的第几天 day_of_year % 输出：太阳赤纬（弧度） delta = 23.45 * sin(2 * pi * (284 + day_of_year) / 365) * pi / 180;

这个公式精度大约在1度以内，适合算法验证和教学。第二种是Spencer级数展开，精度可以做到0.01度级别，工程上够用：

% day_of_year 为一年中第几天，doy 为对应的角参数 doy = 2 * pi * (day_of_year - 1) / 365; delta = 0.006918 - 0.399912 * cos(doy) + 0.070257 * sin(doy) ... - 0.006758 * cos(2 * doy) + 0.000907 * sin(2 * doy) ... - 0.002697 * cos(3 * doy) + 0.00148 * sin(3 * doy);

用Cooper公式做初值、再用Spencer做精算，是多档精度程序里的常见做法。判断一个太阳位置程序靠不靠谱，先看它用的是几阶展开，再看有没有均时差修正。如果两样都没有，那这个程序的精度只有粗匹配的水平。

2.3 方位角计算里的atan2陷阱

有了赤纬delta、当地纬度phi和时角HA，高度角alt和方位角az可以直接用球面三角公式求。但方位角求法有个分支陷阱：

% 输入：纬度 phi、赤纬 delta、时角 HA（均为弧度） % 输出：方位角 az（弧度，从北顺时针）和高度角 alt sin_alt = sin(phi) * sin(delta) + cos(phi) * cos(delta) * cos(HA); alt = asin(sin_alt); cos_az = (sin(delta) - sin(phi) * sin_alt) / (cos(phi) * cos(alt)); az = acos(cos_az); % 如果太阳在正南偏西，acos给出的角度才正确 if sin(HA) > 0 az = 2 * pi - az; end

这个if sin(HA) > 0判断是算对方位角的关键。acos函数的返回值范围是[0, pi]，对应从北经东到南的半圈，但下午太阳在西南方向，方位角应该在[pi, 2*pi]区间。不用atan2修正的话，下午的方位角会全部镜像到东南方向。更稳健的做法是直接用atan2一步到位：

az = atan2(sin(HA), cos(HA) * sin(phi) - tan(delta) * cos(phi)) + pi;

这个写法把所有象限判断交给了atan2，不会再出现镜子方向问题。实际工程里我用第二种写法比较多，因为不需要额外判断逻辑。

3. 拆解压缩包里的MATLAB程序结构

3.1 函数文件各司其职

压缩包里的文件命名很有规律，基本可以推断出每个文件的职责。SolarAngle.m和SolarAngle1.m一主一辅，区别大概率在输入参数格式上；xigongda.m和xigongda11.m有成对出现的嫌疑，可能是不同版本的计算主流程；skew_symmetric.m是计算反对称矩阵的工具函数，这在偏振矢量运算里相当常见；FaiToSouth.m和FaifromSun.m的Fai是角度变量，ToSouth和fromSun暗示这两个函数在做一个方向基准转换。

把这些函数串起来看，整套程序的工作流程应该是这样：xigongda.m作为主入口调用SolarAngle算出太阳的位置参数，再用FaiToSouth.m把天球坐标系里得到的太阳位置转到载体坐标系，最后通过skew_symmetric.m参与偏振矢量的叉乘运算。这里的Fai很可能就是偏振角（Polarization Angle），也就是E矢量振动方向与参考方向的夹角。

3.2 标准输入输出格式

写MATLAB程序最怕函数输入隐式依赖全局变量。判断这套代码能不能直接拿去用，建议先打开SolarAngle.m看看函数声明是下面哪种风格。

风格一（推荐，纯参数传递式）：

function [azimuth, elevation] = SolarAngle(lat, lon, year, month, day, hour, minute, second)

风格二（隐式全局变量式）：

function [azimuth, elevation] = SolarAngle() % 内部依赖 global lat lon ... global lat lon year month day hour;

如果是风格二，可以直接弃用，改成风格一的参数传递写法再跑。实际工程里我踩过不少global变量互相覆盖的坑，尤其是多个文件同时运行时，lat和lon这种通用名很容易被别的脚本误写。

3.3 日期时间输入处理

太阳位置程序的输入时间格式决定了它的应用边界。压缩包自带说明文件说明.txt，如果里面明确写了年份、月份、日期、小时分开传参，那这套程序就是服务于离线计算场景的。如果你需要在线实时计算，建议加一个基于MATLAB内置函数的封装层：

function [azimuth, elevation] = SolarAngleNow(lat, lon) % 获取当前UTC时间并调用原有计算函数 t = datetime('now', 'TimeZone', 'UTC'); [azimuth, elevation] = SolarAngle(lat, lon, ... year(t), month(t), day(t), hour(t), minute(t), second(t)); end

datetime对象支持TimeZone属性的设置，直接得到UTC时间可以避免本地时区带来的误差。在实际使用时，经纬度也要注意正负号约定：东经为正、西经为负，北纬为正、南纬为负。大多数程序默认这个约定，但在不同代码模块拼接时一定要统一。通常我会在函数入口加一个断言来检查经纬度范围。

3.4 skew_symmetric的实际作用

skew_symmetric.m在偏振角算法里扮演了矢量叉乘的角色。对于三维矢量，叉乘操作可以写成矩阵形式，即反对称矩阵乘以另一个矢量。代码结构如下：

function S = skew_symmetric(v) % 输入三维列向量 v，输出其反对称矩阵 S = [0, -v(3), v(2); v(3), 0, -v(1); -v(2), v(1), 0]; end

为什么要用反对称矩阵而不是直接cross(v1, v2)？原因是在偏振导航解算里，你可能需要对同一个矢量反复做叉乘，而如果把反对称矩阵存下来，每次叉乘就变成一次矩阵乘法，计算效率高且便于线性化。在求偏振方位角时，入射光E矢量和散射面法向的叉乘关系用这个函数表达非常自然。

4. 偏振角与太阳位置参数的联动关系

4.1 偏振角如何从太阳位置推导

偏振导航的基本逻辑是：太阳光在大气中散射后，散射光的偏振方向与散射平面垂直。如果你用偏振传感器测量某方向的偏振角，就可以反推出太阳相对于该方向的方位关系。这一步需要太阳位置参数作为参考。在FaiToSouth.m和FaifromSun.m这两个函数里，做的正是从“太阳矢量”到“载体航向角”的转换。

定义太阳在载体坐标系下的单位矢量为s_b，偏振传感器测得的E矢量方向为p_b，那么根据瑞利散射模型，p_b应当垂直于s_b与观测方向o_b构成的散射平面。算法实现如下：

% s_b: 太阳方向单位矢量（载体坐标系） % o_b: 观测方向单位矢量（载体坐标系） % p_b: 测量得到的偏振E矢量方向 s_cross_o = skew_symmetric(s_b) * o_b; % 散射平面法向 p_pred = s_cross_o / norm(s_cross_o); % 归一化得到预测偏振方向 ang = acos(dot(p_b, p_pred)); % 预测与实测偏振角的差

这个差值如果接近0，说明太阳位置计算准确；如果偏差过大，就要检查是不是SolarAngle.m的坐标转换出了问题。偏振角对太阳位置误差的敏感程度比强度信号高得多，太阳方位角差1度，偏振角观测残差可能被放大到3到5度。

4.2 载体有姿态时不能直接套公式

上面流程唯一的限制是——它假设载体坐标系和当地水平坐标系重合。实际应用中载体有横滚和俯仰，所以一定先要把太阳矢量从水平系转换到载体系。这里需要载体的姿态角：

% 水平系下太阳矢量（北东地或北天东，取决于代码约定） s_ned = [cos(alt) * cos(az); cos(alt) * sin(az); sin(alt)]; % 载体坐标系下太阳矢量，使用旋转矩阵 C_b_n s_b = C_b_n * s_ned;

如果你用的是北天东坐标系，公式里的分量排列要做对应调整。MATLAB的航空航天工具箱自带dcmbody2ned函数，但很多老程序是手工旋转矩阵拼出来的，检查时重点看旋转顺序是Z-Y-X还是Z-X-Y。压缩包里的FaiToSouth.m大概率封装的正是这个过程，但它的姿态输入是欧拉角还是四元数，需要打开源码确认一下。

4.3 偏振角解算的边界条件

偏振导航和太阳位置计算有一个几乎必然出现的问题：观测方向与太阳方向重合时，散射平面不唯一，偏振角退化。说白了就是当传感器直视太阳时，偏振信息几乎为零，再用偏振角反算方位就完全失效。好的程序会在这种条件下输出NaN或者置一个标志位而非继续算出一个假值。调试时如果发现偏振角跳变剧烈，先检查是否进入了退化构型。在程序包现有代码基础上加一个质量因子比较简单：

quality = abs(dot(s_b, o_b)); % 接近1表示观测方向与太阳方向接近 if quality > 0.99 polarization_valid = false; % 偏振信息越接近0越不可用 else polarization_valid = true; end

这个质量因子建议作为额外输出带回上层逻辑，用于滤波器的量测噪声自适应调整。skew_symmetric矩阵在这种情况下也会出现病态，算出来的方位角没有意义，加了质量判断后整体算法就多了一层安全保障。

5. 从验算到现场部署的四个实用技巧

5.1 用已知城市数据做基准校验

拿到程序包后的第一步不是直接跑，而是用一组已知数据校验。这里以北京为例：

春分日太阳赤纬接近0度，正午时北京太阳高度角应接近90 - 39.9042 = 50.0958度，方位角近似正南180度。如果程序输出偏差超过0.1度，就要去查均时差修正和经度到时间的换算环节。这个做基准测试的方法比对着表查天文年历高效得多，挑特殊节气日验证是最常用的。

5.2 时区传递中的字符串与数值陷阱

把UTC小时传给函数时，有一个经常出问题的细节——是否做了取余运算。比如UTC时间如果是23点，东经120度的当地太阳时约为23 + 8 = 31时，时角会超过180度，sin和cos计算没问题，但某些程序里时角没有取模就会导致结果不连续。稳妥的做法是在时角计算后统一进行角度归一化：

HA = mod(HA + pi, 2 * pi) - pi; % 时角范围归一化到 [-pi, pi]

这个mod操作对后续的atan2计算非常关键。如果不做归一化，连续跨越两天的航拍数据在午夜前后会出现方位角跳变，直接干扰偏振导航滤波器的收敛。

5.3 用符号运算验证公式正确性

MATLAB的Symbolic Math Toolbox可以帮我们快速验证公式推导有没有错误。

syms phi delta HA sin_alt = sin(phi) * sin(delta) + cos(phi) * cos(delta) * cos(HA); % 检查极端值：赤道正午（phi=0, delta=0, HA=0）时高度角应为90度 subs(sin_alt, {phi, delta, HA}, {0, 0, 0}) % 结果应为 sin(pi/2) 即 1

这种验证方式的优势在于一次检查所有极端情况，比打印数值调试快很多。另一个值得验证的点是：当phi = delta时太阳应该在天顶，方位角在这时是定义不明的，程序输出什么值都算正常，注意在算法层规避。

5.4 部署前把输出结果落盘成标准格式

在仿真链路里，太阳位置只是中间量，最终要和传感器数据一起进解算滤波器。建议把SolarAngle.m的输出格式固定为结构体或MATLAB timetable，便于和其他时间序列数据对齐。

function solar = SolarAngleStruct(lat, lon, utc_time) % 输出包含方位角、高度角和时间戳的结构体 [az, el] = SolarAngle(lat, lon, ... year(utc_time), month(utc_time), day(utc_time), ... hour(utc_time), minute(utc_time), second(utc_time)); solar = struct('time', utc_time, 'azimuth', az, 'elevation', el); end

用timetable存储后，后续可以直接用synchronize和载体姿态数据、偏振传感器数据对齐时间戳，省去大量手动循环索引的代码。对于还要跑批量仿真的人来说，这个封装结构能让整个链路干净很多。实测中我会额外记录一条delta和HA的中间变量，排错时能直接定位是在哪一步出的偏差，而不用从最终方位角反推。

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