CN103605887B - 一种在轨卫星成像波浪水面太阳耀光时空分布仿真模型 - Google Patents

Abstract

本发明涉及一种在轨卫星成像波浪水面太阳耀光时空分布仿真模型，用于判断卫星传感器在观测水面时是否会受到太阳耀光的影响以及影响程度，其步骤如下：根据天文学知识计算全球任意时间与地点的太阳天顶角与方位角；根据卫星观测日期、时间、地面点经纬度、卫星星历与轨道参数，计算卫星传感器的观测几何，进而计算星下点轨迹；根据风速与风向等海面状况参数计算波浪水面的二向反射率，用于描述波浪水面的太阳耀光分布模式；基于大气辐射传输模型，进一步计算大气层顶部的太阳耀光的时空分布模式。本发明分别能计算卫星在对地观测时在全球范围内是否会出现太阳耀光、特定区域出现太阳耀光的时间段与太阳耀光的分布模式等。此外本发明是一种基于几何模型推演的优化算法，具有精度高、运算速度快等优点，且对卫星轨道参数没有特殊要求，具有业务化推广应用的前景。

Description

一种在轨卫星成像波浪水面太阳耀光时空分布仿真模型

技术领域

本发明涉及一种在轨卫星成像波浪水面太阳耀光时空分布模型，属于光学遥感领域，在海洋水色遥感与海面目标遥感探测方面具有重要意义。

背景技术

波浪水面太阳耀光是一种常见的自然现象，是水面镜面反射太阳辐射的结果。静止水面可以看作光滑的表面，当强光源入射时，在与入射线方向关于法线对称的反射线方向上形成强烈的镜面反射。当有风浪存在时，水面形成倾角而造成太阳耀光(或称太阳耀斑)。当水面出现太阳耀光时，反射能量剧增，严重干扰了目标自身的信号，给水环境遥感监测与水下生态系统研究、海面目标识别等带来极大的困难。

如果对遥感(或航摄)时间、太阳高度、太阳方位与航向等相互关系缺乏周全的设计，将会导致水体遥感的失败。因此能否避开镜面反射，使遥感影像能正确地反映水质的光学特性差异，是水体遥感的关键问题之一。对于海面目标监测与跟踪，当海面大型舰船淹没在海面太阳耀光背景当中时，导致舰船目标和背景海面难以区分。由于水面太阳耀光掩盖了水体的真实物理特性，严重干扰了水下与水面特征的获取。因此需要对太阳耀光的空间分布模式进行深入探讨。不同位置的水面状态、观测角度、太阳光入射角度都会影响探测器接收信号的强弱，遥感图像覆盖范围广，不同地点上述因素差异又很大，不同像元受耀光污染的程度并不相同，因此需要对太阳耀光的影响因素进行分析，以确定太阳耀光的时空分布特征，为太阳耀光背景海面目标探测提供理论依据。

本发明通过天文学知识计算了太阳天顶角与方位角，通过给定的卫星轨道参数，计算了全球不同时间与地点的卫星观测天顶角与方位角，进而计算了卫星星下点轨迹分布。建立了波浪水面太阳耀光模型，进而通过大气辐射传输模型，模拟了大气层顶的太阳耀光分布模式，由此可以计算太阳耀光在全球出现的时间、地点与分布模式。

发明内容

本发明涉及一种在轨卫星成像波浪水面太阳耀光时空分布仿真模型，其具体步骤如下：

(1)太阳天顶角与方位角计算；

(2)卫星遥感器观测几何条件计算；

(3)波浪海面太阳耀光反射率分布建模；

(4)大气层顶太阳耀光反射率分布模式计算；

(5)全球太阳耀光出现的时间段与区域分布计算。

1步骤(1)中所述的“太阳天顶角与方位角计算”，其实现步骤如下：

(1)太阳天顶角计算

以地心为原点，建立三维坐标系，根据球面几何知识，太阳天顶角θ_s的余弦为

cosθ_s＝sinδsinφ+cosδcosφcost

其中，δ表示太阳赤纬，φ表示地面观测点纬度，t表示太阳时角，太阳赤纬依赖于计算日期对应的一年中的天数，根据天文学相关知识可以用如下的傅里叶级数近似计算：

δ＝0.006918-0.399912cos(A)+0.070257sin(A)-0.006758cos(2A)

+0.000907sin(2A)-0.002697cos(3A)+0.001480sin(3A)

其中，A＝2πJ/365，J是儒略历天数。

