CN104897130B - 天基遥感分块插值计算太阳高度角的方法 - Google Patents

Abstract

本发明涉及天基遥感分块插值计算太阳高度角的方法，包括确认遥感数据状态；以附带参数数据作为插值参照；遵循理论依据确定数据分块大小；确定数据插值模式；落实数据插值方程式；计算各项插值数据；计算各元数据块太阳高度角参数。适当增加太阳高度角数据，弥补了天基遥感太阳高度角数据不足；建立了确定天基遥感太阳高度角数据密度的方法和计算公式，利于确定合理的太阳高度角数据分布；充分提高天基遥感太阳高度角校正数据精度，利于遥感地物反演和识别；采用本技术方案，一级数据产品不需携带固定间隔点太阳高度角数据；技术方法具有易理解、易实现、实用佳和高精度等特点。

Description

天基遥感分块插值计算太阳高度角的方法

技术领域

本发明涉及天基遥感分块插值计算太阳高度角的方法，属于遥感技术领域。

技术背景

卫星、宇宙飞船和航天飞机等天基对地遥感获取数据过程中受到遥感平台和自然环境等多种因素的制约，在遥感反演之前必须对遥感数据进行定位、辐射和几何等方面的校正。太阳对地球的辐射中，太阳高度角具有日周期和年周期变化，对观测点的辐亮度影响甚大。太阳辐射类地球遥感以辐亮度反演地物则不利于对地物的识别，不符合定量遥感的要求。天基遥感太阳高度角校正有3个主流模式，即表观辐亮度、表观反射率和地表反射率等模式，主流模式在理论上都比较成熟，许多文献对其作了较为严谨的推介。太阳高度角校正实质上只是辐射角度问题，在主流模式上所起到的作用都一样。还有其它模式，如栗琳等学者用关系模型来研究太阳高度角校正。在遥感数据预处理中，太阳高度角校正是太阳辐射类地球遥感数据校正中的一项常态性处理工作。

遥感数据产品提供太阳高度角参数数据不足，遥感数据是大数据，若每个像元或局部像元附带一个太阳高度角数据或者定位信息，会相应增大数据产品数据量，增加数据通信负荷，增加数据存储空间和降低遥感应用时效。李先华等学者指出，中高空间分辨率卫星每景数据仅提供中心点太阳高度角数据，对于数据校正精度则十分勉强，并根据卫星传感器运行原理与定位关系叙述了卫星数据逐个像元的太阳高度角参数计算方法。该方法需掌握传感器在轨运行原理及其相关参数数据，对末端数据用户要求偏高。黄晓园、周靖斐等学者利用一级数据产品携带的固定间隔点太阳高度角数据进行插值，增加了太阳高度角数据，提高了太阳高度角校正精度。由于该方法插值函数与实际关系不一致，不利于充分提高太阳高度角校正精度。当前对太阳高度角数据的插值密度还没有相应研究，不利于对其数据插值密度的合理确定。没有携带充足太阳高度角数据的遥感数据产品则没有较好的实用方式解决太阳高度角精校正问题。

发明内容

本发明目的在于解决天基遥感数据产品太阳高度角数据不足而影响数据校正精度的现实问题，为了增加天基遥感太阳高度角数据而设计一种天基遥感分块插值计算太阳高度角的方法。

本发明为了实现上述目的，天基遥感分块插值计算太阳高度角的方法包括：确认遥感数据状态以便采取相应处理；利用附带参数数据作为数据插值的参照依据；按照数据精度要求、利用太阳高度角数据遵循理论依据确定数据分块大小；根据数据状态确定数据插值模式；根据插值模式落实数据插值方程式；根据插值方程式计算各项插值数据；利用插值数据计算各元数据块太阳高度角参数。

所述确认遥感数据状态以便采取相应处理，即遥感数据产品从校正数据状况可分为0～五级，确认数据校正级别，落实数据产品是否已经完成数据定位、辐射定标、系统几何校正和坐标映射，不同级别的数据产品需要根据附带定位数据分布、采集数据阵形和数据采集时序等状况采取相应的数据插值处理。

所述利用附带参数数据作为数据插值的参照依据，即天基遥感数据产品附带有相关参数数据，利用这些参数数据作为数据插值的参照依据。

所述按照数据精度要求、利用太阳高度角数据遵循理论依据确定数据分块大小，包括数据分块理论依据和数据分块操作方法。所述数据分块理论依据，即通过数据校正偏差率公式推论确立数据校正偏差率与太阳高度角偏差的关系，通过太阳高度角偏差特征规律确立数据块行列宽度与太阳高度角偏差的关系，以太阳高度角偏差特征规律和数据校正偏差率公式推论作为数据分块的理论依据。所述数据分块操作方法，包括数据分块方案和数据分块计算。所述数据分块方案，即对整幅或整景或整轨数据无法按圆形区域分块，为了满足数据精度要求，并使插值过程简单化，采取圆形区域最大内方块为元数据块。所述数据分块计算，即根据数据用途确定校正数据精度大小，调用附带太阳高度角参数数据，或根据公式计算太阳高度角参数，由数据校正偏差率公式计算出元数据块的太阳高度角偏差，依据数据块行列宽度与太阳高度角偏差的关系计算出元数据块行列宽度，确定元数据块大小。

所述根据数据状态确定数据插值模式，包括确定一级或二级数据产品元数据块经纬度插值模式；确定一级或二级数据产品元数据块采集时间插值模式；确定赤纬、时角参数插值模式。

所述根据插值模式落实数据插值方程式，包括根据相应的插值模式，落实一级数据产品元数据块经纬度插值方程式和采集时间插值方程式；或者根据相应的插值模式，落实二级数据产品元数据块经纬度插值方程式和采集时间插值方程式；根据赤纬、时角参数插值模式，落实赤纬、时角参数插值方程式。

所述根据插值方程式计算各项插值数据，即利用附带参数数据根据相应的插值方程式计算各项插值数据。

所述利用插值数据计算各元数据块太阳高度角参数，即根据太阳高度角公式，代入由插值得到的参数数据，计算各元数据块太阳高度角参数。

0级数据产品应完成一级数据产品相应校正处理后才可进行太阳高度角校正，可参照一级数据产品分块插值计算太阳高度角的方法；0级或一级数据产品完成二级数据产品相应的数据校正，可参照二级数据产品分块插值计算太阳高度角的方法；三级以后数据产品一般已进行了太阳高度角校正，可去除太阳高度角粗校正，添加太阳高度角精校正，可参照二级数据产品分块插值计算太阳高度角的方法。

太阳高度角与太阳天顶角之和为90°，对于太阳高度角校正，太阳高度角计算与太阳天顶角计算是等效的。

本发明只适用于太阳辐射类天基遥感增加太阳高度角数据的需要，应用于可见光、近红外等波段，地基、空基遥感由于观测区域较小，使用定位数据计算太阳高度角参数就已经满足数据校正精度要求，不需使用本技术方案；夜间对地遥感(使用红外等长波波段)和主动遥感(星载雷达如合成孔径雷达(SAR)使用微波等波段)属于非太阳辐射类遥感，不涉及太阳辐射因素，不需使用本技术方案。

在本发明中，其核心创新点在太阳高度角数据密度的确定上，因此，本发明还提供了一种确定天基遥感太阳高度角数据密度的方法，包括数据分块理论依据和数据分块操作方法；所述数据分块理论依据，即通过数据校正偏差率公式推论确立数据校正偏差率与太阳高度角偏差的关系，通过太阳高度角偏差特征规律确立数据块行列宽度与太阳高度角偏差的关系，以太阳高度角偏差特征规律和数据校正偏差率公式推论作为数据分块的理论依据；所述数据分块操作方法，包括数据分块方案和数据分块计算；所述数据分块方案，即对整幅或整景或整轨数据无法按圆形区域分块，为了满足数据精度要求，并使插值过程简单化，采取圆形区域最大内方块为元数据块；所述数据分块计算，即根据数据用途确定校正数据精度大小，调用或计算太阳高度角参数数据，由数据校正偏差率公式计算元数据块的太阳高度角偏差，依据数据块行列宽度与太阳高度角偏差的关系计算出元数据块行列宽度，确定元数据块大小；在遥感数据处理中元数据块大小确定，并且每个元数据块共同使用一个太阳高度角数据，也就确定了遥感数据的太阳高度角数据密度。

本发明得到一些有益效果：本发明以分块插值计算太阳高度角参数，为天基遥感增加太阳高度角数据，解决了天基遥感太阳高度角数据不足的实际问题；建立了确定天基遥感太阳高度角数据密度的方法和计算公式，利于确定合理的太阳高度角数据分布；充分提高天基遥感太阳高度角校正数据精度，利于遥感地物反演和识别；采用本技术方案，一级数据产品不需携带固定间隔点太阳高度角数据；采用插值法降低了末端数据用户增加太阳高度角数据的处理难度，技术方法具有易理解、易实现、实用佳和高精度等特点。

