静止卫星遥感图像太阳及卫星天顶、方位角(SOZ/SOA/SAZ/SAA)计算方法

各位b友,请直接跳转至最新的总结,见:https://blog.csdn.net/qq_37970770/article/details/109278770


原文

题目:静止卫星遥感图像太阳及卫星天顶、方位角(SOZ/SOA/SAZ/SAA)计算方法
作者:白艺亭
日期:2020.02.25
摘要

遥感图像角度数据(太阳天顶角、太阳方位角、卫星天顶角、卫星方位角)是遥感图像基础数据的重要组成部分,在遥感图像预处理和遥感应用如气溶胶光学厚度反演中是必不可少的参数之一。国产遥感数据经常缺失角度数据支持,大多数产品都不给出成像时对应的逐像元角度数据,需要用户手动计算,而关于角度数据的定义及相应计算方法,网络上很多纷繁复杂的公式,且很多都没有参考来源,本文从遥感专业的角度,在查阅且推导了很多国际英文文献中角度计算公式后,给出四个角度的正确计算公式,为遥感图像的角度数据计算提供支持。

关键词 : SOZ、SOA、SAZ、SAA计算; 遥感图像

1引言

一般遥感数据分为静止卫星数据和极轨卫星数据。国际上的静止卫星数据大多都带有角度数据(如H8 H9 GK2A ),位于波段数据之后,而例如我国的风云四号静止卫星则没有在0.5-2km空间分辨率数据后设置角度数据,因此如果要对其进行预处理(如大气校正)和气溶胶光学厚度反演时则需要逐像元的角度数据。对于小幅宽(小于250km)的极轨卫星数据,一般在进行预处理和气溶胶光学厚度反演应用时,我们一般选择成像时中心像元的太阳天顶角、太阳方位角、卫星天顶角、卫星方位角代表所有像元的该角度值,而在大幅宽的极轨卫星数据时,由于太阳及卫星观测角度大有不同,且受像元上空的大气空间异质性特征影响,无法利用中心像元代表所有像元,此时则需要计算极轨卫星逐像元角度数据。

本文计算太阳的天顶和方位角在极轨和静止卫星都是可以用的,但卫星的天顶和方位角只适用于静止卫星,由于涉及到星下点经纬度,在极轨卫星上是否能用还有待我研究一下,目前还不敢确切说,以免有任何不必要的错误,因为本文卫星天顶和方位角所有公式的参考文献都是出自静止卫星的应用。

2计算方法

网络上有关四个角度数据计算的方法成出不群,但大都没有参考文献来源,且对于方位角的定义千奇百怪,导致很多公式计算下来的结果并不可靠。本文从遥感专业角度,就以上问题于数月前开展相关研究,在查阅大量英文文献后,经推导筛选验证,最终选择确定了可靠的计算方法,现给出其定义及计算公式。如果小伙伴对推导感兴趣,我们可以就文章中的内容进行探讨。
目前对于一些概念,遥感界有太多纷繁复杂、各式各样奇奇怪怪没有统一的定义了,或许我们需要将一些概念统一一下,使其成为经典概念而为后世所学。

2.1太阳天顶角

太阳天顶角(Solar zenith angle, SOZ):像元位置处该时入射太阳光线与该像元地平面法线之间的夹角。取值范围:0°- 90°,有的地方写为0-180°,注意我们不采用,依然使用0-90°.

正确的计算公式为:

计算公式在[1] [2]文章中:
【[1]中位于第一章15页,[2]中有详细推导过程,其计算本质上是利用纳皮尔准则Napier Rules在球形三角学中的应用推导出来的),关于纳皮尔准则球形三角学的知识参考[3] [4]仅供阅读了解 .】
在这里插入图片描述
式中,在这里插入图片描述为太阳天顶角,为像元纬度(Latitude),为像元太阳赤纬角,在这里插入图片描述为太阳时角。
太阳赤纬角为该时间太阳直射点与地心的连线 和 赤道平面所成的夹角,每年在-23.45°和23.45°之间变化,3.21 6.22 9.23 12.22分别为北半球一年中的春分,夏至,秋分和冬至日,直射点分别位于赤道,北回归线,赤道,南回归线。太阳赤纬角计算公式[2](参考文献的第二章56页)为:
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
AST为像元当地时,一般如果知道成像时间北京时或者其他地方时的话可以根据经度差进行逐像元AST计算,例如为以北京时中心,成像为BT10:30,则计算AST在这里插入图片描述

2.2太阳方位角