太阳时角的计算如下：首先根据格林威治标准时GMT计算经度为λ处地面点的平太阳时MST：

MST＝GMT+λ/15

之后将平太阳时修正为真太阳时(True Solar Time,TST)

TST＝MST+ET

其中，ET为真太阳时MST与平太阳时TST的差值，表达式如下：

$\begin{array}{c}ET=\frac{12}{\mathrm{\π}}(0.000075+0.001868\cos (A)-0.032077\sin \left(A\right)\\ -0.014615\cos \left(2A\right)-0.040849\sin \left(2A\right))\end{array}$

其中，A的定义同上。可以得到太阳时角t为

$t=15\frac{\mathrm{\π}}{180}(TST-12)$

这样就实现了太阳天顶角θ_s计算。

(2)太阳方位角计算

太阳方位角φ_s通过求解天顶(N)-太阳(S)-地面观测点(P)组成的球面三角形NSP计算得到，

$\frac{sin\mathrm{\χ}}{sin(\frac{\mathrm{\π}}{2}-\mathrm{\δ})}=\frac{\sin t}{{\mathrm{sin\θ}}_s}$

其中，χ是由南向西沿顺时针得到的太阳方位角。如果以正北方向为参考方向定义太阳方位角φ_s，可以得到

φ_s＝π+χ

这样计算得到的方位角取值范围为呈对称分布。可以根据太阳与地面观测点的经纬度关系将方位角范围变换为0～2π。

2步骤(2)中所述的“卫星遥感器观测几何条件计算”，其计算过程如下：

为了得到卫星传感器的观测天顶角与方位角，首先需要计算一定时刻卫星星下点轨迹的经纬度。星下点轨迹经纬度是根据目标卫星给定或假设的轨道六要素及时间信息解算得出的，计算中认为目标卫星绕地球作理想圆轨道运行，忽略了地球章动及星体摄动等因素对卫星轨道的影响，其计算过程如下：

目标卫星的轨道要素包括升交点赤经Ω、轨道面倾角i、卫星轨道长半轴a、卫星椭圆轨道离心率e、近地点角距ω、真近点角f。在上述轨道六要素中只有真近点角为关于时间的函数，必须明确某一时刻卫星的真近点角才能写出f(T)的函数，因此，假设卫星某次经过升交点的时刻。

根据轨道运动学定律，求解出卫星角速度W：

$W=\sqrt{\frac{G\·Me}{a^3}}$

其中G为万有引力常数，Me为地球质量，真近点角可表示为：

f(T)＝W·T

卫星在轨道直角坐标系中的位置可表示为：

$\left[\begin{array}{c}{x}_s\\ y_s\\ z_s\end{array}\right]=a\left[\begin{array}{c}cos\left(f\right(T\left)\right)\\ sin\left(f\right(T\left)\right)\\ 0\end{array}\right]$

轨道直角坐标系到协议天球坐标系的转换矩阵为：

$C_s^I=\left(\begin{array}{ccc}\cos \mathrm{\Ω}& -\sin \mathrm{\Ω}& 0\\ \sin \mathrm{\Ω}& \cos \mathrm{\Ω}& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos i& -\sin i\\ 0& \sin i& \cos i\end{array}\right)\left(\begin{array}{ccc}\cos \mathrm{\ω}& -\sin \mathrm{\ω}& 0\\ \sin \mathrm{\ω}& \cos \mathrm{\ω}& 0\\ 0& 0& 1\end{array}\right)$

再将卫星位置转换到地球直角坐标系下，协议天球坐标系到地球直角坐标系的转换矩阵为：

$C_I^t=\left(\begin{array}{ccc}cos\left(GAST\right)& sin\left(GAST\right)& 0\\ -sin\left(GAST\right)& cos\left(GAST\right)& 0\\ 0& 0& 1\end{array}\right)$

其中GAST为春分点格林威治恒星时，从而得到卫星在地球坐标系下的空间位置坐标：

$\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=C_I^t{C}_s^I\left[\begin{array}{c}{x}_s\\ y_s\\ z_s\end{array}\right]$

最后根据球面几何关系确定出卫星的星下点经纬度数据。

求出卫星星下点的经纬度后，可以进一步计算卫星的观测天顶角与方位角。设O为地心，S′为卫星星下点，P为地面观测点，卫星距地面高度为H，地球半径为R。由此可以计算OS′和OP′的夹角α，即卫星与地心连线、地面观测点与地心连线的夹角。