附图说明

图1为本发明实现流程图；

图2、图3分别为TERRA/MODIS、AQUA/MODIS卫星遥感一级数据产品某轨采集数据阵形示意图；

图4、图5分别为图2、图3对应的二级数据产品采集数据阵形示意图；

图6为HJ-1/CCD卫星遥感二级数据产品某景降轨采集数据阵形示意图。

具体实施方式

天基遥感分块插值计算太阳高度角的方法参见图1，是基于插值法以分块形式得到元数据块中心点经纬度值和采集时间，根据公式计算元数据块太阳高度角参数，其方法包括：确认遥感数据状态以便采取相应处理；利用附带参数数据作为数据插值的参照依据；按照数据精度要求、利用太阳高度角数据遵循理论依据确定数据分块大小；根据数据状态确定数据插值模式；根据插值模式落实数据插值方程式；根据插值方程式计算各项插值数据；利用插值数据计算各元数据块太阳高度角参数。选择IDL或FORTRAN等语言编写程序处理数据，先以太阳同步极轨卫星为例阐述如下：

⒈确认遥感数据状态以便采取相应处理遥感数据产品从校正数据状况可分为0～五级，太阳高度角校正适宜在一级或二级数据产品实施。0级数据产品未完成数据定位、辐射定标等处理，未达到本技术方案插值处理条件，一级和二级数据产品已完成数据定位、辐射定标，并且二级数据产品也完成了系统几何校正和坐标映射。太阳高度角校正可安排在系统几何校正之前或之后，应安排在几何精校正和地形辐射校正等之前。不同数据产品级别与数据插值模式有关。一级数据产品采集数据阵形与0级相同(参见图2、图3)，为方形，虽然完成了数据定位，由于未完成系统几何校正，数据定位信息在采集数据阵形上尚未得到体现。数据采集时序与0级相同，文件数据行与数据采集顺序一致。二级数据产品已完成系统几何校正，数据定位信息在采集数据阵形上已得到体现，并且由坐标映射建立了采集数据与坐标系的关系(参见图4、图5)。受地球自转、在轨道飞行方向及其速度和传感器推扫方向等影响，降轨时采集数据阵形扭摆西南向，升轨时采集数据阵形扭摆西北向，文件数据(含边空数据)顶角和图像数据(不含边空数据)顶角不一致，采集数据阵形接近菱形。数据采集行与文件数据行存在一定大小夹角，数据采集时序在降轨时为顺序，在升轨时为逆序。受地球曲率影响，幅宽较大的数据两侧出现校正拉伸。在系统几何校正中不同的投影方式对采集数据阵形的改变则有一些差别，等经纬度投影对高纬度区域拉伸尺度大于低纬度区域。幅宽较大的一级数据产品其太阳高度角校正往往与局地数据投影一起处理。不同级别的数据产品需要根据附带定位数据分布、采集数据阵形和数据采集时序等状况采取相应的数据插值处理。

⒉利用附带参数数据作为数据插值的参照依据天基遥感数据产品附带有相关参数数据，利用这些参数数据作为数据插值的参照依据，包括升降轨标记、数据空间分辨率L_m、数据采集日期D_ymd(或者采集年及积日)、数据采集开始时间t_s、采集结束时间t_e、数据总行数N_l和总列数N_s，以及一些关键点像元的定位、太阳高度角等参数数据。一级数据产品附带有固定间隔像元行数N_rl及列数N_rs、间隔点像元的定位和太阳高度角等参数数据；二级数据产品附带有整幅或整景或整轨数据中心点、四周顶角点像元的定位、太阳高度角等参数数据，如文件数据顶角经纬度，即左上角纬度经度(λ_vlu)，右上角纬度经度(λ_vru)，左下角纬度经度(λ_vld)，右下角纬度经度(λ_vrd)，以及整幅或整景或整轨数据中心点太阳高度角数据(θ_SEF)等，有的还附带有图像数据顶角参数数据，即左上角纬度经度(λ_lu)，右上角纬度经度(λ_ru)，左下角纬度经度(λ_ld)，右下角纬度经度(λ_rd)。在二进制数据文件中附带参数数据，常以整数表示，数据调用时要注意小数点移位，采用数据须符合国际单位要求，比如太阳高度角原数据已放大100倍，以整数表示，数据调用时要乘以0.01，数据复原后其单位为度。

⒊按照数据精度要求、利用太阳高度角数据遵循理论依据确定数据分块大小太阳高度角在较小的空间上和较短的时间上变化较小，在误差范围内根据校正数据精度可按一定数量的行宽列宽分块插值计算太阳高度角，同一分块使用同一个太阳高度角数据。数据分块大小是涉及插值密度的问题，是本技术方案的核心内容。

⑴数据分块理论依据以太阳高度角偏差特征规律和数据校正偏差率公式推论作为确定数据分块大小的理论依据。

①太阳高度角偏差特征规律太阳高度角偏差特征规律包括一般偏差特征规律和最大偏差特征规律。就一般偏差特征规律和最大偏差特征规律阐述如下：

㈠一般偏差特征规律太阳高度角计算可使用成熟通用公式

式中θ_SE为观测点太阳高度角，δ为观测日太阳赤纬，为观测点纬度，τ为观测点时角。规定赤道纬度为0°，北纬0～+90°，南纬0～-90°；格林尼治天文台经度为0°，西经0～-180°，东经0～+180°。

观测点与中心点太阳高度角偏差

θ_SEV＝θ_SE－θ_SE0 (2)

式中θ_SE0为中心点太阳高度角。

从空间几何关系分析，θ_SEV分布规律：对于同一时刻，太阳垂直照射点在观测区域内，以垂直照射点为圆心、圆心到中心点距离为半径，圆形区域内θ_SEV为正值，圆心点θ_SEV值最大，圆周上θ_SEV值为零，其它区域θ_SEV为负值，远离圆心点θ_SEV值较小；太阳垂直照射点在观测区域外，从中心点看，面向太阳这边为正值，正面对太阳方向值最大；背向太阳这边为负值，正背对太阳方向值最小；垂直照射点到中心点距离为半径，经过中心点的弧线上θ_SEV值为零；从中心点到正值负值区向外逐渐放大，出现周边效应；θ_SEV正值负值分布与观测时太阳所处位置有关，出现偏转现象。θ_SEV与θ_SE0关系：对于同一时刻，θ_SE0增大，其周围相对应点θ_SE增加，反之，θ_SE0减少，相对应点θ_SE减少。

㈡最大偏差特征规律太阳高度角偏差空间分布具有一定特征规律性。为了考察误差极值，分析太阳高度角最大偏差：

a)计算方法计算太阳高度角最大偏差采取如下方法：

I考察圆形区域数据整幅或整景或整轨遥感数据受地球自转、在轨道飞行方向及其速度等影响而导致观测区域为近似菱形，为适应太阳偏转，利于衡量太阳高度角最大偏差，以半径R_A(km)的圆形区域数据为考察对象。设太阳垂直照射在该圆形区域外(含圆边线)，根据太阳高度角偏差一般特征规律，方位角等于太阳方位角的圆形区域边上点其θ_SEV为正极值θ_SEV_M+，其镜像点θ_SEV为负极值θ_SEV_M-。太阳高度角最大偏差θ_SEV_M＝max(θ_SEV_M+，|θ_SEV_M-|)，式中max表示取最大值。

II计算太阳方位角跟踪太阳偏转需计算太阳方位角。太阳方位角计算可使用成熟通用公式

式中θ_A为太阳方位角，以正南为0°，上午为0～-180°，下午为0～180°。根据式(3)可计算某日某时刻中心点太阳方位角θ_A0。

III计算太阳方位角区域边上点经纬度求解方位角等于太阳方位角的圆形区域边上点经纬度，需计算地球表面两点之间的距离。设地球表面两点A(A_LON，A_LAT)和B(B_LON，B_LAT)，下标LON表示经度，下标LAT表示纬度。经纬度衡量尺度规定同上。设地球为正圆体，取地球平均半径R_E＝6371.004km，地球表面A、B两点距离公式

L_AB＝R_Earccos(C_M)π/180

C_M＝cosA_LATcosB_LATcos(A_LON－B_LON)+sinA_LATsinB_LAT (4)

式中L_AB(km)为A、B点距离，C_M为计算过程中间参数，π为圆周率。

设圆形区域中心点O(O_LON，O_LAT)和某日某时刻方位角等于太阳方位角的圆形区域边上点M(M_LON，M_LAT)及其镜像点N(N_LON，N_LAT)。在地球表面两点在纬度方向上距离与纬度差成正比，约111.193km等于1度纬度差宽；两点在经度方向上距离与经度差却不成正比，越接近两极其等距离的经度差则越大，越接近赤道其等距离的经度差则越小。求解方位角等于太阳方位角的圆形区域边上点及其镜像点经纬度，可根据几何关系由中心点纬度和方位角余弦函数关系计算纬度值，确定其纬度值后，由地球表面两点距离公式计算其经度值。计算边上点M(M_LON，M_LAT)及其镜像点N(N_LON，N_LAT)纬度和经度

