太阳影子定位：基于Matlab的数学建模与优化求解全解析

1. 从一道经典赛题说起：太阳影子定位的工程魅力

2015年高教社杯全国大学生数学建模竞赛的A题，题目是“太阳影子定位”。这道题当年难倒了不少队伍，但也让很多同学第一次真切地感受到，数学和编程是如何联手解决一个看似“玄学”的实际问题的。简单来说，题目给了你一段视频，视频里有一根直杆，太阳照在杆子上会投下影子，影子随着时间在动。题目要求你，仅仅通过分析视频中影子长度的变化，反推出这根杆子所在的地理位置（经纬度）以及拍摄视频的日期。

听起来是不是有点像侦探破案？没错，这就是数学建模的魅力——把现实世界模糊的、连续的现象，抽象成精确的、可计算的数学模型。我当时带学生做这道题，最大的感触不是最后的答案有多准，而是整个从物理原理到代码实现，再到结果优化的过程，充满了工程实践的乐趣和挑战。今天，我就以一个过来人的视角，掰开揉碎了讲讲这道题的Matlab求解思路与核心代码实现。无论你是正在备赛的同学，还是对天文地理算法感兴趣的朋友，这篇文章都能给你一套可以直接“抄作业”的完整方案。

2. 问题本质拆解：从影子到坐标的数学桥梁

拿到这种问题，最忌讳的就是一头扎进代码里。我们先得把问题看清楚，拆明白。太阳影子定位的核心，其实是一个“正问题”和“反问题”的结合。

正问题（Forward Problem）：如果我知道拍摄地点的经纬度、日期时间、杆子高度，我一定能计算出任何时刻影子的长度和方向。这个过程是确定的，有成熟的公式（太阳高度角、方位角计算公式）。

反问题（Inverse Problem）：题目给的是结果（影子长度随时间变化的序列），要求我们反推原因（地点、日期）。这是问题的难点，因为同一个影子变化曲线，理论上可能对应多个不同的（经纬度，日期）组合，这就是所谓的“多解性”或“不适定性”。

所以，我们的解题策略就清晰了：

1. 建立正问题模型：用数学公式严谨地描述“太阳-地球-杆子-影子”这个系统。
2. 将反问题转化为优化问题：我们猜测一个地点和日期，用正问题模型算出这个猜测下的影子变化曲线，然后和题目给出的真实影子曲线进行比较。两者越接近，说明我们的猜测越准。那么，寻找最准的猜测，就变成了一个“让计算曲线和真实曲线差异最小”的优化问题。
3. 用算法求解优化问题：利用Matlab强大的优化工具箱，让计算机自动去搜索那个让误差最小的（经纬度，日期）。

下面，我们就沿着这三个步骤，一步步展开。

2.1 正问题模型：太阳位置计算是基石

一切的基础是计算太阳在天空中的位置，具体表现为两个角：太阳高度角（Altitude）和太阳方位角（Azimuth）。高度角决定了影子有多长，方位角决定了影子指向哪里。

计算这两个角需要一系列参数和公式，涉及天文、地理和时间的转换。这里我给出最核心的计算步骤和对应的Matlab函数思维，省略一些过于繁琐的中间推导。

关键参数：

• 儒略日（Julian Day, JD）：这是天文学中一个连续计时系统，是计算所有天文现象的时间基准。把我们的普通日期时间（年、月、日、时、分、秒）转换成儒略日是第一步。
• 平太阳时（Mean Solar Time）与真太阳时（Apparent Solar Time）：我们手表的时间是“平太阳时”，但地球公转轨道是椭圆且自转轴有倾斜，导致真太阳在天空中的运动并不均匀。它们的差值叫做“时差（Equation of Time）”。计算太阳位置需要用真太阳时。
• 太阳赤纬（Declination, δ）：太阳直射点纬度，随日期变化。
• 太阳时角（Hour Angle, ω）：相对于当地正午的时间角度。

核心计算公式（简化版）：

对于一个给定的地点（经度Long，纬度Lat）和世界时（UTC）t，太阳高度角α和方位角A的计算遵循以下关系：