cosα＝sinδsinφ+cosδcosφcost

其中，δ、φ分别为卫星星下点和地面观测点的纬度，Δ为两者的经度差值。求出角α后解三角形OS'P可以得到角θ，根据

${\mathrm{\θ}}_v=asin\left(\right(1+\frac{H}{R}\left)sin\mathrm{\θ}\right)$

可以计算出卫星的观测天顶角θ_v。卫星观测方位角的计算与太阳方位角的计算相似。3步骤(3)中所述的“波浪海面太阳耀光分布建模”，其计算步骤如下：

波浪水面的太阳耀光反射率为

$\mathrm{\ρ}({\mathrm{\θ}}_s,{\mathrm{\θ}}_v,{\mathrm{\φ}}_s-{\mathrm{\φ}}_v)=\frac{\mathrm{\π}r}{4{\mathrm{cos\θ}}_s{\mathrm{cos\θ}}_v{\cos }^4\mathrm{\β}}P(z_x,z_y)$

其中，r为波浪面元菲涅尔反射率，β为波浪表面法线方向与天顶之间的夹角，计算如下：

$\sec \mathrm{\β}=\frac{\sqrt{2(1+{\mathrm{cos\θ}}_s\·{\mathrm{cos\θ}}_v+{\mathrm{sin\θ}}_s\·{\mathrm{sin\θ}}_v\·cos({\mathrm{\φ}}_s-{\mathrm{\φ}}_v))}}{{\mathrm{cos\θ}}_s+{\mathrm{cos\θ}}_v}$

波浪面元坡度概率分布P(z_x,z_y)的计算如下：

$P(Z_x,Z_y)=\frac{1}{2{\mathrm{\π\σ}}_x{\mathrm{\σ}}_y}\exp (-\frac{{\mathrm{\ξ}}^2+{\mathrm{\η}}^2}{2})\·f_{GC}$

$f_{GC}=\left\{\begin{array}{c}1-\frac{1}{2}{C}_{21}({\mathrm{\ξ}}^2-1)-\frac{1}{6}{C}_{03}({\mathrm{\η}}^2-3\mathrm{\η})+\frac{1}{24}{C}_{40}({\mathrm{\ξ}}^4-6{\mathrm{\ξ}}^2+3)\\ +\frac{1}{4}{C}_{22}({\mathrm{\ξ}}^2-1)({\mathrm{\η}}^2-1)+\frac{1}{24}{C}_{04}({\mathrm{\η}}^4-6{\mathrm{\η}}^2+3)\end{array}\right\}$

其中，σ_x、σ_y分别为沿风方向与垂直风向坡度方差：