M_LAT＝O_LAT－R_Acosθ_A0/111.193

N_LAT＝O_LAT+R_Acosθ_A0/111.193 (5)

M_LON＝O_LON+acos((cos(180R_A/(πR_E))－sinO_LATsinM_LAT)/(cosO_LATcosM_LAT))

N_LON＝O_LON－acos((cos(180R_A/(πR_E))－sinO_LATsinM_LAT)/(cosO_LATcosM_LAT))

或者M_LON＝O_LON－acos((cos(180R_A/(πR_E))－sinO_LATsinM_LAT)/(cosO_LATcosM_LAT))

N_LON＝O_LON+acos((cos(180R_A/(πR_E))－sinO_LATsinM_LAT)/(cosO_LATcosM_LAT)) (6)

已知O点经纬度值和M点纬度值，由地球表面两点距离公式计算M点经度值时得到两个经度值，刚好是目标点M及其对称点经度值(等于镜像点N点经度值)，可结合O点太阳方位角θ_A0能区分那个是M点、那个是N点经度值。

IV最大偏差计算已知某日某时刻及经纬度值，根据式(1)可计算中心点、方位角等于太阳方位角的圆形区域边上点、及其镜像点太阳高度角，进而计算得到太阳高度角偏差正极值θ_SEV_M+和负极值θ_SEV_M-。

b)特征规律据上分析，太阳高度角偏差θ_SEV_M+在太阳方位角方向上，θ_SEV_M-在太阳方位角反方向上。以太阳垂直照射在圆形区域外为例，经不同地域、不同季节、各时点大量资料验算，得到表1太阳高度角最大偏差θ_SEV_M与圆形区域半径R_A统计关系，有如下结论：

表1太阳高度角最大偏差与区域半径统计关系

半径RA/km	254.558	127.279	63.640	31.820	15.910	7.955
							最大偏差θSEVM/°	2.28924	1.14464	0.57232	0.28616	0.14308	0.07154

同一区域的太阳高度角偏差正极值θ_SEV_M+与负极值θ_SEV_M-绝对值相等；太阳高度角最大偏差θ_SEV_M与圆形区域半径R_A为线性关系，即

θ_SEV_M＝C₁R_A (7)

式中C₁＝8.993×10^-3(°/km)。虽然太阳高度角θ_SE与监测地经纬度、监测日期和监测时间有关，但是太阳高度角最大偏差θ_SEV_M与监测地经纬度、监测日期和监测时间为无相关关系，反映出监测地经纬度、监测日期和监测时间对于太阳高度角最大偏差θ_SEV_M不是实质性相关因子，与偏差圆形区域半径R_A才是实质性相关因子。

若太阳垂直照射在圆形区域里(不含圆边线)，设太阳垂直照射点A，其镜像点B，根据太阳高度角偏差一般特征规律，可得到太阳垂直照射点A与中心点O太阳高度之差θ_SEV_A＝θ_SEV_M+，并且镜像点B与中心点O太阳高度之差绝对值|θ_SEV_B|＝θ_SEV_A，由上述统计规律可得|θ_SEV_B|＜|θ_SEV_M-|，因此θ_SEV_M＝|θ_SEV_M-|，θ_SEV_M+＜θ_SEV_M。

②数据校正偏差率公式推论基于同质性像元假设，只考虑太阳高度角因素差异，由太阳高度角校正公式(注：可使用表观辐亮度公式或表观反射率公式，地表反射率公式因大气因素复杂化，不选用)套用偏差率公式，推导出太阳高度角校正偏差率

RV％＝(sinθ_SE0/sinθ_SE－1)×100 (8)

表明太阳高度角校正偏差率RV％与中心点太阳高度角θ_SE0、观测点太阳高度角θ_SE因子有关。由式(2)代入式(8)，得到

RV％＝(sinθ_SE0/sin(θ_SE0+θ_SEV)－1)×100 (9)

表明太阳高度角校正偏差率RV％与中心点太阳高度角θ_SE0、观测点太阳高度角偏差θ_SEV因子有关。经过计算比较，当θ_SEV一定时，θ_SE0减少，RV％增大，反之，θ_SE0增大，RV％减少。

一般所言数据精度满足应用需求，若评估数据误差极值，原则上所有点数据精度都要小于某一要求数值。太阳高度角校正最大偏差率

RV％_M±＝(sinθ_SE0/sin(θ_SE0±θ_SEV_M)－1)×100 (10)

太阳垂直照射的遥感数据区域是个例，其θ_SEV_M+＜θ_SEV_M，θ_SEV_M＝|θ_SEV_M-|，式中统一使用θ_SEV_M，可简化公式表达。对比发现，θ_SEV_M-对RV％_M影响大于θ_SEV_M+，在θ_SE0较大时，θ_SEV_M-和θ_SEV_M+对RV％_M影响较为接近；在θ_SE0较小时，θ_SEV_M-对RV％_M影响明显大于θ_SEV_M+，因此，数据校正最大偏差率公式优先采用偏差负极值θ_SEV_M-

RV％_M＝(sinθ_SE0/sin(θ_SE0－θ_SEV_M)－1)×100 (11)

对于数据分块，数据精度采用最大偏差率衡量，中心点太阳高度角采用附带参数最低的太阳高度角数据，由式(11)可确定数据块的太阳高度角最大偏差，由式(7)可计算出数据块圆形区域半径。可见，由上述推导的太阳高度角最大偏差特征规律及其校正最大偏差率公式推论能够确定数据分块大小。

⑵数据分块操作方法上述太阳高度角最大偏差θ_SEV_M与圆形区域半径R_A为线性关系是假设在某一时刻下推导的结论，地球遥感整幅或整景或整轨数据采集一般发生在一个较短的时间段内，虽然时间变化会影响太阳方位角改变，但是该结论不受时间变化影响，在遥感数据分析中仍然成立。可见，太阳高度角偏差特征规律和数据校正偏差率公式推论，成为了数据分块的理论依据。

①数据分块方案太阳高度角最大偏差特征规律是通过圆形区域考察得到的，但是，对整幅或整景或整轨数据进行几何分块却无法以圆形区域作为单元划分。为了保证分块区域小于该圆形区域，满足数据精度要求，并使插值过程简单化，采取圆形区域最大内方块为元数据块，块行列宽度L_p(km)，L_p＝2^1/2R_A。由式(7)推导得到太阳高度角最大偏差与方形数据块行列宽度的线性关系

θ_SEV_M＝C₂L_p (12)

式中C₂＝6.359×10^-3(°/km)。

②数据分块计算根据数据用途对校正数据精度的要求，判断数据有效位需要精确到哪一位数，确定校正数据精度大小。调用该幅或该景或该轨数据太阳高度角参数数据，或者根据公式计算四周顶角太阳高度角参数，若有多点参数数据则以最小值为准参与计算。由式(11)可计算出元数据块的太阳高度角最大偏差，由式(12)可计算出元数据块行列宽度L_p，L_p＝L_mN_p，式中L_m(km)为数据空间分辨率，N_p为元数据块像元行列数，可确定元数据块像元行列数N_p。为了插值计算简洁，取元数据块像元行列数为同步行数的倍数

N_p＝INT(L_p/(L_mN_t))N_t (13)

式中INT表示取整，仅取整数部分，小数不进位，N_t为数据采集的同步行数。框幅式传感器则需要特例处理，由于N_t＝N_l，导致违背设想初衷，按N_t＝1计算式(13)。若出现N_p＝0，小于空间分辨率则取N_p＝1；若N_p≥N_r，式中N_r＝max(N_rl，N_rs)，没有插值必要则终止插值操作。

⒋根据数据状态确定数据插值模式不同的数据状态决定了采用怎样的数据插值模式。

⑴确定一级或二级数据产品元数据块经纬度插值模式从太阳高度角校正对影响数据精度效果看，数据位置畸变尺度对于元数据块而言是小级别，可忽视数据位置畸变。一级数据产品未进行系统几何校正，数据位置畸变较大，然而，可参照邻近定位点经纬度插值就不影响数据位置插值精度；二级数据产品已进行系统几何校正，数据位置畸变较小，参照整幅或整景或整轨数据四周顶角或中心的定位点经纬度插值可满足数据位置插值精度要求。数据位置插值可采取参照点之间插值方式或者参照点外插值方式。宏观时空以等尺度量度，适宜使用线性插值算法。参照点之间插值方式以定位点位置为参照，参照点之间的经度、纬度以等尺度规律按经向、纬向线性插值；参照点外经度、纬度以等尺度规律按经向、纬向外推线性插值。遥感数据产品附带参数数据可满足参照点之间插值方式，比如地理纬度(或经度)线性插值：A_P＝A_A+(A_B－A_A)L_PA/L_BA，式中A_P为插值点P纬度(或经度)，A_A为参照点A纬度(或经度)，A_B为参照点B纬度(或经度)，L_PA为P点到A点的纬向(或经向)距离，L_BA为B点到A点的纬向(或经向)距离。

