簡介這份資源提供太陽高度角與方位角的高精度MATLAB計(jì)算代碼面向氣象學(xué)、建筑采光、地理信息系統(tǒng)及太陽能工程等領(lǐng)域的科研人員與工程師用于精確掌握太陽位置、優(yōu)化光伏板朝向或分析特定地點(diǎn)的日照條件。壓縮包共3個(gè)文件約1.28MB包含核心計(jì)算腳本、理論說明文檔與授權(quán)文件腳本負(fù)責(zé)結(jié)合日期時(shí)間與經(jīng)緯度、海拔等地球坐標(biāo)完成太陽位置解算文檔則解釋算法原理與公式背景便于理解實(shí)現(xiàn)邏輯并擴(kuò)展應(yīng)用。目前已有3714人學(xué)習(xí)下載代碼基于天文學(xué)常用公式實(shí)現(xiàn)讀者可據(jù)此快速搭建太陽位置計(jì)算模塊用于日光照射分析、建筑遮陽設(shè)計(jì)或太陽能系統(tǒng)效率評(píng)估同時(shí)加深對(duì)MATLAB科學(xué)計(jì)算與天文算法的理解。1. 太陽高度角方位角高精度計(jì)算一份能直接跑通的 MATLAB 工程代碼做光伏追光、建筑日照分析或者天文觀測(cè)的朋友大概率都遇到過同一個(gè)問題網(wǎng)上找的太陽位置公式算出來的高度角方位角總差那么一兩度短時(shí)間看不出來一旦用來控制追光支架或者做全年輻照度仿真誤差就被放大到?jīng)]法接受。這份 MATLAB 工程代碼就是沖著這個(gè)痛點(diǎn)來的它把太陽高度角和方位角的高精度計(jì)算封裝成了可直接調(diào)用的函數(shù)輸入時(shí)間、經(jīng)緯度就能拿到結(jié)果不需要額外的工具箱依賴。代碼結(jié)構(gòu)清晰注釋到位適合做新能源、氣象、測(cè)繪方向的學(xué)生和工程師直接拿去用也適合想搞明白天文算法底層邏輯的人逐行拆解。下面我從工程結(jié)構(gòu)、核心算法、運(yùn)行步驟到二次開發(fā)完整走一遍。2. 工程結(jié)構(gòu)與算法選型為什么不用工具箱自帶的函數(shù)拿到一份計(jì)算代碼第一件事不是急著跑而是先看它怎么組織的。這份工程沒有堆在一個(gè)腳本里而是按職責(zé)拆成了幾個(gè)文件這種結(jié)構(gòu)對(duì)后續(xù)改參數(shù)、換算法、做批量仿真都友好得多。2.1 目錄結(jié)構(gòu)與文件職責(zé)工程根目錄下大致是這么幾類文件SolarPosition/ ├── main_demo.m # 演示腳本跑一遍看結(jié)果 ├── solar_position.m # 核心計(jì)算函數(shù)輸入時(shí)間/經(jīng)緯度輸出高度角方位角 ├── julian_day.m # 儒略日轉(zhuǎn)換時(shí)間計(jì)算的基礎(chǔ) ├── solar_declination.m # 太陽赤緯角計(jì)算 ├── equation_of_time.m # 時(shí)差修正真太陽時(shí)與平太陽時(shí)的差 ├── local_solar_time.m # 本地真太陽時(shí)換算 ├── validate_against_data.m # 用已知天文數(shù)據(jù)做校驗(yàn) └── data/ └── sample_input.csv # 示例輸入時(shí)間序列 經(jīng)緯度main_demo.m是入口跑它就能看到一組輸出。solar_position.m是真正干活的地方其他幾個(gè)函數(shù)是它調(diào)用的子模塊。這種拆法的好處是你想換赤緯角算法只改solar_declination.m就行不用動(dòng)主邏輯想批量處理一年的數(shù)據(jù)直接循環(huán)調(diào)solar_position即可。2.2 算法選型SPA 還是簡化公式太陽位置計(jì)算有好幾個(gè)精度檔次。最簡單的用 Cooper 公式算赤緯角誤差能到 1 度以上好一點(diǎn)的用 Spencer 公式或者 PSA 算法最高精度的當(dāng)屬 NREL 的 SPASolar Position Algorithm誤差在 ±0.0003 度以內(nèi)。這份代碼走的是中間偏上的路線——用天文歷書里的經(jīng)典公式組合赤緯角和時(shí)差都做了多階諧波修正實(shí)際精度能壓到 0.01 度量級(jí)對(duì)絕大多數(shù)工程場(chǎng)景夠用了。為什么不上 SPA因?yàn)?SPA 實(shí)現(xiàn)復(fù)雜參數(shù)多代碼量大而且對(duì)普通追光或日照分析來說0.01 度和 0.0003 度的差別在實(shí)際工程里根本體現(xiàn)不出來。反過來如果只用最簡單的 Cooper 公式夏天正午可能差出 1 度多追光支架就會(huì)明顯偏。所以這個(gè)精度檔位是權(quán)衡后的選擇代碼量可控精度又足夠。提示如果你做的是天文臺(tái)級(jí)別的指向控制建議在solar_declination.m里換成 SPA 的赤緯角計(jì)算主框架不用動(dòng)。2.3 輸入輸出接口設(shè)計(jì)核心函數(shù)solar_position.m的簽名大致是這樣function [altitude, azimuth] solar_position(datetime_vec, latitude, longitude, timezone) % SOLAR_POSITION 計(jì)算太陽高度角和方位角 % 輸入 % datetime_vec - N×6 矩陣每行 [年 月 日 時(shí) 分 秒]本地時(shí)間 % latitude - 緯度度北正南負(fù) % longitude - 經(jīng)度度東正西負(fù) % timezone - 時(shí)區(qū)小時(shí)如東八區(qū)為 8 % 輸出 % altitude - N×1 高度角度地平線以上為正 % azimuth - N×1 方位角度從北順時(shí)針這個(gè)接口設(shè)計(jì)有幾個(gè)講究。時(shí)間用矩陣而不是 MATLAB 的datetime類型是為了兼容老版本也方便從 CSV 直接讀進(jìn)來。經(jīng)緯度用十進(jìn)制度不搞度分秒省得調(diào)用方還要轉(zhuǎn)換。方位角從北順時(shí)針起算這是光伏和氣象領(lǐng)域的通用約定跟羅盤讀數(shù)一致。3. 核心算法拆解從儒略日到高度角方位角這一章把計(jì)算鏈路拆開講。太陽位置計(jì)算本質(zhì)上是一條流水線時(shí)間 → 儒略日 → 太陽幾何參數(shù) → 本地坐標(biāo)轉(zhuǎn)換 → 高度角方位角。每一環(huán)都有坑下面逐個(gè)說。3.1 儒略日轉(zhuǎn)換一切時(shí)間計(jì)算的起點(diǎn)儒略日J(rèn)ulian Day是天文學(xué)里表示時(shí)間的標(biāo)準(zhǔn)方式從公元前 4713 年某個(gè)中午起算的天數(shù)。所有太陽位置公式都建立在儒略日或者儒略世紀(jì)的基礎(chǔ)上所以第一步必須把本地時(shí)間轉(zhuǎn)成儒略日。function jd julian_day(year, month, day, hour, minute, second) % 計(jì)算儒略日支持格里高利歷 if month 2 year year - 1; month month 12; end A floor(year / 100); B 2 - A floor(A / 4); % 格里高利歷修正 jd floor(365.25 * (year 4716)) floor(30.6001 * (month 1)) ... day B - 1524.5; jd jd (hour minute / 60 second / 3600) / 24; end這段代碼的關(guān)鍵在B那一項(xiàng)。儒略歷和格里高利歷的轉(zhuǎn)換差別就在這里1582 年 10 月 15 日之前用儒略歷之后用格里高利歷。工程上我們處理的都是現(xiàn)代日期所以直接按格里高利歷算就行。1524.5是基準(zhǔn)偏移量0.5對(duì)應(yīng)儒略日從中午起算的約定。如果你把0.5漏了算出來的時(shí)間會(huì)差 12 小時(shí)太陽位置直接反了。3.2 太陽赤緯角和時(shí)差修正有了儒略日接下來算兩個(gè)核心角度赤緯角太陽直射點(diǎn)緯度和時(shí)差真太陽時(shí)與平太陽時(shí)的差。function [decl, eot] solar_declination(jd) % 計(jì)算太陽赤緯角和時(shí)差 % 輸入 jd儒略日 % 輸出 decl赤緯角度eot時(shí)差分鐘 n jd - 2451545.0; % 距 J2000.0 的天數(shù) L mod(280.460 0.9856474 * n, 360); % 平黃經(jīng) g mod(357.528 0.9856003 * n, 360); % 平近點(diǎn)角 lambda L 1.915 * sind(g) 0.020 * sind(2 * g); % 黃經(jīng) epsilon 23.439 - 0.0000004 * n; % 黃赤交角 decl asind(sind(epsilon) * sind(lambda)); % 時(shí)差計(jì)算 alpha atan2d(cosd(epsilon) * sind(lambda), cosd(lambda)); eot 4 * (L - alpha); % 轉(zhuǎn)換為分鐘 eot mod(eot 720, 1440) - 720; % 歸一化到 [-720, 720] end這里用了sind、cosd、asind這些以度為單位的三角函數(shù)比用弧度再轉(zhuǎn)換清爽得多。赤緯角公式里1.915 * sind(g)和0.020 * sind(2*g)是黃經(jīng)修正項(xiàng)少了這兩項(xiàng)精度會(huì)掉到 0.1 度以上。時(shí)差那一步的mod歸一化很容易被忽略如果不做時(shí)差可能落在 [-720, 720] 之外導(dǎo)致真太陽時(shí)算錯(cuò)。3.3 本地真太陽時(shí)與地平坐標(biāo)轉(zhuǎn)換赤緯角和時(shí)差拿到后要結(jié)合經(jīng)度和時(shí)區(qū)算出本地真太陽時(shí)再轉(zhuǎn)成高度角和方位角。function [alt, az] solar_position(datetime_vec, lat, lon, tz) N size(datetime_vec, 1); alt zeros(N, 1); az zeros(N, 1); for i 1:N dt datetime_vec(i, :); jd julian_day(dt(1), dt(2), dt(3), dt(4), dt(5), dt(6)); [decl, eot] solar_declination(jd); % 本地真太陽時(shí)小時(shí) solar_time dt(4) dt(5)/60 dt(6)/3600 ... (lon - 15 * tz) / 15 eot / 60; % 時(shí)角度正午為 0下午為正 H 15 * (solar_time - 12); % 高度角 alt(i) asind(sind(lat) * sind(decl) ... cosd(lat) * cosd(decl) * cosd(H)); % 方位角 az(i) atan2d(sind(H), ... cosd(H) * sind(lat) - tand(decl) * cosd(lat)) 180; end end逐段說。solar_time那一行是核心本地時(shí)間加上經(jīng)度與時(shí)區(qū)的修正再加上時(shí)差。(lon - 15*tz)/15這一項(xiàng)把經(jīng)度差轉(zhuǎn)成小時(shí)差東八區(qū) tz8中央經(jīng)線 120 度如果你在 116 度這一項(xiàng)就是負(fù)的太陽時(shí)比北京時(shí)間晚。時(shí)角H以正午為 0每小時(shí) 15 度下午為正。高度角公式是球面三角的標(biāo)準(zhǔn)形式方位角用atan2d保證象限正確最后加 180 是因?yàn)閍tan2d算出來的是從南起算加 180 轉(zhuǎn)成從北起算。注意方位角加 180 后可能超過 360實(shí)際使用時(shí)建議再做一次mod(az, 360)。3.4 精度校驗(yàn)?zāi)靡阎獢?shù)據(jù)對(duì)一遍代碼里帶了validate_against_data.m用幾個(gè)已知的天文數(shù)據(jù)點(diǎn)做校驗(yàn)。比如春分日正午太陽赤緯角應(yīng)該接近 0高度角應(yīng)該等于 90 減去緯度。跑一遍校驗(yàn)?zāi)_本如果誤差超過 0.05 度說明某個(gè)環(huán)節(jié)有問題優(yōu)先查儒略日的0.5和時(shí)差的歸一化。4. 運(yùn)行與調(diào)試從零跑通這份代碼代碼拿到手怎么跑起來、怎么確認(rèn)結(jié)果對(duì)這一章說清楚。4.1 環(huán)境準(zhǔn)備與依賴檢查這份代碼對(duì) MATLAB 版本要求不高R2016b 以上都能跑因?yàn)橛玫搅穗[式擴(kuò)展再老的版本可能需要改循環(huán)。不需要任何工具箱sind、cosd、atan2d都是基礎(chǔ)函數(shù)。如果你用的是 Octave大部分也能兼容但datetime相關(guān)的部分要留意。跑之前先確認(rèn)路徑把SolarPosition文件夾加到 MATLAB 的搜索路徑里或者直接cd進(jìn)去。然后運(yùn)行main_demo.m。4.2 跑通演示腳本main_demo.m大致長這樣% 演示計(jì)算某地一天內(nèi)每小時(shí)太陽位置 lat 39.9; % 緯度 lon 116.4; % 經(jīng)度 tz 8; % 東八區(qū) hours (6:18); datetime_vec [repmat([2024, 6, 21], length(hours), 1), ... hours, zeros(length(hours), 2)]; [alt, az] solar_position(datetime_vec, lat, lon, tz); % 輸出結(jié)果 fprintf(時(shí)間 高度角 方位角\n); for i 1:length(hours) fprintf(%2d:00 %6.2f %6.2f\n, hours(i), alt(i), az(i)); end跑完你會(huì)看到從早上 6 點(diǎn)到下午 6 點(diǎn)的高度角和方位角變化。夏至日正午高度角應(yīng)該在 73 度左右90 減去緯度和赤緯角之差方位角從東北方向轉(zhuǎn)到西北方向。如果正午高度角明顯不對(duì)先查緯度和赤緯角的正負(fù)號(hào)。4.3 批量處理與結(jié)果導(dǎo)出實(shí)際工程里往往要算一整年、多個(gè)地點(diǎn)。把solar_position放進(jìn)循環(huán)就行% 批量計(jì)算一年逐時(shí)太陽位置 dates datetime(2024,1,1,0,0,0) : hours(1) : datetime(2024,12,31,23,0,0); datetime_vec [year(dates), month(dates), day(dates), ... hour(dates), minute(dates), second(dates)]; [alt, az] solar_position(datetime_vec, lat, lon, tz); % 導(dǎo)出 CSV T table(dates, alt, az, VariableNames, {Time, Altitude, Azimuth}); writetable(T, solar_position_2024.csv);這里把datetime類型轉(zhuǎn)成矩陣再傳入是為了保持接口統(tǒng)一。導(dǎo)出 CSV 后可以直接拿去給追光控制器或者輻照度模型用。4.4 常見報(bào)錯(cuò)與定位跑的過程中可能遇到幾類問題。一是Undefined function sind說明 MATLAB 版本太老或者路徑?jīng)]加對(duì)。二是結(jié)果全為 NaN通常是輸入矩陣?yán)镉?NaN 或者日期格式不對(duì)。三是高度角在夜間算出正值檢查時(shí)角H的范圍夜間時(shí)角可能超出 ±180 度公式本身會(huì)給出負(fù)高度角如果沒做過濾就會(huì)看到異常值。5. 避坑與常見問題那些年踩過的太陽位置計(jì)算坑這一章列幾個(gè)實(shí)際開發(fā)中反復(fù)出現(xiàn)的問題每條按現(xiàn)象、原因、解決來說。5.1 正午高度角差 1 度以上現(xiàn)象算出來的正午高度角跟天文年歷對(duì)不上差 1 度多。原因赤緯角用了簡化公式比如只取23.44 * sind(...)一項(xiàng)沒做諧波修正。解決換成solar_declination.m里的多階修正版本或者直接上 SPA 的赤緯角。5.2 方位角在正午跳變現(xiàn)象上午方位角從 90 度慢慢減小正午突然跳到 270 度附近。原因atan2d的象限判斷在高度角接近 90 度時(shí)不穩(wěn)定分母趨近于零。解決在高度角大于 89.9 度時(shí)直接令方位角為 0 或 180或者用向量法替代atan2d。5.3 時(shí)差修正漏了導(dǎo)致早晚偏差大現(xiàn)象早晨和傍晚的高度角偏差比正午大。原因時(shí)差eot沒加進(jìn)真太陽時(shí)計(jì)算或者歸一化寫錯(cuò)。解決確認(rèn)solar_time那一行里eot/60是加不是減且eot已經(jīng)歸一化到 [-720, 720] 分鐘。5.4 跨年計(jì)算時(shí)儒略日溢出現(xiàn)象算 12 月 31 日之后的數(shù)據(jù)結(jié)果突然跳到上一年的值。原因儒略日公式里year和month的調(diào)整邏輯在跨年時(shí)沒處理好。解決在調(diào)用julian_day之前先把日期標(biāo)準(zhǔn)化確保月份在 1 到 12 之間。5.5 批量計(jì)算速度慢現(xiàn)象算一年逐時(shí)數(shù)據(jù)要跑好幾分鐘。原因solar_position里用了逐行循環(huán)每次調(diào)用都重新算儒略日。解決把儒略日、赤緯角、時(shí)差的計(jì)算向量化一次性算完再進(jìn)循環(huán)速度能快十倍以上。6. 二次開發(fā)與改造把這份代碼變成你自己的工具代碼跑通只是第一步真正有價(jià)值的是按自己的場(chǎng)景改。這一章說幾個(gè)改造方向和一個(gè)我常用的驗(yàn)證習(xí)慣。6.1 替換赤緯角算法提升精度如果你做的是聚光太陽能或者天文指向0.01 度不夠可以把solar_declination.m換成 SPA 的實(shí)現(xiàn)。SPA 的赤緯角計(jì)算涉及更多項(xiàng)但接口不變主邏輯不用動(dòng)。改完后用validate_against_data.m重新校驗(yàn)確認(rèn)誤差降到 0.001 度以內(nèi)。6.2 增加大氣折射修正太陽高度角在地平線附近受大氣折射影響實(shí)際看到的位置比幾何位置高約 0.5 度。如果你的應(yīng)用涉及日出日落時(shí)刻或者低空追光需要加折射修正。常見做法是用 Bennett 公式function alt_corrected refract_correction(alt) % 大氣折射修正輸入幾何高度角度 alt_corrected alt; mask alt -1 alt 90; alt_corrected(mask) alt(mask) ... 1.02 ./ (60 * tand(alt(mask) 10.3 ./ (alt(mask) 5.11))); end這個(gè)修正只在高度角低于 5 度時(shí)明顯正午可以忽略。6.3 封裝成 Simulink 模塊或 App如果你要做實(shí)時(shí)仿真可以把solar_position封裝成 MATLAB Function 模塊放進(jìn) Simulink輸入用時(shí)鐘信號(hào)輸出直接驅(qū)動(dòng)追光模型?；蛘哂?App Designer 做個(gè)界面輸入經(jīng)緯度和日期實(shí)時(shí)顯示太陽軌跡。這兩種改造都不需要?jiǎng)雍诵乃惴ㄖ皇菗Q個(gè)調(diào)用方式。6.4 我的驗(yàn)證習(xí)慣改完任何參數(shù)我都會(huì)強(qiáng)制走一遍validate_against_data.m再手動(dòng)對(duì)幾個(gè)特殊日期春分、夏至、秋分、冬至的正午高度角。春分正午高度角應(yīng)該等于 90 減緯度夏至等于 90 減緯度加 23.44冬至等于 90 減緯度減 23.44。這幾個(gè)數(shù)對(duì)上了基本就不會(huì)有大問題。從那以后我每次改完太陽位置相關(guān)的代碼都強(qiáng)制走一遍這個(gè)校驗(yàn)流程省了很多返工。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取