$\left\{\begin{array}{c}{\mathrm{\σ}}_x^2=0.003+0.00192w\±0.002\\ {\mathrm{\σ}}_y^2=0.000+0.00316w\±0.004\\ C_{21}=0.01-0.0086w\±0.03\\ C_{03}=0.04-0.033w\±0.12\\ C_{40}=0.40\±0.23\\ C_{22}=0.12\±0.06\\ C_{04}=0.23\±0.41\end{array}\right.$

$z_x=\frac{{\mathrm{sin\θ}}_{v}sin({\mathrm{\φ}}_v-{\mathrm{\φ}}_s)}{{\mathrm{cos\θ}}_s+{\mathrm{cos\θ}}_v}$

$z_y=\frac{{\mathrm{sin\θ}}_{v}cos({\mathrm{\φ}}_v-{\mathrm{\φ}}_s)+{\mathrm{sin\θ}}_s}{{\mathrm{cos\θ}}_s+{\mathrm{cos\θ}}_v}$

$r=0.5\frac{{\sin }^2({\mathrm{\θ}}_a-{\mathrm{\θ}}_w)}{{\sin }^2({\mathrm{\θ}}_a+{\mathrm{\θ}}_w)}+0.5\frac{{\tan }^2({\mathrm{\θ}}_a-{\mathrm{\θ}}_w)}{{\tan }^2({\mathrm{\θ}}_a+{\mathrm{\θ}}_w)}$

其中，θ_a为水面入射光的入射角，θ_w为折射光的折射角，w为海面风速。

4步骤(4)中所述的“大气层顶太阳耀光分布模式计算”，具体计算过程如下：根据步骤(1)计算得到的太阳天顶角与方位角，步骤(2)计算得到的卫星观测天顶角与方位角，给定海面风速与风向，输入到步骤(3)建立的波浪水面太阳耀光反射率模型中，计算得到波浪水面太阳耀光的反射率，将其输入到大气辐射传输模型中，计算得到卫星高度经大气衰减后的海面太阳耀光分布模式。

5步骤(5)中所述的“计算全球太阳耀光出现的时间段与区域分布”，具体计算过程如下：根据步骤(4)计算得到的卫星平台太阳耀光分布模式的全球投影图中，叠加全球水陆分界线，剔除太阳耀光出现在陆地表面情形，即得到实际太阳耀光分布位置；设定不同的卫星成像时间，即可模拟在全球任意位置太阳耀光是否会出现，以及出现的空间位置。

本发明与现有技术相比的优点在于：

(1)本发明可以模拟任意卫星轨道平台成像的太阳耀光分布模式，适用于任意太阳同步轨道、地球同步轨道或其它倾斜轨道，只要给定卫星轨道6个参数，即可模拟任意时刻，全球范围在特定的观测几何条件下，海面区域是否会出现太阳耀光以及出现的位置与覆盖范围。即可模拟卫星在轨运行过程中，太阳耀光出现与消失的动态变化过程，可为海洋遥感监测提供新的技术平台。

(2)本发明得到的卫星平台太阳耀光分布模式基于精确的理论计算，得到的是定量的反射率，物理意义明确，可用于海洋水色遥感离水辐亮度或遥感反射率的反演中太阳耀光掩膜计算。本发明的精度高，速度快，可操作增强。

附图说明

图1为本发明的技术流程，图2为某一颗在轨卫星在2013年7月1日北京时间8:00、12:00与17:00三个时刻在全球60S—60N区域观测到的太阳耀光所在位置及分布模式。

具体实施方式

为了更好地说明本发明涉及的一种在轨卫星成像波浪水面太阳耀光时空分布模型，利用本发明的模型方法，进行了测试分析，取得了良好的效果，实验结果如图2所示，具体实施方法如下：

(1)根据给定的格林威治时间，计算太阳天顶角与方位角；

(2)根据卫星轨道参数，计算卫星遥感器观测几何条件；

(3)建立波浪海面太阳耀光分布模型，计算得到在特定时刻太阳耀光分布模式；

(4)耦合大气辐射传输模型，计算大气层顶太阳耀光分布模式；

(5)计算全球太阳耀光出现的时间段与区域分布。以2013年7月1日北京时间8时、12时与17时为例，计算了全球60S—60N区域水面太阳耀光的分布。其中，08:00、仿真图像中海面没有颜色表示不会观测到太阳耀光；12:00时图像中呈椭圆形区域表示太阳耀光的分布模式，其中S与Sat分别表示此时太阳与卫星所在位置；17:00图像中陆地出现等值线，表示计算得到的耀光位置，但实际地面为大陆，所以不会观测到耀光。

Claims (1)

1.一种在轨卫星成像波浪水面太阳耀光时空分布模型，其特征在于包含以下步骤：

(1)太阳天顶角与太阳方位角计算，太阳天顶角的具体计算过程如下：

以地心为原点，建立三维坐标系，根据球面几何知识，太阳天顶角θ_s的余弦为

cosθ_s＝sinδsinφ+cosδcosφcost

其中，δ表示太阳赤纬，φ表示地面观测点纬度，t表示太阳时角，太阳赤纬依赖于计算日期对应的一年中的天数，根据天文学相关知识可以用如下的傅里叶级数近似计算：

δ＝0.006918-0.399912cos(A)+0.070257sin(A)-0.006758cos(2A)

+0.000907sin(2A)-0.002697cos(3A)+0.001480sin(3A)

其中，A＝2πJ/365，J是儒略历天数；

太阳时角的计算如下：首先根据格林威治标准时GMT计算经度为λ处地面观测点的平太阳时MST：

MST＝GMT+λ/15

之后将平太阳时修正为真太阳时TST：

TST＝MST+ET

其中，ET为真太阳时MST与平太阳时TST的差值，表达式如下：

$\begin{array}{c}ET=\frac{12}{\mathrm{\π}}(0.000075+0.001868\cos \left(A\right)-0.032077\sin \left(A\right)\\ -0.014615\cos \left(2A\right)-0.040849\sin \left(2A\right))\end{array}$

其中，A的定义同上，可以得到太阳时角t为

$t=15\frac{\mathrm{\π}}{180}(TST-12)$

这样就实现了太阳天顶角θ_s计算；

太阳方位角的计算过程具体如下：

太阳方位角φ_s通过求解天顶N-太阳S-地面观测点P组成的球面三角形NSP计算得到，

$\frac{sin\mathrm{\χ}}{sin(\frac{\mathrm{\π}}{2}-\mathrm{\δ})}=\frac{\sin t}{{\mathrm{sin\θ}}_s}$

其中，χ是由南向西沿顺时针得到的太阳方位角，如果以正北方向为参考方向定义太阳方位角φ_s，可以得到

φ_s＝π+χ

这样计算得到的方位角取值范围为呈对称分布，可以根据太阳与地面观测点的经纬度关系将方位角范围变换为0～2π；

(2)卫星遥感器观测几何条件计算，具体计算过程如下：

为了得到卫星传感器的观测天顶角与方位角，首先需要计算一定时刻卫星星下点轨迹的经纬度；星下点轨迹经纬度是根据目标卫星给定或假设的轨道六要素及时间信息解算得出的，计算中认为目标卫星绕地球作理想圆轨道运行，忽略了地球章动及星体摄动等因素对卫星轨道的影响，其计算过程如下：

目标卫星的轨道要素包括升交点赤经Ω、轨道面倾角i、卫星轨道长半轴a、卫星椭圆轨道离心率e、近地点角距ω、真近点角f，在上述轨道六要素中只有真近点角f为关于时间的函数，必须明确某一时刻卫星的真近点角才能写出f(T)的函数，因此，假设卫星某次经过升交点的时刻；

根据轨道运动学定律，求解出卫星角速度W：

$W=\sqrt{\frac{G\·Me}{a^3}}$

其中G为万有引力常数，Me为地球质量，真近点角可表示为：

f(T)＝W·T

卫星在轨道直角坐标系中的位置可表示为：

$\left[\begin{array}{c}{x}_s\\ y_s\\ z_s\end{array}\right]=a\left[\begin{array}{c}cos\left(f\right(T\left)\right)\\ sin\left(f\right(T\left)\right)\\ 0\end{array}\right]$

轨道直角坐标系到协议天球坐标系的转换矩阵为：

$C_s^I=\left(\begin{array}{ccc}cos\mathrm{\Ω}& -sin\mathrm{\Ω}& 0\\ sin\mathrm{\Ω}& \cos \mathrm{\Ω}& 0\\ 0& 0& 1\end{array}\right)\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos i& -\sin i\\ 0& \sin i& \cos i\end{array}\right)\left(\begin{array}{ccc}cos\mathrm{\ω}& -sin\mathrm{\ω}& 0\\ sin\mathrm{\ω}& \cos \mathrm{\ω}& 0\\ 0& 0& 1\end{array}\right)$