⑵确定一级或二级数据产品元数据块采集时间插值模式一幅或一景或一轨数据采集所需时间较短，并且太阳高度角在时空上缓慢变化，对采集时间插值没有必要区分到像元，可只区分到数据行。数据采集时间变化规律与传感器采集模式有关，各种传感器，例如框幅式传感器、逐点逐行扫描式传感器、逐点多行扫描式传感器、线阵列推扫式传感器、甚密面阵推扫式传感器和甚密线阵列扫描式传感器等采集机制各不相同，同时获取的同步行数也不一样，采集时间插值应根据不同模式传感器其同步行数、数据采集行的时序进行推算。宏观时空以等尺度量度，适宜使用线性插值算法。一级数据产品文件数据行与数据采集行时序仍然相同，可按照文件数据行推算，采集时间插值公式较简单；二级数据产品已完成系统几何校正和坐标映射，文件数据行与数据采集行已不再相同，存在一定大小夹角，即倾斜状态，降轨时采集数据阵形扭摆西南向，升轨时采集数据阵形扭摆西北向。数据采集时序在降轨时为顺序，在升轨时为逆序。按照倾斜的数据采集行及采集时序推算，采集时间插值公式相对较繁琐。倾斜状态可参照采集数据阵形，数据采集行斜率可由图像数据四周顶角定位数据计算，或者由文件数据的像元排列位置粗略计算，文件数据中各元数据块采集时间需要根据数据采集行斜率、采集开始时间和结束时间等数据插值计算。

⑶确定赤纬、时角参数插值模式受万有引力作用，卫星绕行星转，行星绕恒星转，恒星绕银河系转，其运行轨道表现出非常明显的简单规律性，但是，由于其它星球引力牵引，我们的星球运行轨迹也出现扰动现象，因此，也带来天文参数计算的复杂化。式(1)太阳高度角公式和式(3)太阳方位角公式在一般性应用的相关参数计算中较为多见，但是，太阳赤纬和时差等天文参数的计算，有数值模拟法和理论展开式法，有多种公式；时角参数计算引入时差订正，赤纬和时角插值则套用有关天文学者研究设计或改进的、精度较高又较为简洁的公式，其它相应的常见公式也可以套用。

⒌根据插值模式落实数据插值方程式根据太阳高度公式计算观测点太阳高度角参数，需要观测点纬度、观测日太阳赤纬和观测点时角，计算观测点时角需要观测点时间和观测点经度。确定了元数据块大小，即确定N_p值，可计算出元数据块行列数，当N_l－INT(N_l/N_p)N_p＝0，表示末行元数据块足行，元数据块行数为INT(N_l/N_p)+1，元数据块最大行号M_l＝INT(N_l/N_p)(注：从0号编起)，当N_l－INT(N_l/N_p)N_p＞0，表示末行元数据块不足行，元数据块行数为INT(N_l/N_p)+2，元数据块最大行号M_l＝INT(N_l/N_p)+1；当N_s－INT(N_s/N_p)N_p＝0，表示末列元数据块足列，元数据块行数为INT(N_s/N_p)+1，元数据块最大列号M_s＝INT(N_s/N_p)(注：从0号编起)，当N_s－INT(N_s/N_p)N_p＞0，表示末列元数据块不足列，元数据块列数为INT(N_s/N_p)+2，元数据块最大列号M_s＝INT(N_s/N_p)+1。设m、n分别为元数据块行列号，无符整数，0≤m≤M_l，0≤n≤M_s，需落实观测点纬度观测点经度λ_m,n＝λ(m,n)、观测点时间t_m,n＝t(m,n)、观测日太阳赤纬δ_m,n＝δ(m,n)、观测点时角τ_m,n＝τ(m,n)的具体形式。经纬度衡量尺度规定同上，以国际时计时，行增序方向与纬度增序方向相反。一级与二级数据产品的插值模式不同，其插值方程式也不一样。对一级、二级数据产品元数据块经纬度、采集时间和赤纬、时角的插值方程式介绍如下：

⑴一级数据产品元数据块经纬度、采集时间插值方程式元数据块中心点经纬度和采集时间插值方程式具体落实如下：

①元数据块中心点经纬度插值方程式地理经纬度可采取线性插值或双线性插值，以线性插值为例。根据一级数据产品相应的插值模式，元数据块中心点经纬度插值是依靠邻近定位点数据为参照。设固定间隔点像元行列数分别为N_rl、N_rs，第1个间隔点起始像元行列数分别为N_rl0、N_rs0，元数据块(m，n)中心点位于数据像元位置(N_p(m+0.5)，N_p(n+0.5))，设N_rlm＝N_p(m+0.5)－N_rl0，N_rsn＝N_p(n+0.5)－N_rs0，元数据块中心点与相邻定位点一般有4种状况：

㈠元数据块中心点与定位点重合，满足INT(N_rlm/N_rl)N_rl＝N_rlm和INT(N_rsn/N_rs)N_rs＝N_rsn条件，重合参照点C_{_m,n}位于定位点位置(N_rlm/N_rl，N_rsn/N_rs)，读取定位点数据，得到元数据块中心点经纬度

λ_m,n＝λ_{C_m,n} (14)

㈡元数据块中心点在两个同行相邻定位点之间，满足INT(N_rlm/N_rl)N_rl＝N_rlm和INT(N_rsn/N_rs)N_rs＜N_rsn条件，左参照点L_{_m,n}位于定位点位置(N_rlm/N_rl，INT(N_rsn/N_rs))，右参照点R_{_m,n}位于定位点位置(N_rlm/N_rl，INT(N_rsn/N_rs)+1)。读取定位点数据：L_{_m,n}纬度经度λ_{L_m,n}；R_{_m,n}纬度经度λ_{R_m,n}。根据邻近定位点位置和距离线性插值，得到元数据块中心点经纬度

λ_m,n＝λ_{L_m,n}+(λ_{R_m,n}－λ_{L_m,n})×(N_rsn－INT(N_rsn/N_rs)N_rs)/N_rs (15)

㈢元数据块中心点在两个同列相邻定位点之间，满足INT(N_rlm/N_rl)N_rl＜N_rlm和INT(N_rsn/N_rs)N_rs＝N_rsn条件，上参照点U_{_m,n}位于定位点位置(INT(N_rlm/N_rl)，N_rsn/N_rs)，下参照点D_{_m,n}位于定位点位置(INT(N_rlm/N_rl)+1，N_rsn/N_rs)。读取定位点数据：U_{_m,n}纬度经度λ_{U_m,n}；D_{_m,n}纬度经度λ_{D_m,n}。根据邻近定位点位置和距离线性插值，得到元数据块中心点经纬度

λ_m,n＝(λ_{D_m,n}+λ_{U_m,n})/2 (16)

㈣元数据块中心点与相邻定位点不同行不同列，满足INT(N_rlm/N_rl)N_rl＜N_rlm和INT(N_rsn/N_rs)N_rs＜N_rsn条件，左上参照点LU_{_m,n}位于定位点位置(INT(N_rlm/N_rl)，INT(N_rsn/N_rs))，右上参照点RU_{_m,n}位于定位点位置(INT(N_rlm/N_rl)，INT(N_rsn/N_rs)+1)，左下参照点LD_{_m,n}位于定位点位置(INT(N_rlm/N_rl)+1，INT(N_rsn/N_rs))，右下参照点RD_{_m,n}位于定位点位置(INT(N_rlm/N_rl)+1，INT(N_rsn/N_rs)+1)。读取定位点数据：LU_{_m,n}纬度经度λ_{LU_m,n}；RU_{_m,n}纬度经度λ_{RU_m,n}；LD_{_m,n}纬度经度λ_{LD_m,n}；RD_{_m,n}纬度经度λ_{RD_m,n}。根据邻近定位点位置和距离线性插值，得到元数据块中心点经纬度

λ_m,n＝(λ_{LU_m,n}+λ_{LD_m,n}+(λ_{RU_m,n}－λ_{LU_m,n}+λ_{RD_m,n}－λ_{LD_m,n})×

(N_rsn－INT(N_rsn/N_rs)N_rs)/N_rs)/2 (17)

特殊情况，在首行、末行、首列和末列等临近的元数据块中心点经纬度插值中，若元数据块中心点不在两个定位点之间，也不与某个定位点重合，则以邻近定位点经纬度值为参照插值。

②元数据块中心点采集时间插值方程式根据一级数据产品相应的插值模式，元数据块中心点采集时间插值按数据采集行推算，且数据采集行与文件数据行同序。设数据采集同步行数为N_t，对于框幅式传感器，t_e＝t_s，N_t＝N_l；对于逐行扫描式传感器，N_t＝1。N_t间隔时间t_t＝(t_e－t_s)/(N_l/N_t)，得到元数据块中心点采集时间