太阳方位角(Solar azimuth angle, SOA):由正北方向起,顺时针转至某个时刻太阳在该像元的入射光线在地表的投影所在的射线,顶点为地表像元点,这个旋转的角度称为太阳方位角,或者,由正北方向起,顺时针转至太阳与像元的连线(将此时的太阳看做一个小球)在地表的投影所在的射线,顶点为地表像元点,取值范围为0° - 360°. 值得注意的是南半球太阳方位角定义一样,不会因为南半球而重新定义。
有的地方会定义为像元正南正北连线以东的太阳方位角角度为负(与上午时角为负有关),以西的太阳方位角角度为正(与下午时角为正有关),正负在公式里利用符号函数进行判断,
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
取值范围为-180° - 180°,此种定义常见于很多非遥感专业文献中,例如此种表达也出现在参考文献【1】中,不建议采用这种计算方法。

正确的计算公式为:

在这里插入图片描述
为太阳方位角,在这里插入图片描述为太阳高度角,取值范围为0-90,太阳入射光线与像元地平面的夹角。其他变量的含义同2.1.

2.3卫星天顶角

卫星天顶角(satellite zenith angle,SAZ):卫星传感器与像元的连线 和 像元地平面法线所成的夹角,取值范围为0°到90°.有的地方卫星天顶角定义如下,取值范围为-90°到90°:
在这里插入图片描述
我们仍不采用,仍然采用前述0-90°的定义方式。

正确的计算公式为参考文献6中的系列公式(1)-(4):

见以下截图:

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述即为卫星天顶角SAZ.
式中,在这里插入图片描述为纬度,等于Lat.在这里插入图片描述为卫星星下点所在经度,在这里插入图片描述为像元经度,等于Lon(Longitude).
R为地球半径,r为卫星到地心的距离,d为卫星到像元的距离.
在这里插入图片描述为卫星高度角.

2.4卫星方位角

卫星方位角(satellite azimuth angle,SAA):由正北方向起,顺时针转至某个时刻卫星和像元连线在地表的投影所在的射线,顶点为地表像元点,这个旋转的角度称为像元的卫星方位角,取值范围为0° - 360°.
有的地方会定义为当地正南正北连线(子午线)以东的角度为负,以西的角度为正,取值范围为-180° - 180°,此种定义也不建议采用。

正确的计算公式为:

根据纳皮尔准则,计算公式为:
(仍见参考文献[6])
在这里插入图片描述
Lat为像元纬度。

在这里插入图片描述为球面角在这里插入图片描述
原文文章内容截图:

在这里插入图片描述

和SAA的关系如原文中所述,要看像元和卫星星下点的相对位置,分以下四种情况,下面详细介绍求解过程:
在这里插入图片描述,式中,
在这里插入图片描述为像元经度与星下点经度差。为像元纬度与星下点纬度差。

(1)卫星星下点位于像元P’的东北NE方向,包含正北这一方向时,即像元经度小于等于星下点经度,像元纬度小于星下点纬度。
在这里插入图片描述<=0,且在这里插入图片描述<0时,在这里插入图片描述
(2)卫星星下点位于像元P’的东南SE方向,包括正东方向,即像元经度小于星下点经度,像元纬度大于等于星下点纬度。

在这里插入图片描述<0,且在这里插入图片描述>=0时,在这里插入图片描述
(3)卫星星下点位于像元P’的西南SW方向,包括正南方向,即像元经度大于等于星下点经度,纬度大于星下点纬度。
在这里插入图片描述>=0,且在这里插入图片描述>0时,在这里插入图片描述
(4)卫星星下点位于像元P’的西北NW方向,包括正西方向,即像元经度大于星下点经度,纬度小于等于星下点纬度。
在这里插入图片描述>0,且在这里插入图片描述<=0时,在这里插入图片描述

NOTE:

1.以上所有计算均假设地球是个球体,若是把地球当作椭球体,计算公式比较繁琐,但其计算精度会更高;
2.像元所代表地平面视为水平面(遥感当中一般都这么假设),若是视为垂直面或者有一定倾角的面,则上述天顶角和方位角的计算公式不普遍成立;
3.每一个角度求反函数前,若反函数为acos、asin等时,注意被反数组取值范围为[-1,1],否则就会有错误产生。

3结果

以国产风云四号静止卫星为例,利用上述公式计算结果与官方给出4km的角度数据结果一致。
在这里插入图片描述

4总结

本文是作者2019年夏天完成的工作,一直没有总结整理,加上对于很多计算公式的正确与否一直处于判断状态,正直疫情在家时期,将其整理成篇,希望帮助到有需要的人。值得指出的是,在最终求解出来的角度值如有Nan值,则需要检查是否在上一步反函数的被反数组中出现了三角函数最大最小值之外的值,将其修正即可。

创作不易,耗费了整整几天的功夫整理、撰写,如蒙喜欢,麻烦点赞关注一键二联,不准白嫖~

《晴雯歌》红楼梦乐曲