再将卫星位置转换到地球直角坐标系下，协议天球坐标系到地球直角坐标系的转换矩阵为：

$C_I^t=\left(\begin{array}{ccc}\cos \left(GAST\right)& \sin \left(GAST\right)& 0\\ -\sin \left(GAST\right)& \cos \left(GAST\right)& 0\\ 0& 0& 1\end{array}\right)$

其中GAST为春分点格林威治恒星时，从而得到卫星在地球坐标系下的空间位置坐标：

$\left[\begin{array}{c}x\\ y\\ z\end{array}\right]=C_I^t{C}_s^I\left[\begin{array}{c}{x}_s\\ y_s\\ z_s\end{array}\right]$

最后根据球面几何关系确定出卫星的星下点经纬度数据；

求出卫星星下点的经纬度后，可以进一步计算卫星的观测天顶角与方位角；设O为地心，S′为卫星星下点，P为地面观测点，卫星距地面高度为H，地球半径为R；由此可以计算OS′和OP的夹角α，即卫星与地心连线、地面观测点与地心连线的夹角：

cosα＝sinδ′sinφ+cosδ′cosφcosΔ

其中，δ、φ分别为卫星星下点和地面观测点的纬度，Δ为两者的经度差值；求出角α后解三角形OS'P可以得到角θ，根据

${\mathrm{\θ}}_v=asin\left(\right(1+\frac{H}{R}\left)sin\mathrm{\θ}\right)$

可以计算出卫星的观测天顶角θ_v，卫星观测方位角的计算与太阳方位角的计算相似；(3)波浪海面太阳耀光分布建模，具体实现过程如下：

波浪水面的太阳耀光反射率为