t_m,n＝t_s+(t_e－t_s)×(m+0.5)N_p/N_l (18)

框幅式传感器t_m,n＝t_s。

⑵二级数据产品元数据块经纬度、采集时间插值方程式根据二级数据产品相应的插值模式，元数据块中心点经纬度插值以整幅或整景或整轨数据四周顶角定位数据为依据，元数据块中心点采集时间插值按倾斜的数据采集行和采集时序推算。

①元数据块中心点经纬度插值方程式对元数据块中心点经度、纬度分别插值：

㈠元数据块中心点纬度：左边行密度右边行密度 在实际数据中，左边行密度与右边行密度一般略有差异，行密度取值为线性行密度

可简单取值为平均行密度，即得到元数据块中心点纬度

㈡元数据块中心点经度：上端列密度ρ_λu＝(λ_vru－λ_vlu)/N_s，下端列密度ρ_λd＝(λ_vrd－λ_vld)/N_s，在实际数据中，上端列密度与下端列密度也略有差异，列密度ρ_λ取值为线性列密度

ρ_λ＝(λ_vru－λ_vlu)/N_s+n(λ_vrd－λ_vld－λ_vru+λ_vlu)/(M_sN_s) (21)

ρ_λ可简单取值为平均列密度，即ρ_λ＝(λ_vru－λ_vlu+λ_vrd－λ_vld)/2N_s，得到元数据块中心点经度

λ_m,n＝λ_vlu+ρ_λ(n+0.5)N_p (22)

②元数据块中心点采集时间插值方程式落实采集时间插值方程式包括数据采集行斜率计算和按采集时序插值：

㈠数据采集行斜率计算数据采集行斜率可由整幅或整景或整轨四周顶角定位数据计算，或由文件数据的像元排列位置粗略计算:

a)顶角定位数据计算若附带有整幅或整景或整轨图像数据顶角参数数据，可计算图像数据采集行斜率：数据采集行上端斜率下端斜率 在实际数据中，数据采集行上端斜率与下端斜率一般略有差异，数据采集行斜率k_m,n取值为线性斜率

k_m,n＝k_u+m(k_d－k_u)/M_l (23)

k_m,n可简单取值为平均斜率，即k_m,n＝(k_u+k_d)/2。

b)像元排列位置计算调用数据文件中图像数据顶角像元排列位置，得到上左点位置(C_mlu，C_nlu)、上右点位置(C_mru，C_nru)、下左点位置(C_mld，C_nld)和下右点位置(C_mrd，C_nrd)。行增序方向与纬度增序方向相反，数据采集行上端斜率k_u＝(C_mlu－C_mru)/(C_nru－C_nlu)，下端斜率k_d＝(C_mld－C_mrd)/(C_nrd－C_nld)。在实际数据中，数据采集行上端斜率与下端斜率一般略有差异，数据采集行斜率k_m,n取值为线性斜率

k_m,n＝k_u+m(k_d－k_u)/M_l (24)

k_m,n可简单取值为平均斜率，即k_m,n＝(k_u+k_d)/2。

㈡按采集时序插值数据采集时序有顺序和逆序之分，顺序斜率为负，逆序斜率为正。根据数据采集时序选择顺序插值方式，或逆序插值方式：

a)采集时间顺序插值对文件数据处理，一般从其左侧上端开始，到右侧下端结束，由于数据采集行倾斜，即使采集时序是顺序，文件数据行列开始端和结束端并不对应于数据采集时间的开始和结束。求解元数据块中心点采集时间t_m,n，需要计算数据采集行斜率k_m,n、数据采集开始时间到结束时间所对应文件数据行数N_se、每行平均占用时间t_l、以及文件数据左侧上端开始时间t_vs。按顶角定位数据计算，倾斜列前增列数N_su＝N_s(λ_lu－λ_vlu)/(λ_vru－λ_vlu)，倾斜列后空列数N_sd＝N_s(λ_rd－λ_vld)/(λ_vrd－λ_vld)；按像元排列位置计算，N_su＝C_nlu，N_sd＝C_nrd－C_nld。数据采集开始时间t_s到结束时间t_e所对应文件数据行数

N_se＝N_l－k_uN_su+k_dN_sd (25)

每行平均占用时间

t_l＝(t_e－t_s)/N_se (26)

文件数据左侧上端开始时间

t_vs＝t_s－k_uN_sut_l (27)

元数据块中心点采集时间

t_m,n＝t_vs+t_l((m+0.5)+k_m,n(n+0.5))N_p (28)

b)采集时间逆序插值同上分析，由于数据采集行倾斜，即使采集时序是逆序，文件数据行列开始端和结束端也不对应于数据采集时间的结束和开始。求解元数据块中心点采集时间t_m,n，需要计算数据采集行斜率k_m,n、数据采集开始时间到结束时间所对应文件数据行数N_se、每行平均占用时间t_l、以及文件数据左侧上端结束时间t_ve。按顶角定位数据计算，倾斜列前空列数N_su＝N_s(λ_ru－λ_vlu)/(λ_vru－λ_vlu)，倾斜列后增列数N_sd＝N_s(λ_ld－λ_vld)/(λ_vrd－λ_vld)；按像元排列位置计算，N_su＝C_nru－C_nlu，N_sd＝C_nld。数据采集开始时间t_s到结束时间t_e所对应文件数据行数

N_se＝N_l－k_uN_su+k_dN_sd (29)

每行平均占用时间

t_l＝(t_e－t_s)/N_se (30)

文件数据左侧上端结束时间

t_ve＝t_e+k_uN_sut_l (31)

元数据块中心点采集时间

t_m,n＝t_ve－t_l((m+0.5)+k_m,n(n+0.5))N_p (32)

⑶赤纬、时角参数数据插值方程式赤纬和时角插值则套用有关天文学者研究设计或改进的、精度较高又较为简洁的公式，其它相应的常见公式也可以套用。以王炳忠等学者改进的太阳赤纬、时角等天文参数公式为例，这些公式融入了年度、经度和时刻等订正。

①赤纬参数数据插值方程式太阳赤纬是地球赤道平面与太阳和地球中心的连线之间的夹角，太阳赤纬公式

δ＝0.3723+23.2567sinθ+0.1149sin2θ－0.1712sin3θ－0.758cosθ+0.3656cos2θ+0.0201cos3θ

θ＝2π(J_d－(79.6764+0.2422(J_y－1985)－INT((J_y－1985)/4))+ΔJ)/365.2422 (33)

式中θ称日角(rad)，J_d为积日，积日为日期在年内的天数，由采集日期D_ymd计算，J_y表示监测年份，ΔJ为修正值(d)，ΔJ＝ΔJ₁+ΔJ₂，式中ΔJ₁为经度修正，ΔJ₂为时刻修正。ΔJ₁＝-0.00278λ_C，式中λ_C为监测地经度(°)，ΔJ₂＝T/24，式中T为监测地国际时时间(h)，若计时采用地方时，需根据时差将地方时转换为国际时。ΔJ值只需保留小数点后两位有效数。太阳赤纬δ即当地当日太阳直射地球的纬度，为年周期函数，在-23°26′～+23°26′范围内移动，从采集数据幅宽和采集时间等因素看，对于某幅或某景或某轨数据之中的各元数据块会使δ产生变化，但是，变化却甚小，可视作δ与元数据块行列无关，即δ_m,n＝δ。折中，取λ_C为整幅或整景或整轨数据中心点经度，取T为整幅或整景或整轨数据中心点国际时时间，计算修正值。

②时角参数数据插值方程式观测点时角即从观测点地球子午圈至太阳正对地球子午圈的夹角，以当地真太阳时计算，时角正午时为0°，每小时相差15°，午前为0～-180°，午后0～+180°，以元数据块参数表示，时角公式

τ_m,n＝15(t_m,n+Δt/60)+λ_m,n－180

Δt＝0.0028－1.9857sinθ+9.9059sin2θ–9.0924cosθ－0.6882cos2θ (34)

式中τ_m,n为元数据块时角，t_m,n为元数据块国际时采集时间(h)，λ_m,n为元数据块经度(°)，Δt为真平太阳时差(min)，θ按式(33)计算。

⒍根据插值方程式计算各项插值数据需要计算的各项插值数据有：λ_m,n、t_m,n、τ_m,n和δ，其中0≤m≤M_l，0≤n≤M_s。先计算δ，可一边对λ_m,n、t_m,n、τ_m,n插值，一边计算太阳高度角θ_SEm,n，或者全部对λ_m,n、t_m,n、τ_m,n完成插值后，再统一计算θ_SEm,n。

⑴边插值边计算编制程序代码处理数据时，建立λ、t、δ、τ变量，先计算δ，其它变量以元数据块行列号为控制，根据相应的插值方程式计算每个元数据块插值，赋值给λ_m,n赋值给λ、t_m,n赋值给t、τ_m,n赋值给τ，每个元数块插值完成，则计算相应θ_SEm,n；

