本文共 5214 字，大约阅读时间需要 17 分钟。

IDL制作Himawari-8太阳角度数据集（太阳方位角、天顶角）

计算卫星方位角和方位角与计算太阳方位角的原理相同，根据《卫星轨道坐标系及坐标转换》（向汝建）https://max.book118.com/html/2015/0312/13226794.shtm一文中的总结，计算方位的步骤为以下几步: 计算太阳角度，首先计算儒略日JD，需要从HSD数据中读入observation start time

openr,h8_lun,inputfile,/get_lun;获取Observation start timepoint_lun,h8_lun,46imgtime=dblarr(1);R8readu,h8_lun,imgtime

由于HSD数据中读入的observation start time为MJD,又知MJD=JD-2400000.5，故而

jd = imgtime[0]*1.0d + 2400000.5

imgtime要时一个后面带d的double型数，乘以1.0d就可以使数据转为double类型，否则IDL会出现搞笑一幕：

如果是UTC时间，在IDL,Python或者matlab中都有将UTC时间转为儒略日时间的函数，可直接转换。 JD计算完成后是计算观测点的ECI坐标，即笛卡尔坐标系下的坐标。

计算ECI坐标，这里在之前的生成Himawari-8几何校正文件一文中https://blog.csdn.net/qq_33339770/article/details/102957857，对影像上每个像元的经纬度都进行过计算，故可形成像元坐标文件，这里的坐标是lon,lat,alt（可用ENVI自带DEM数据提取高程，也可直接定义为0）。 也可以自己定义左上角坐标80.0，60.0，行列数6001或12001，像元宽高0.02或0.01来书写lon，lat文件，但是最终计算结果会在纬度为0.0的一行结果反向，不想单独处理这一行，故而选择使用通过HSD文件中的行列偏移计算出的坐标文件。 需要将lon,lat,alt转为ECI坐标系下的x0,y0,z0坐标，根据论文中所述，地理经度在ECI坐标中是时间的函数，故而首先计算地方恒星时:地方恒星时=Greenwich恒星时+lon Greenwich恒星时的查找表比较麻烦，可以选择自己计算，计算的方法在论文中描述如下： 首先计算UTC时间为0：00的恒星时： ut是这一日的儒略日小数，也就是上图公式中的ΔT，由于儒略日的起始时间是公元前4713年1月1日中午12点，故而回溯半天（+0.5）再对1取余。 根据论文计算UTC 0：00的恒星时为（单位为秒）：

UT = (jd+0.5) mod 1.0jd_2000 = 2451545.0t_cen = (jd - UT - jd_2000) / 36525.0GMST = 24110.54841 + t_cen * (8640184.812866 + t_cen *(0.093104 - t_cen * 6.2E-6))

然后计算ΔT时间的恒星时，并将恒星时化为正值：

omega_e = 1.00273790934; Earth rotations per sidereal dayseconds_per_day = 86400.0GMST = ( GMST + seconds_per_day * omega_E * UT) MOD seconds_per_dayGMST[where(GMST LT 0.0)] += seconds_per_day

omega_e为地球每恒星日转过的弧度

之后将恒星时转为弧度：

rval = 2.0 * pi * GMST / seconds_per_day