sin(α) = sin(Lat)*sin(δ) + cos(Lat)*cos(δ)*cos(ω) cos(A) = (sin(δ) - sin(Lat)*sin(α)) / (cos(Lat)*cos(α)) // 注意方位角象限的判断

其中，时角ω与当地真太阳时有关。有了高度角α，一根高度为H的直杆，其影子长度L就非常简单了：

L = H / tan(α)

只要α > 0（白天），这个公式就成立。

在Matlab里，我们会把这一系列计算封装成一个函数，比如叫calculateShadowLength。这个函数的输入是（经纬度，日期时间，杆高），输出就是影子长度。

function L = calculateShadowLength(lat, lon, date_vec, H) % lat, lon: 纬度和经度（度） % date_vec: 日期时间向量 [年, 月, 日, 时, 分, 秒]， 注意时区！通常题目数据是北京时间（东八区），需要先转UTC。 % H: 杆高（米） % L: 影子长度（米） % 1. 将北京时间转换为UTC时间（如果题目数据是北京时间） UTC_vec = date_vec; UTC_vec(4) = UTC_vec(4) - 8; % 东八区转UTC，减8小时 % 2. 计算儒略日 JD = getJulianDay(UTC_vec); % 3. 计算太阳赤纬δ、时差EoT等参数（这里调用子函数） [delta, EoT] = calculateSolarParameters(JD); % 4. 计算本地真太阳时 LST = calculateLocalApparentSolarTime(lon, UTC_vec, EoT); % 5. 计算太阳时角ω omega = 15 * (LST - 12); % 每小时15度，正午时角为0 % 6. 计算太阳高度角α lat_rad = deg2rad(lat); delta_rad = deg2rad(delta); omega_rad = deg2rad(omega); sin_alpha = sin(lat_rad)*sin(delta_rad) + cos(lat_rad)*cos(delta_rad)*cos(omega_rad); alpha = asin(sin_alpha); % 高度角，弧度 % 7. 计算影子长度L L = H / tan(alpha); end

注意：上面的代码是一个高度简化的框架，getJulianDay,calculateSolarParameters等子函数需要你根据完整的天文算法实现。网络上可以找到诸如“Meeus天文算法”的Matlab实现，这是最权威的之一。在竞赛中，使用经过验证的代码片段是明智的。

2.2 误差函数的构建：衡量“猜得准不准”

正问题模型是我们的武器。现在，假设我们从视频里提取出了一系列数据：在时间点t1, t2, ..., tn，测量到的影子长度是L1_meas, L2_meas, ..., Ln_meas。

我们的目标是找到一组参数p = [纬度， 经度， 日期（通常用年积日表示）]，使得由这组参数通过正问题模型计算出的影子长度L_calc，与测量值L_meas的总体误差最小。

最常用的误差函数是残差平方和（Sum of Squared Residuals, SSR）：

SSR(p) = Σ [ L_calc(t_i; p) - L_meas(t_i) ]^2

我们的任务就是寻找参数p，使得SSR(p)这个值达到最小。

在Matlab中，我们将其写成一个函数，供优化器调用：

function error = shadowError(params, time_list, measured_lengths, H) % params: 待优化的参数向量 [latitude, longitude, day_of_year] % time_list: 测量时间点列表（已转换为datetime或序列化格式） % measured_lengths: 对应的测量影子长度 % H: 杆高 % error: 残差平方和 lat = params(1); lon = params(2); day_of_year = round(params(3)); % 日期参数通常取整 % 将“年积日”转换为具体的年月日（需要知道年份，年份有时也是未知参数或可假设） % 假设我们知道年份，例如2015年 year = 2015; date_vec = datevec(datenum(year, 1, day_of_year)); % 获取该年该天的0时0分0秒 calculated_lengths = zeros(size(time_list)); for i = 1:length(time_list) % 组合完整的日期时间：日期 + 当天的时间 current_time_vec = date_vec; current_time_vec(4:6) = [hour(time_list(i)), minute(time_list(i)), second(time_list(i))]; % 计算该时刻的影子长度 calculated_lengths(i) = calculateShadowLength(lat, lon, current_time_vec, H); end % 计算残差平方和 residuals = calculated_lengths - measured_lengths; error = sum(residuals.^2); end