⑵全部完成插值再计算建立δ变量和行为M_l+1、列为M_s+1的λ、t、τ数组，根据相应的插值方程式先计算δ，其它变量每个元数据块插值，赋值给数组、λ_m,n赋值给λ数组、t_m,n赋值给t数组、τ_m,n赋值给τ数组，统一完成插值再调用插值数组数据计算各项θ_SEm,n。

⒎利用插值数据计算各元数据块太阳高度角参数太阳高度角是指某地太阳光线与该地作垂直于地心的地表切线的夹角，其随着观测点、地方时和太阳赤纬的变化而变化。增加太阳高度角数据需要依靠插值数据支撑，根据太阳高度角公式，代入由插值得到的参数数据计算各元数据块太阳高度角参数。由太阳高度角公式，得到太阳高度角插值方程式

可一边插值λ_m,n、t_m,n、τ_m,n，一边计算相应θ_SEm,n，或者全部完成插值λ_m,n、t_m,n、τ_m,n后，再统一计算各项θ_SEm,n。建立行为M_l+1、列为M_s+1的θ_SE数组，将插值数据λ_m,n、t_m,n、τ_m,n和δ代入式(35)，计算每个元数据块太阳高度角θ_SEm,n，赋值给θ_SE数组，或者将θ_SE数组转存为数据文件，便于太阳高度角校正软件调用。

极轨遥感卫星考虑获取资料的需要，普遍采用太阳同步轨道。地球同步轨道卫星(简称静止卫星)由于其沿赤道绕地球转，并在赤道上空与地球保持静止状态，虽然轨道状态没有升轨和降轨的说辞，对于一级数据产品采集数据陈形和数据采集时序与极轨卫星一致；对于二级数据产品采集数据陈形及数据采集时序也是由地球自转、在轨道飞行方向及其飞行速度、传感器推扫方向等因素决定的，只是卫星在轨道抵消了地球自转因素，简化只看传感器推扫方向因素。采集数据陈形为方形，按传感器推扫方向不同，或倾斜或水平或垂直等状态，若传感器推扫方向从北向南，数据采集时序为顺序，若从南向北则为逆序；采集数据陈形为倾斜状态，采集时间与行、列都有关，为水平状态则只与行有关，与列无关，为垂直状态则只与列有关，与行无关。可见，一级数据产品元数据块经纬度插值及采集时间插值、二级数据产品元数据块经纬度插值等可直接参照极轨卫星方式的相应处理，只是二级数据产品元数据块采集时间插值则更复杂，不仅有扭摆西南向顺序插值和扭摆西北向逆序插值，还有扭摆东南向顺序插值、扭摆东北向逆序插值、水平顺序插值、水平逆序插值、垂直左向右推扫插值和垂直右向左推扫插值。扭摆西南向顺序插值和扭摆西北向逆序插值可直接参照极轨卫星方式的相应处理，其它形式的元数据块采集时间插值可借鉴极轨卫星方式的思路进行处理。宇宙飞船和航天飞机一级数据产品采集数据陈形和数据采集时序与极轨卫星一致；虽然在轨状态可能更复杂，简化归类，二级数据产品采集数据陈形和数据采集时序与静止卫星类似。形成怎样的采集数据陈形，主要在系统几何校正、坐标映射中考虑完成的工作，采集时间插值只需参考此前校正中得出的采集数据陈形结论。因此，宇宙飞船和航天飞机一级或二级数据产品元数据块经纬度插值、采集时间插值同样可参照极轨卫星、静止卫星等方式的相应处理。静止卫星、宇宙飞船和航天飞机一级或二级数据产品其它步骤可直接参照极轨卫星方式的相应处理。

以NOAA/AVHRR卫星遥感一级数据产品和HJ-1/CCD卫星遥感二级数据产品两个实施例作进一步说明。

实施例一：NOAA/AVHRR卫星遥感一级数据产品分块插值计算太阳高度角的方法

NOAA/AVHRR卫星遥感一级数据产品分块插值计算太阳高度角的方法，包括：

⒈确认遥感数据状态以便采取相应处理使用卫星数据是NOAA/AVHRR卫星遥感一级数据产品，已完成数据定位、辐射定标，适宜实施太阳高度角校正。采集数据阵形、数据采集时序与0级相同，文件数据行与数据采集顺序一致。数据产品为5通道1B数据格式的二进制数据文件，多为40～100MB，幅宽2800km。空间分辨率较低，整幅数据较宽，太阳高度角校正往往与局地数据投影一起处理。数据产品状态与数据插值模式有关。

⒉利用附带参数数据作为数据插值的参照依据数据产品附带相关参数数据有：升降轨标记U_d、数据空间分辨率L_m(星下点1.1km)、数据采集年及其积日、数据采集开始时间t_s、采集结束时间t_e、数据总行数N_l和总列数N_s(2048)、固定间隔点像元的定位和太阳高度角等参数数据。还附带有每行扫描线时间，每行安排51个固定间隔点，第一个间隔点在每行25列，间隔列数N_rs为40，间隔行数N_rl为1。参数数据附带在二进制数据文件记录头里，每个太阳天顶角数据占1个字节，前7位为整数，后1位为小数，该数据单位为1/2度；地理经纬度每个数据占2个字节.按纬度、经度顺序排列，前9位为整数，后7位为小数，该数据单位为1/128度。在调用这些附带参数数据时，应该准确读取并作相应处理。利用这些参数数据作为数据插值的参照依据。

⒊按照数据精度要求、利用太阳高度角数据遵循理论依据确定数据分块大小包括数据分块理论依据和数据分块操作方法。数据分块理论包括太阳高度角偏差特征规律和数据校正偏差率公式推论，数据校正偏差率公式推论确立数据校正偏差率与太阳高度角偏差的关系，太阳高度角偏差特征规律确立数据块行列宽度与太阳高度角偏差的关系；数据分块操作方法包括数据分块方案和数据分块计算。数据分块采取圆形区域最大内方块为元数据块，数据精度采用最大偏差率衡量，根据需要确定最大偏差率，比如确定最大偏差率为1％，采用附带参数最低的太阳高度角数据，由式(11)计算出元数据块的太阳高度角最大偏差，由式(12)计算出元数据块行列宽度L_p，由式(13)计算出元数据块像元行列数N_p。

⒋根据数据状态确定数据插值模式不同的数据状态决定了采用怎样的数据插值模式。⑴一级数据产品元数据块经纬度插值模式：一级数据产品未进行系统几何校正，数据位置畸变较大，可参照邻近定位点经纬度插值，就不影响数据位置插值精度；宏观时空以等尺度量度，适宜使用线性插值算法。参照点之间插值方式以定位点位置为参照，参照点之间的经度、纬度以等尺度规律按经向、纬向线性插值；遥感数据产品附带参数数据可满足参照点之间插值方式。⑵一级产品元数据块采集时间插值模式：一幅数据采集时间较短，采集时间插值只区分到数据行。使用逐行扫描式传感器，逐行推进。宏观时空以等尺度量度，适宜使用线性插值算法。采集时间插值应按照数据采集行的时序进行推算。一级数据产品文件数据行与数据采集行时序仍然相同，可按照文件数据行推算，采集时间插值公式较简单。⑶确定赤纬、时角参数插值模式：计算太阳赤纬、时差等天文参数有数值模拟法和理论展开式法，有多种公式，时角参数计算引入时差订正，赤纬和时角插值则套用有关天文学者研究设计或改进的、精度较高又较为简洁的公式，其它相应的常见公式也可以套用。

⒌根据插值模式落实数据插值方程式NOAA/AVHRR卫星遥感一级数据产品幅宽较宽，通常通过局地投影裁剪一定范围、需要用到的局部数据。在局地投影之前先安排太阳高度角校正，并借局地投影设置参数之机对投影范围使用到的固定间隔点定位、太阳高度角、采集开始时间、结束时间等参数数据从新整理，与投影区域数据建立对应关系，便于数据准确调用；在新设定的投影区域数据范围内对其像元、元数据块分别独立行列编号，以该区域像元行列号和元数据块行列号为控制进行数据处理。需要落实一级数据产品元数据块中心点经纬度、采集时间、赤纬、时角的插值方程式。元数据块中心点经纬度插值对于一级数据产品是依赖固定间隔点定位数据为参照，元数据块中心点与参照点如上所述一般有4种情况，但是，间隔点行数N_rl＝1，出现2种情况，其一为元数据块中心点与定位点重合，其二为元数据块中心点在两个同行定位点之间。元数据块中心点与定位点重合，按式(14)处理，元数据块中心点在两个同行定位点之间，按式(15)处理。一级数据产品数据采集时序与文件数据行一致，采集时间插值按式(18)处理，或者直接调用每行扫描时间。赤纬参数插值套用公式按式(33)处理，时角参数插值套用公式按式(34)处理。特殊情况，在首列和末列等临近的元数据块中心点经纬度插值中，若元数据块中心点不在两个定位点之间，也不与某个定位点重合，则以邻近定位点经纬度值为参照插值。