整段计算恒星时的代码如下(此段theta_JD函数参考https://github.com/larrylart/Unimap/blob/db5c5b503d441bb784c7ed551741233ac7350e08/src/sky/sat_code/observe.cpp中的函数ThetaG):

;计算格林威治时间并转为角度Function thetaG_JD,jd  COMPILE_OPT idl2  pi=3.14159265358979323846  omega_e = 1.00273790934; Earth rotations per sidereal day    UT = (jd+0.5) mod 1.0  seconds_per_day = 86400.0  jd_2000 = 2451545.0  t_cen = (jd - UT - jd_2000) / 36525.0  GMST = 24110.54841 + t_cen * (8640184.812866 + t_cen *(0.093104 - t_cen * 6.2E-6))  GMST = ( GMST + seconds_per_day * omega_E * UT) MOD seconds_per_day  GMST[where(GMST LT 0.0)] += seconds_per_day  rval = 2.0 * pi * GMST / seconds_per_day  return,rvalEND

将地理经度带入公式，其后即可计算地方恒星时

d2r=0.01745329252;pi/180theta=(thetaG_JD(jd)+lon*d2r) MOD (2*pi)

有了地方恒星时，则可以将lon,lat,alt转为ECI坐标下的x0,y0,z0:

;经纬度转球面坐标(WGS_72椭球)Function ll2xyz_ellips,jd,lon,lat,alt  COMPILE_OPT idl2  pi=3.14159265358979323846  d2r=0.01745329252;pi/180  r2d=57.295779513;180/pi  Re=6378.14;赤道半径6378.140km(长半轴)  theta=(thetaG_JD(jd)+lon*d2r) MOD (2*pi)  ;转换到ECI坐标系下  f=1.0/298.256;地球扁率  c=1/sqrt(1+f*(f-2)*sin(lat*d2r)*sin(lat*d2r))  S=(1-f)*(1-f)*C  ;计算ECI坐标系下的坐标x,y,z  x=(Re+alt/1000.0)*C*cos(lat*d2r)*cos(theta)  y=(Re+alt/1000.0)*C*cos(lat*d2r)*sin(theta)  z=(Re+alt/1000.0)*S*sin(lat*d2r)  return,hash('x0',x,'y0',y,'z0',z)END

上述代码将经纬度坐标转为椭球体上的xyz坐标，也可以将地球视为球体，转为球面的坐标：

;经纬度转球面坐标(球体)Function ll2xyz_sphere,jd,lon,lat,alt  COMPILE_OPT idl2  pi=3.14159265358979323846  d2r=0.01745329252;pi/180  r2d=57.295779513;180/pi  Re=6378.14;赤道半径6378.140km(长半轴)  theta=(thetaG_JD(jd)+lon*d2r) MOD (2*pi)  x=Re*cos(lat*d2r)*cos(theta)  y=Re*cos(lat*d2r)*sin(theta)  z=Re*sin(lat*d2r)  return,hash('x0',x,'y0',y,'z0',z)END

对于ECI坐标计算存疑的，可以带入论文中的一套数据来验证计算出的x0,y0,z0是否正确,作者验证过椭球ECI坐标为正确的。

有了观测点的ECI坐标（x0,y0,z0）,同时可从HSD数据集中读到Sun’s position(ECI坐标)

;获取Sun's positionpoint_lun,h8_lun,510sun_pos=dblarr(3);R8readu,h8_lun,sun_pos

则可以计算观测点与太阳的向量差[rx，ry，rz]=[xs-x0，ys-y0，zs-z0]

这个向量差为ECI坐标系下的差向量，需要转为测站坐标系（地心坐标系）下的极坐标形式的向量差，就可以得到方位角、高度角（与天顶角互余）。 这里θ为天顶角（高度角的余角），φ为方位角 计算太阳方位角、天顶角的代码如下：

;太阳天顶角、方位角Function obslook_solar,jd,sun_pos,lon,lat,alt  COMPILE_OPT idl2  pi=3.14159265358979323846  d2r=0.01745329252;pi/180  r2d=57.295779513;180/pi  Re=6378.14;赤道半径6378.140km  ;计算经纬度对应的ECI坐标  xyz=ll2xyz_ellips(jd,lon,lat,alt)  theta=(thetaG_JD(jd)+lon*d2r) MOD (2*pi)  x0=xyz['x0']  y0=xyz['y0']  z0=xyz['z0']  OBJ_DESTROY,xyz  ;计算相对距离  xs=sun_pos[*,*,0]  ys=sun_pos[*,*,1]  zs=sun_pos[*,*,2]  ;卫星的位置是ECI坐标系，定义为[xs，ys，zs]，观察者是[x0，y0，z0]  ;那么距离向量就是[rx，ry，rz]=[xs-x0，ys-y0，zs-z0]  rx=xs-x0  ry=ys-y0  rz=zs-z0   ;[rx，ry，rz]是ECI坐标,要生成视角,需要转为地心地平线系统  ;这个系统的z轴指向天顶,x轴指向南,y轴指向西(地心坐标),地心坐标系随地球旋转  ls=sin(lat*d2r)*cos(theta)*rx+sin(lat*d2r)*sin(theta)*ry-cos(lat*d2r)*rz;南方,极坐标下：x=r*sin(el)*cos(az)  lse=-sin(theta)*rx+cos(theta)*ry;西方,极坐标下：y=r*sin(el)*sin(az)  lz=cos(lat*d2r)*cos(theta)*rx+cos(lat*d2r)*sin(theta)*ry+sin(lat*d2r)*rz;天顶,极坐标下：z=r*cos(el)  ;az:AzimuthAngle,方位角  az=atan(-lse/ls);atan(y/x)  az[where(ls gt 0)] += pi  az[where(az lt 0)] +=pi*2   ;相对距离  rg=sqrt(rx*rx+ry*ry+rz*rz)   ;el:ZenithAngle,天顶角  el=asin(lz/rg)  ;弧度转角度  el=el*r2d  az=az*r2d  ;输出SOA,SOZ  az[where(az GT 180.0)] -= 360.0  return,hash('SOA',az,'SOZ',90.0-el)END

对角度的计算验证也可以带入论文中的一组数据，这里作者未带入数据，直接拿结果与NC文件中的角度数据集对比正确。

代码传入数据说明，此文中的imgtime）,lon,lat,alt均与矫正裁剪后影像同等大小，为60016001或1200112001大小的数组，sun_pos需要存储三维的太阳坐标（xs,ys,zs）故数组大小为600160013或12001120013 imgtime数组的数值必须带d的double型数，包装成数组时可以采取下面的方法：

obstime=make_array(cols,rows,value=imgtime[0]*1.0d)

这样的话JD的计算就不需要*1.0d

jd = obstime+ 2400000.5

由于HSD每个波段分为10个分块，每个分块的Observation start time和Sun’s position都有变化，为准确的计算太阳角度，需要将10块imgtime和sun_pos读取之后作为此波段的time和sunpos数组，之后与影像一起进行几何校正和裁剪，即可作为角度数据集的输入数组。

转载地址：https://blog.csdn.net/qq_33339770/article/details/104710136 如侵犯您的版权，请留言回复原文章的地址，我们会给您删除此文章，给您带来不便请您谅解！