参考文献:
[1]Duffie J A, Beckman W A, Blair N. Solar Engineering of Thermal Processes, Photovoltaics and Wind[M]. John Wiley & Sons, 2020.
[2]William B. Stine, Michael Geyer. Power From The Sun[M].NET,2001.
[3]百度百科.球面三角学.
[4]百度文库.球面三角学.
[5]Kalogirou S A. Solar energy engineering: processes and systems[M]. Academic Press, 2013.
[6]Soler T, Eisemann D W. Determination of look angles to geostationary communication satellites[J]. Journal of surveying engineering, 1994, 120(3): 115-127.

以下是用Python实现计算太阳天顶角、方位角卫星天顶角、方位角的代码: ```python import math # 计算太阳天顶角 def solarZenithAngle(latitude, longitude, year, month, day, hour, minute, second): # 计算儒略日 J0 = 367 * year - int((7 * (year + int((month + 9) / 12.0))) / 4.0) + int(275 * month / 9.0) + day + 1721013.5 UT = hour + minute / 60.0 + second / 3600.0 T = (J0 - 2451545.0 + UT / 24.0) / 36525.0 # 计算太阳平均位置 L0 = 280.46645 + 36000.76983 * T + 0.0003032 * T**2 M = 357.52910 + 35999.05030 * T - 0.0001559 * T**2 - 0.00000048 * T**3 C = (1.914600 - 0.004817 * T - 0.000014 * T**2) * math.sin(math.radians(M)) + (0.019993 - 0.000101 * T) * math.sin(math.radians(2 * M)) + 0.000290 * math.sin(math.radians(3 * M)) sunLongitude = L0 + C # 计算太阳赤纬角 epsilon = 23.439 - 0.0000004 * T alpha = math.atan2(math.cos(math.radians(epsilon)) * math.sin(math.radians(sunLongitude)), math.cos(math.radians(sunLongitude))) delta = math.asin(math.sin(math.radians(epsilon)) * math.sin(math.radians(sunLongitude))) # 计算太阳时角 hourAngle = math.radians(15.0 * (12.0 - longitude / 15.0 - UT)) # 计算太阳高度角 solarZenithAngle = math.acos(math.sin(math.radians(latitude)) * math.sin(delta) + math.cos(math.radians(latitude)) * math.cos(delta) * math.cos(hourAngle)) return math.degrees(solarZenithAngle) # 计算太阳方位角 def solarAzimuth(latitude, longitude, year, month, day, hour, minute, second): # 计算儒略日 J0 = 367 * year - int((7 * (year + int((month + 9) / 12.0))) / 4.0) + int(275 * month / 9.0) + day + 1721013.5 UT = hour + minute / 60.0 + second / 3600.0 T = (J0 - 2451545.0 + UT / 24.0) / 36525.0 # 计算太阳平均位置 L0 = 280.46645 + 36000.76983 * T + 0.0003032 * T**2 M = 357.52910 + 35999.05030 * T - 0.0001559 * T**2 - 0.00000048 * T**3 C = (1.914600 - 0.004817 * T - 0.000014 * T**2) * math.sin(math.radians(M)) + (0.019993 - 0.000101 * T) * math.sin(math.radians(2 * M)) + 0.000290 * math.sin(math.radians(3 * M)) sunLongitude = L0 + C # 计算太阳赤纬角 epsilon = 23.439 - 0.0000004 * T alpha = math.atan2(math.cos(math.radians(epsilon)) * math.sin(math.radians(sunLongitude)), math.cos(math.radians(sunLongitude))) delta = math.asin(math.sin(math.radians(epsilon)) * math.sin(math.radians(sunLongitude))) # 计算太阳时角 hourAngle = math.radians(15.0 * (12.0 - longitude / 15.0 - UT)) # 计算太阳高度角 solarZenithAngle = math.acos(math.sin(math.radians(latitude)) * math.sin(delta) + math.cos(math.radians(latitude)) * math.cos(delta) * math.cos(hourAngle)) # 计算太阳方位角 solarAzimuth = math.atan2(math.sin(hourAngle), math.cos(hourAngle) * math.sin(math.radians(latitude)) - math.tan(delta) * math.cos(math.radians(latitude))) return math.degrees(solarAzimuth) # 计算卫星天顶角 def satelliteZenithAngle(satelliteElevation): return 90.0 - satelliteElevation # 计算卫星方位角 def satelliteAzimuth(satelliteElevation, satelliteAzimuth, latitude, longitude): az = math.atan2(math.sin(math.radians(satelliteAzimuth)), math.cos(math.radians(satelliteAzimuth)) * math.sin(math.radians(latitude)) - math.tan(math.radians(satelliteElevation)) * math.cos(math.radians(latitude))) if az < 0: az += 2.0 * math.pi return math.degrees(az) ``` 以上代码,`latitude`表示观测点的纬度,`longitude`表示观测点的经度,`year`、`month`、`day`、`hour`、`minute`、`second`表示观测时间,`satelliteElevation`表示卫星仰角,`satelliteAzimuth`表示卫星方位角。函数返回值分别为太阳天顶角、太阳方位角卫星天顶角和卫星方位角
评论 19
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值