3. 优化求解策略：如何让计算机找到最优解

构建好误差函数后，我们就把它扔给Matlab的优化器。但这里有几个非常关键的技巧，直接决定了你能否找到正确答案，以及找得快不快。

3.1 优化算法的选择

这是一个多参数、非线性、可能存在多个局部极小值的优化问题。Matlab的fmincon（约束优化）或lsqnonlin（非线性最小二乘）是常用选择。我个人更倾向于lsqnonlin，因为我们的误差本质就是最小二乘形式，这个函数专门为此设计，效率更高。

% 假设已有数据：times, measured_L, H % 设置初始猜测值 p0 = [初始纬度， 初始经度， 初始年积日] p0 = [30, 120, 180]; % 例如：猜在北纬30度，东经120度，年中左右 % 设置参数边界（非常重要！） lb = [-90, -180, 1]; % 下限：纬度-90到90，经度-180到180，年积日1到365 ub = [90, 180, 365]; % 调用 lsqnonlin % 注意：我们需要调整误差函数，使其返回残差向量而非标量和 options = optimoptions('lsqnonlin', 'Display', 'iter', 'Algorithm', 'trust-region-reflective'); [p_opt, resnorm, residual, exitflag, output] = lsqnonlin(@(p) myResidualFunc(p, times, measured_L, H), p0, lb, ub, options); function residuals = myResidualFunc(params, time_list, measured_lengths, H) % 这个函数返回残差向量，而不是平方和 % ... 内部计算 calculated_lengths ... residuals = calculated_lengths - measured_lengths; end

优化得到的p_opt就是我们认为最可能的 [纬度， 经度， 年积日]。

3.2 初始值与搜索范围的玄学

这是本题最大的坑，也是体现经验的地方。优化算法像是一个蒙着眼睛的登山者，你告诉它要找最低点（最小误差）。如果你把它放在青藏高原（一个糟糕的初始点），它可能很快找到旁边的青海湖（一个局部最优点）就停了，却永远找不到真正的太平洋海底（全局最优点）。

如何设置好的初始值？

1. 利用常识：题目视频通常在中国境内拍摄，所以纬度大概在20°N到50°N，经度在80°E到130°E之间。日期可以根据影子变化快慢、正午影子长短有个大致判断（夏季影子短，冬季影子长；中午影子短，变化慢）。
2. 网格搜索（Brute Force Grid Search）：在可能的经纬度、日期范围内，按一定间隔（如经纬度5度一格，日期10天一格）遍历所有组合，计算每个点的误差。选出误差最小的几个点作为fmincon或lsqnonlin的初始值。这个方法计算量大，但非常稳妥，能有效避免陷入局部最优。可以用parfor进行并行计算加速。
3. 利用影子方向：如果视频还能提取影子方向（而不仅仅是长度），那约束力就强太多了。方向信息对经度非常敏感。你可以先只用方向数据做一个粗略定位，再用这个结果作为长度优化的初始值。

搜索范围（边界）的设置同样关键：

• 纬度必须限定在[-90, 90]。
• 经度通常限定在[-180, 180]或[0, 360]，注意一致性。
• 日期（年积日）限定在[1, 365]（或366）。
• 强烈建议：如果你大致能判断在北半球，可以把纬度下限定为0。如果你知道是东经，可以把经度下限定为0。这能极大地缩小搜索空间，提高优化效率和准确性。

3.3 多解性与结果验证

由于影子长度变化曲线可能相似，优化结果可能存在多解。例如，在春分/秋分前后，南北半球对称纬度上的太阳高度角变化可能很像。如何增加确定性？

1. 加入先验信息：如果题目暗示或你知道视频拍摄于中国某城市，那么结果应该接近该城市的坐标。
2. 利用日期约束：如果你能从视频背景、植被、人物衣着等非数学模型信息中推断出大致季节，就可以极大地缩小日期搜索范围，从而排除掉另一个季节的“镜像解”。
3. 敏感性分析：在得到最优解后，轻微扰动参数（如纬度±1度），看误差是否急剧增大。如果误差变化平缓，说明这个方向上的解可能不唯一；如果变化陡峭，说明解比较稳定。
4. 可视化验证：将最优参数代入模型，计算出全天的影子长度变化曲线，与测量数据点画在同一张图上。肉眼观察拟合程度。一个好的拟合，点应该紧密分布在曲线两侧。