⒍根据插值方程式计算各项插值数据利用附带参数数据根据相应的插值方程式计算λ_m,n、t_m,n、τ_m,n和δ各项插值数据。

⒎利用插值数据计算各元数据块太阳高度角参数根据太阳高度角公式，代入由插值得到的参数数据，计算各元数据块太阳高度角参数。

实施例二：HJ-1/CCD卫星遥感二级数据产品分块插值计算太阳高度角的方法

HJ-1/CCD卫星遥感二级数据产品分块插值计算太阳高度角的方法,，包括：

⒈确认遥感数据状态以便采取相应处理使用卫星数据是HJ-1/CCD卫星数据二级数据产品，已完成数据定位、辐射定标、系统几何校正和坐标映射，适宜实施太阳高度角校正。受地球自转、在轨道飞行方向及其速度和传感器推扫方向等影响，二级数据产品降轨时采集数据阵形扭摆西南向(参见图6)，升轨时采集数据阵形扭摆西北向。数据采集行与文件数据行存在一定大小夹角，数据采集时序在降轨时为顺序，在升轨时为逆序。数据产品为2级4通道4个TIF数据文件，多为790～898MB，幅宽360km。数据产品状态与数据插值模式有关。

⒉利用附带参数数据作为数据插值的参照依据数据产品附带相关参数数据有：升降轨标记U_d、数据空间分辨率L_m(30km)、采集日期D_ymd、采集开始时间t_s、采集结束时间t_e、数据总行数N_l和总列数N_s，以及整景数据中心点、四周顶角点像元的定位、太阳高度角等参数数据，如文件数据顶角经纬度，即左上角纬度经度(λ_vlu)，右上角纬度经度(λ_vru)，左下角纬度经度(λ_vld)，右下角纬度经度(λ_vrd)，以及整幅数据中心点太阳高度角数据(θ_SEF)等，还附带有图像数据顶角参数数据，即左上角纬度经度(λ_lu)，右上角纬度经度(λ_ru)，左下角纬度经度(λ_ld)，右下角纬度经度(λ_rd)。参数数据记录在与TIF数据文件同名、后缀为XML的文件中。利用这些参数数据作为数据插值的参照依据。

⒊按照数据精度要求、利用太阳高度角数据遵循理论依据确定数据分块大小包括数据分块理论依据和数据分块操作方法。数据分块理论包括太阳高度角偏差特征规律和数据校正偏差率公式推论，数据校正偏差率公式推论确立数据校正偏差率与太阳高度角偏差的关系，太阳高度角偏差特征规律确立数据块行列宽度与太阳高度角偏差的关系；数据分块操作方法包括数据分块方案和数据分块计算。数据分块采取圆形区域最大内方块为元数据块，数据精度采用最大偏差率衡量，根据需要确定最大偏差率，比如确定最大偏差率为0.1％，采用附带景数据中心点太阳高度角数据，由式(11)计算出元数据块的太阳高度角最大偏差，由式(12)计算出元数据块行列宽度L_p，由式(13)计算出元数据块像元行列数N_p。

⒋根据数据状态确定数据插值模式不同的数据状态决定了采用怎样的数据插值模式。⑴二级数据产品元数据块经纬度插值模式：二级数据产品已进行系统几何校正，数据位置畸变较小，参照整景数据四周顶角或中心的定位点经纬度插值，可满足数据位置插值精度要求。宏观时空以等尺度量度，适宜使用线性插值算法。数据位置插值可采取参照点之间插值方式，以定位点位置为参照，参照点之间的经度、纬度以等尺度规律按经向、纬向线性插值。遥感数据产品附带参数数据可满足参照点之间插值方式。⑵二级产品元数据块采集时间插值模式：一景数据采集所需时间较短，对采集时间插值只区分到数据行。使用线阵列推扫式传感器，逐行推进。二级数据产品降轨时采集数据阵形扭摆西南向，升轨时采集数据阵形扭摆西北向。数据采集行与文件数据行存在一定大小夹角，数据采集时序在降轨时为顺序，在升轨时为逆序。宏观时空以等尺度量度，适宜使用线性插值算法。按照倾斜的数据采集行及采集时序推算，采集时间插值公式相对较繁琐。倾斜状态可参照采集数据阵形，数据采集行斜率可由四周顶角定位数据计算，或者由文件数据的像元排列位置粗略计算，文件数据中各元数据块采集时间需要根据数据采集行斜率、采集开始时间和结束时间等数据插值计算。⑶确定赤纬、时角参数插值模式：赤纬和时角插值则套用有关天文学者研究设计或改进的、精度较高又较为简洁的公式，其它相应的常见公式也可以套用。

⒌根据插值模式落实数据插值方程式根据二级数据产品相应的插值模式，元数据块中心点经纬度插值以整景数据四周顶角定位数据为依据，元数据块中心点采集时间插值按倾斜的数据采集行和采集时序推算。以降轨为例，受地球自转、在轨道飞行方向及其速度等影响，采集数据阵形扭摆西南向，采集时序为顺序。需要落实二级数据产品元数据块中心点经纬度、采集时间、赤纬、时角的插值方程式。求解元数据块中心点纬度是通过先计算纬度行密度，以开始行纬度为依据，参照式(20)处理；同理，求解元数据块中心点经度λ_m,n，先计算经度列密度，以开始列经度为依据，参照式(22)处理；求解元数据块中心点采集时间t_m,n，需要计算数据采集行斜率k_m,n、数据采集开始时间到结束时间所对应文件数据行数N_se、每行平均占用时间t_l、以及文件数据左侧上端开始时间t_vs，采集时序为顺序，计算k_m,n、N_se、t_l、t_vs分别参照式(23)、式(25)、式(26)、式(27)处理，计算t_m,n参照式(28)处理。赤纬参数插值套用公式按式(33)处理，时角参数插值套用公式按式(34)处理。升轨时，采集时序为逆序，与顺序插值类似，参照前述相应公式处理。

⒍根据插值方程式计算各项插值数据利用附带参数数据根据相应的插值方程式计算λ_m,n、t_m,n、τ_m,n和δ各项插值数据。

⒎利用插值数据计算各元数据块太阳高度角参数根据太阳高度角公式，代入由插值得到的参数数据，计算各元数据块太阳高度角参数。

卫星等天基遥感数据产品多为多波段(或称多通道)遥感数据，比如NOAA/AVHRR为5波段和HJ-1/CCD为4波段等，某幅或某景或某轨产品各个波段数据一般为同地域同幅宽同时采集，使用某幅或某景或某轨产品附带参数为插值依据进行分块插值计算得出的太阳高度角数据只适用于该幅或该景或该轨产品各个波段数据。

采用数据校正最大偏差率衡量数据精度则较为苛刻，也可采用平均偏差率等衡量。由数据校正平均偏差率公式确立数据校正平均偏差率与太阳高度角平均偏差的关系，根据太阳高度角偏差特征规律建立太阳高度角平均偏差与数据块行列宽度的关系，可确定元数据块大小。由数据校正平均偏差率建立的数据插值密度方法同样达到较好的数据校正效果。从工程实现来看，数据分块方案采用圆形区域最大内方块为单元分块最理想，其实，采用与圆形区域等面积方块为单元分块，或者采用圆形区域最小外方块为单元分块等，其效果同样很好，因为同样有圆形区域参照。类似这些更改都是在本技术方案中指导下实现的。

本文三角函数角的默认单位以度(°)表示，用计算机语言编程要注意三角函数角的单位，比如IDL语言三角函数角的默认单位是弧度(rad)，在使用中若为默认状态或弧度方式,需将角度值以弧度表示，式(4)中π/180、式(6)中180/π应相应去掉。

Claims (6)

1.一种确定天基遥感太阳高度角数据密度的方法，包括数据分块理论依据和数据分块操作方法；其特征在于：

以太阳高度角偏差特征规律和数据校正偏差率公式推论作为确定数据分块大小的理论依据；

所述数据分块理论依据，即通过数据校正偏差率公式推论确立数据校正偏差率与太阳高度角偏差的关系，通过太阳高度角偏差特征规律确立数据块行列宽度与太阳高度角偏差的关系，以太阳高度角偏差特征规律和数据校正偏差率公式推论作为数据分块的理论依据；具体为：

①太阳高度角偏差特征规律

一般偏差特征规律太阳高度角计算可使用成熟通用公式

式中θ_SE为观测点太阳高度角，δ为观测日太阳赤纬，为观测点纬度，τ为观测点时角；规定赤道纬度为0°，北纬0～+90°，南纬0～-90°；格林尼治天文台经度为0°，西经0～-180°，东经0～+180°；