$\mathrm{\ρ}({\mathrm{\θ}}_s,{\mathrm{\θ}}_v,{\mathrm{\φ}}_s-{\mathrm{\φ}}_v)=\frac{\mathrm{\π}r}{4{\mathrm{cos\θ}}_s{\mathrm{cos\θ}}_v{\cos }^4\mathrm{\β}}P(z_x,z_y)$

其中，r为波浪面元菲涅尔反射率，β为波浪表面法线方向与天顶之间的夹角，计算如下：

$\sec \mathrm{\β}=\frac{\sqrt{2(1+{\mathrm{cos\θ}}_s\·{\mathrm{cos\θ}}_v+{\mathrm{sin\θ}}_s\·{\mathrm{sin\θ}}_v\·cos({\mathrm{\φ}}_s-{\mathrm{\φ}}_v\left)\right)}}{{\mathrm{cos\θ}}_s+{\mathrm{cos\θ}}_v}$

波浪面元坡度概率分布P(z_x,z_y)的计算如下：

$P(Z_x,Z_y)=\frac{1}{2{\mathrm{\π\σ}}_x{\mathrm{\σ}}_y}\exp (-\frac{{\mathrm{\ξ}}^2+{\mathrm{\η}}^2}{2})\·f_{GC}$

$f_{GC}=\left\{\begin{array}{c}1-\frac{1}{2}{C}_{21}({\mathrm{\ξ}}^2-1)-\frac{1}{6}{C}_{03}({\mathrm{\η}}^2-3\mathrm{\η})+\frac{1}{24}{C}_{40}({\mathrm{\ξ}}^4-6{\mathrm{\ξ}}^2+3)\\ +\frac{1}{4}{C}_{22}({\mathrm{\ξ}}^2-1)({\mathrm{\η}}^2-1)+\frac{1}{24}{C}_{04}({\mathrm{\η}}^4-6{\mathrm{\η}}^2+3)\end{array}\right\}$

其中，σ_x、σ_y分别为沿风方向与垂直风向坡度方差，

$\left\{\begin{array}{c}{\mathrm{\σ}}_x^2=0.003+0.00192w\±0.002\\ {\mathrm{\σ}}_y^2=0.000+0.00316w\±0.004\\ C_{21}=0.01-0.0086w\±0.03\\ C_{03}=0.04-0.033w\±0.12\\ C_{40}=0.40\±0.23\\ C_{22}=0.12\±0.06\\ C_{04}=0.23\±0.41\end{array}\right.$

$z_x=\frac{{\mathrm{sin\θ}}_{v}sin({\mathrm{\φ}}_v-{\mathrm{\φ}}_s)}{{\mathrm{cos\θ}}_s+{\mathrm{cos\θ}}_v}$

$z_y=\frac{{\mathrm{sin\θ}}_{v}cos({\mathrm{\φ}}_v-{\mathrm{\φ}}_s)+{\mathrm{sin\θ}}_s}{{\mathrm{cos\θ}}_s+{\mathrm{cos\θ}}_v}$

$r=0.5\frac{{\sin }^2({\mathrm{\θ}}_a-{\mathrm{\θ}}_w)}{{\sin }^2({\mathrm{\θ}}_a+{\mathrm{\θ}}_w)}+0.5\frac{{\tan }^2({\mathrm{\θ}}_a-{\mathrm{\θ}}_w)}{{\tan }^2({\mathrm{\θ}}_a+{\mathrm{\θ}}_w)}$

其中，θ_a为水面入射光的入射角，θ_w为折射光的折射角，w为海面风速；

(4)大气层顶太阳耀光分布模式计算，具体计算过程如下：根据步骤(1)计算得到的太阳天顶角与方位角，步骤(2)计算得到的卫星观测天顶角与方位角，给定海面风速与风向，输入到步骤(3)建立的波浪水面太阳耀光反射率模型中，计算得到波浪水面太阳耀光的反射率，将其输入到大气辐射传输模型中，计算得到卫星高度经大气衰减后的海面太阳耀光分布模式；

(5)计算全球太阳耀光出现的时间段与区域分布计算，具体计算过程如下：根据步骤(4)计算得到的卫星平台太阳耀光分布模式的全球投影图中，叠加全球水陆分界线，剔除太阳耀光出现在陆地表面情形，即得到实际太阳耀光分布位置；设定不同的卫星成像时间，即可模拟在全球任意位置太阳耀光是否会出现，以及出现的空间位置。