4. 完整实现流程与代码框架

结合以上所有分析，我给出一个更贴近实战的、完整的代码执行流程框架。请注意，以下代码需要你填充具体的太阳位置计算细节。

%% 主程序：太阳影子定位 clear; clc; close all; % 步骤1：数据准备（这里需要你从题目附件或自己模拟数据） % 假设我们已经从视频中提取出以下数据： % time_str_list: 时间字符串单元格数组，如 {'14:00:00', '14:10:00', ...} (北京时间) % measured_L: 对应的影子长度测量值（米） % H: 已知的杆子高度（米） % 示例模拟数据（用于测试流程） H = 3; % 杆高3米 % 假设真实位置：北纬39.9度，东经116.4度（北京附近），日期：2015年6月1日（年积日152） [sim_times, sim_L] = simulateShadowData(39.9, 116.4, 152, H); % 你需要实现这个模拟函数 time_str_list = sim_times; % 用模拟数据代替真实提取数据 measured_L = sim_L; % 将北京时间字符串转换为datetime数组，并考虑时区转换 zone_diff = 8; % 东八区 times_utc = datetime(time_str_list, 'InputFormat', 'HH:mm:ss') - hours(zone_diff); % 步骤2：定义优化问题 % 误差函数（返回残差向量） resid_func = @(params) calcResiduals(params, times_utc, measured_L, H, 2015); % 初始猜测（基于网格搜索得到的最佳点，这里手动设置为例） initial_guess = [35, 110, 150]; % [lat, lon, day_of_year] % 参数边界 lb = [0, 70, 1]; % 中国境内大致范围：北纬，东经，年积日 ub = [55, 140, 365]; % 步骤3：执行优化 options = optimoptions('lsqnonlin', 'Display', 'final', ... 'MaxFunctionEvaluations', 3000, ... 'MaxIterations', 1000, ... 'FunctionTolerance', 1e-10, ... 'StepTolerance', 1e-10); [opt_params, resnorm, residuals, exitflag, output] = ... lsqnonlin(resid_func, initial_guess, lb, ub, options); fprintf('优化结果：\n'); fprintf('纬度: %.4f°N\n', opt_params(1)); fprintf('经度: %.4f°E\n', opt_params(2)); fprintf('年积日: %.0f (大约对应日期: %s)\n', opt_params(3), ... datestr(datenum(2015,1,round(opt_params(3))), 'yyyy-mm-dd')); fprintf('残差平方和: %.6e\n', resnorm); % 步骤4：结果可视化 % 4.1 绘制拟合曲线 figure(1); [calc_times, calc_L] = simulateShadowData(opt_params(1), opt_params(2), round(opt_params(3)), H); plot(sim_times, measured_L, 'bo', 'DisplayName', '测量数据（模拟）'); hold on; plot(sim_times, calc_L, 'r-', 'LineWidth', 1.5, 'DisplayName', '模型拟合'); xlabel('时间（北京时间）'); ylabel('影子长度 (m)'); title('太阳影子长度变化：测量 vs. 模型拟合'); legend('Location', 'best'); grid on; % 4.2 绘制残差图 figure(2); plot(sim_times, residuals, 'ks-', 'MarkerFaceColor', 'k'); xlabel('时间（北京时间）'); ylabel('残差 (m)'); title('拟合残差图'); hline = refline(0,0); hline.Color = 'r'; hline.LineStyle = '--'; grid on; %% 辅助函数：计算残差向量 function residuals = calcResiduals(params, times_utc, measured_L, H, year) lat = params(1); lon = params(2); day_of_year = round(params(3)); % 将年积日转换为该天的0时 base_date = datetime(year, 1, 1) + days(day_of_year - 1); base_date_vec = datevec(base_date); n = length(times_utc); calc_L = zeros(n, 1); for i = 1:n % 组合日期和时间 current_time_utc = times_utc(i); [~, ~, ~, hh, mm, ss] = datevec(current_time_utc); full_date_vec = base_date_vec; full_date_vec(4:6) = [hh, mm, ss]; % 调用正问题模型计算影子长度 calc_L(i) = calculateShadowLength(lat, lon, full_date_vec, H); end residuals = calc_L - measured_L; end %% 辅助函数：模拟生成影子数据（用于测试，替代真实数据提取） function [time_str_list, shadow_lengths] = simulateShadowData(lat, lon, day_of_year, H) year = 2015; % 生成一天内从上午10点到下午14点，间隔10分钟的数据 start_time = datetime(year,1,day_of_year,10,0,0); end_time = datetime(year,1,day_of_year,14,0,0); interval = minutes(10); time_vec = start_time:interval:end_time; n = length(time_vec); shadow_lengths = zeros(n,1); time_str_list = cell(n,1); base_date_vec = datevec(datetime(year,1,day_of_year)); for i = 1:n t = time_vec(i); [~,~,~,h,m,s] = datevec(t); full_vec = base_date_vec; full_vec(4:6) = [h,m,s]; % 注意：模拟时我们假设已知位置，所以用同样的calculateShadowLength函数 % 但在真实解题中，这个函数是待求的逆过程。 shadow_lengths(i) = calculateShadowLength(lat, lon, full_vec, H); time_str_list{i} = datestr(t, 'HH:MM:SS'); end end