观测点与中心点太阳高度角偏差

θ_SEV＝θ_SE－θ_SE0 (2)

式中θ_SE0为中心点太阳高度角；θ_SE为观测点太阳高度角；

为适应太阳偏转，利于衡量太阳高度角最大偏差，以半径R_A(km)的圆形区域数据为考察对象，太阳高度角最大偏差θ_SEV_M＝max(θ_SEV_M+，|θ_SEV_M-|)，式中max表示取最大值；

最大偏差计算

已知某日某时刻及经纬度值，根据式(1)可计算中心点、方位角等于太阳方位角的圆形区域边上点、及其镜像点太阳高度角，进而计算得到太阳高度角偏差正极值θ_SEV_M+与负极值θ_SEV_M-；

以太阳垂直照射在圆形区域外为例，经不同地域、不同季节、各时点大量资料验算，得到表1太阳高度角最大偏差θ_SEV_M与圆形区域半径R_A统计关系，有如下结论：

表1太阳高度角最大偏差与区域半径统计关系

同一区域的太阳高度角偏差正极值θ_SEV_M+与负极值θ_SEV_M-绝对值相等；太阳高度角最大偏差θ_SEV_M与圆形区域半径R_A为线性关系，即

θ_SEV_M＝C₁R_A (7)

式中C₁＝8.993×10^-3(°/km)；

若太阳垂直照射在圆形区域里，设太阳垂直照射点A，其镜像点B，根据太阳高度角偏差一般特征规律，可得到太阳垂直照射点A与中心点O太阳高度之差θ_SEV_A＝θ_SEV_M+，并且镜像点B与中心点O太阳高度之差绝对值|θ_SEV_B|＝θ_SEV_A，由上述统计规律可得|θ_SEV_B|＜|θ_SEV_M-|，因此θ_SEV_M＝|θ_SEV_M-|，θ_SEV_M+＜θ_SEV_M；

②数据校正偏差率公式推论

只考虑太阳高度角因素差异，由太阳高度角校正公式套用偏差率公式，推导出太阳高度角校正偏差率

RV％＝(sinθ_SE0/sinθ_SE－1)×100 (8)

表明太阳高度角校正偏差率RV％与中心点太阳高度角θ_SE0、观测点太阳高度角θ_SE因子有关；由式(2)代入式(8)，得到

RV％＝(sinθ_SE0/sin(θ_SE0+θ_SEV)－1)×100 (9)

表明太阳高度角校正偏差率RV％与中心点太阳高度角θ_SE0、观测点太阳高度角偏差θ_SEV因子有关；经过计算比较，当θ_SEV一定时，θ_SE0减少，RV％增大，反之，θ_SE0增大，RV％减少；

一般所言数据精度满足应用需求，若评估数据误差极值，原则上所有点数据精度都要小于某一要求数值；太阳高度角校正最大偏差率

RV％_M±＝(sinθ_SE0/sin(θ_SE0±θ_SEV_M)－1)×100 (10)

太阳垂直照射的遥感数据区域是个例，其θ_SEV_M+＜θ_SEV_M，θ_SEV_M＝|θ_SEV_M-|，式中统一使用θ_SEV_M，可简化公式表达；对比发现，θ_SEV_M-对RV％_M影响大于θ_SEV_M+，在θ_SE0较大时，θ_SEV_M-和θ_SEV_M+对RV％_M影响较为接近；在θ_SE0较小时，θ_SEV_M-对RV％_M影响明显大于θ_SEV_M+，因此，数据校正最大偏差率公式优先采用偏差负极值θ_SEV_M-

RV％_M＝(sinθ_SE0/sin(θ_SE0－θ_SEV_M)－1)×100 (11)

对于数据分块，数据精度采用最大偏差率衡量，中心点太阳高度角采用附带参数最低的太阳高度角数据，由式(11)可确定数据块的太阳高度角最大偏差，由式(7)可计算出数据块圆形区域半径；可见，由上述推导的太阳高度角最大偏差特征规律及其校正最大偏差率公式推论能够确定数据分块大小；

所述数据分块操作方法，包括数据分块方案和数据分块计算；

所述数据分块方案，即对整幅或整景或整轨数据无法按圆形区域分块，为了满足数据精度要求，并使插值过程简单化，采取圆形区域最大内方块为元数据块；元数据块的块行列宽度L_p(km)，L_p＝2^1/2R_A；由式(7)推导得到太阳高度角最大偏差与方形数据块行列宽度的线性关系

θ_SEV_M＝C₂L_p (12)

式中C₂＝6.359×10^-3(°/km)；

所述数据分块计算，即根据数据用途确定校正数据精度大小，调用或计算太阳高度角参数数据，由数据校正偏差率公式计算元数据块的太阳高度角偏差，依据数据块行列宽度与太阳高度角偏差的关系计算出元数据块行列宽度，确定元数据块大小；具体为：

由式(11)可计算出元数据块的太阳高度角最大偏差，由式(12)可计算出元数据块行列宽度L_p，L_p＝L_mN_p，式中L_m(km)为数据空间分辨率，N_p为元数据块像元行列数，可确定元数据块像元行列数N_p；为了插值计算简洁，取元数据块像元行列数为同步行数的倍数

N_p＝INT(L_p/(L_mN_t))N_t (13)

式中INT表示取整，仅取整数部分，小数不进位，N_t为数据采集的同步行数；

在遥感数据处理中元数据块大小确定，并且每个元数据块共同使用一个太阳高度角数据，也就确定了遥感数据的太阳高度角数据密度。

2.一种天基遥感分块插值计算太阳高度角的方法，其特征在于，是基于插值法以分块形式得到元数据块中心点经纬度值和采集时间，根据公式计算元数据块太阳高度角的方法，包括下列步骤：

(1)确认遥感数据状态以便采取相应处理；

(2)利用附带参数数据作为数据插值的参照依据；

(3)按照数据精度要求、利用太阳高度角数据遵循理论依据确定数据分块大小；

(4)根据数据状态确定数据插值模式；

(5)根据插值模式落实数据插值方程式；

(6)根据插值方程式计算各项插值数据；

(7)利用插值数据计算各元数据块太阳高度角参数；

在步骤(1)中，所述确认遥感数据状态以便采取相应处理，即遥感数据产品从校正数据状况可分为0～五级，确认数据校正级别，落实数据产品是否已经完成数据定位、辐射定标、系统几何校正和坐标映射，不同级别的数据产品需要根据附带定位数据分布、采集数据阵形和数据采集时序的状况采取相应的数据插值处理；

在步骤(2)中，所述利用附带参数数据作为数据插值的参照依据，即天基遥感数据产品附带有相关参数数据，利用这些参数数据作为数据插值的参照依据；

在步骤(3)中，所述按照数据精度要求、利用太阳高度角数据遵循理论依据确定数据分块大小，包括数据分块理论依据和数据分块操作方法；所述数据分块理论依据，即通过数据校正偏差率公式推论确立数据校正偏差率与太阳高度角偏差的关系，通过太阳高度角偏差特征规律确立数据块行列宽度与太阳高度角偏差的关系，以太阳高度角偏差特征规律和数据校正偏差率公式推论作为数据分块的理论依据；所述数据分块操作方法，包括数据分块方案和数据分块计算；所述数据分块方案，即对整幅或整景或整轨数据无法按圆形区域分块，为了满足数据精度要求，并使插值过程简单化，采取圆形区域最大内方块为元数据块；所述数据分块计算，即根据数据用途确定校正数据精度大小，调用或计算太阳高度角参数数据，由数据校正偏差率公式计算出元数据块的太阳高度角偏差，依据数据块行列宽度与太阳高度角偏差的关系计算出元数据块行列宽度，确定元数据块大小。

3.如权利要求2所述的天基遥感分块插值计算太阳高度角的方法，其特征在于，所述根据数据状态确定数据插值模式，包括：确定一级或二级数据产品元数据块经纬度插值模式；确定一级或二级数据产品元数据块采集时间插值模式；确定赤纬、时角参数插值模式。

4.如权利要求3所述的天基遥感分块插值计算太阳高度角的方法，其特征在于，所述根据插值模式落实数据插值方程式，包括：根据相应的插值模式，落实一级数据产品元数据块经纬度插值方程式和采集时间插值方程式；或者根据相应的插值模式，落实二级数据产品元数据块经纬度插值方程式和采集时间插值方程式；根据赤纬、时角参数插值模式，落实赤纬、时角参数插值方程式。

5.如权利要求4所述的天基遥感分块插值计算太阳高度角的方法，其特征在于，所述根据插值方程式计算各项插值数据，即利用附带参数数据根据相应的插值方程式计算各项插值数据。

6.如权利要求5所述的天基遥感分块插值计算太阳高度角的方法，其特征在于，所述利用插值数据计算各元数据块太阳高度角参数，即根据太阳高度角公式，代入由步骤(6)计算所得的各项插值数据，计算各元数据块太阳高度角参数。