5. 实战中的坑与经验之谈

代码框架有了，但真正做的时候，你会发现一堆“坑”。下面是我总结的几个关键点：

坑1：时间系统的混淆这是最常见、最致命的错误。题目给的时间是北京时间（东八区），而天文计算通常使用世界协调时（UTC）或力学时。你必须进行时区转换。更精细的，计算真太阳时还要考虑“时差”和“经度修正”。一个简单的处理流程是：

北京时间 -> 减去8小时 -> UTC时间 -> 加上时差(EoT) -> 真太阳时

在calculateShadowLength函数开头做这个转换。很多队伍结果偏差几度，根源就在这里。

坑2：杆高H的未知处理原题中，杆高H可能是已知的，也可能是未知的。如果未知，你需要把它也作为一个待优化参数，加入params向量中，比如params = [lat, lon, day, H]。这时，模型的自由度增加，优化难度也会加大，对初始值更加敏感。

坑3：测量数据的噪声与提取误差从视频中提取影子长度，本身就有误差。像素比例尺的标定、影子端点的判断都会引入噪声。你的模型误差（SSR）不可能为0。因此，在优化时，可以适当放宽停止条件（如FunctionTolerance），避免优化器在噪声中过度拟合。另外，考虑使用鲁棒性更好的误差函数，比如Huber损失，而不是简单的平方和，可以减少个别异常数据点对整体结果的影响。

坑4：优化陷入局部最优这是反问题的固有难题。除了前面提到的网格搜索找初始值，还可以：

• 多起点优化：用不同的初始值多跑几次优化，比较最终的结果和残差。如果多个差异很大的初始点都收敛到同一组参数附近，那这组参数就很可能是全局最优。
• 使用全局优化算法：如GlobalSearch或MultiStart，它们会在多个初始点启动局部优化器（如fmincon），增加找到全局解的概率。但计算成本较高。
• 分步优化：先固定日期，只优化经纬度；或者先固定经纬度，只优化日期。通过降低维度来简化问题，再用得到的结果作为全参数优化的初始值。

坑5：结果的地理合理性算出来的经纬度，一定要放到地图上看一眼！如果定位到太平洋中心或西伯利亚荒原，那大概率是错的。结合地理常识进行判断，是最后一道，也是最重要的一道检验。

最后，我想说，这道“太阳影子定位”题之所以经典，是因为它完美地诠释了数学建模的全过程：物理建模 -> 数学抽象 -> 算法实现 -> 数值求解 -> 结果分析。它不要求你发明新算法，但极其考验你将理论知识转化为代码、并处理各种实际细节（时间、坐标、优化）的工程能力。希望这篇超详细的拆解，能帮你打通任督二脉。代码是骨架，背后的原理和思考才是灵魂。祝你建模顺利！