ASTER多角度影像云参数模型解算与云迹风反演

张亚琳，高海博，何永健，邱新法，张国宏，王晓琼

(1.山西省气候中心，太原 030006；2.山西华冶勘测工程技术有限公司，太原 030002；3.南京信息工程大学地理科学学院，南京 210044；4.南京信息工程大学应用气象学院，南京 210044)

风场资料是天气和灾害预报研究中的重要资料，常规获取手段来自地面和高空观测气象站资料，但对于缺乏测站的地区(如高山、海洋、沙漠和极地等)获取风场资料难度较高，目前卫星资料反演风场成为重要的信息源[1]。卫星资料具有覆盖范围广、时空密度大、监测成本低、处理速度快、可动态监测等特点，能有效地弥补缺测地区观测信息不足的现状，利用卫星影像监测云系演变和提取云参数特征并展开各类研究，可为天气分析和数值预报模式提供大量有用资料。卫星云迹风(也称云导风)，即大气运动矢量(atmospheric motion vectors，AMVs)[2]，一般是指用连续几幅静止气象卫星云图追踪图像上示踪云或目标云的位移，并计算示踪云所代表的云或水汽特征所在的层次，以获得这些层次上风的估计值[3]。

气象卫星资料反演云迹风是从20世纪70年代开始的[4]，主要研究有早期的通过电影动画技术人工识别、模式识别法[5]、红外亮温相关系数法[6]等。近年来，学者在前述方法的基础上，通过新型算法或技术，改善了云迹风的算法及可视化显示效率。刘逸洁[7]运用葵花8号卫星实现了基于云导风技术的观测台风风场提取，基于红外、水汽双通道法提取云高的风场的修正，台风中心的自动提取和最大风速半径的估算，得到观测台风风场；何丽莉等[8]采用最大相关系数法等图像识别技术，实现了基于模糊聚类预处理云系的并行云迹风反演；张庆[9]利用长波红外辐射传输方程，系统地分析并完善了静止卫星红外云图晴空区微弱示踪信号提取算法——二阶差分法，为获取静止卫星晴空区低层大气运动提供了一种新方法，弥补了云导风和水汽导风的不足；贾静荣等[10]采用修正全变分(total variation，TV)范数光流法，实现了高分辨率云导风反演，依靠质量标识对所求导风进行质量控制，去掉不合理风矢，方法适用于400 hPa高度以下的中低云云导风计算；王昌帅等[11]提出一种基于OpenMP框架的云导风反演并行算法，对多组云图数据在多核CPU上进行高效反演；张军等[12]提出一种快速风场可视化算法，使用卫星云图序列直接生成便于人眼观察的大范围易于辨识的、时间一致性强的云导风动态可视化视频，能应对红外、可见光和水汽等各种波段的卫星图像数据。

立体图像由融合相同场景的左右视点图像形成[13]，双像对地立体观测是云迹风算法的基础。单星观测为同一卫星对目标连续观测两次，从而构成近似立体观测；双星联合立体观测指两颗静止卫星同时从不同角度建立立体观测；单星多角度观测是指极轨卫星上搭载的多个角度的对地观测仪器对目标进行准同时立体观测。极轨卫星影像较静止卫星精度更高，获取的云迹风风场精度更高。现将实现一种新的云迹风算法，把云顶风速作为一项未知云参数，选用分辨率较高的极轨Terra-ASTER单星多角度(3N、3B)影像，将数字摄影测量理论与遥感技术引入云参数解算中，通过构建立体像对云参数解算模型法和常规的立体像对云迹风反演法计算云顶风场，最后对获取的云顶风场进行检验，以期为天气分析和数值预报模式等研究提供基础资料。

1 研究资料

1.1 研究区概况

研究区位于西藏阿里地区普兰县附近的拉昂错湖和玛旁雍错湖以北，平均海拔4 000～5 000 m，多小型湖泊水体，并穿插有数条小山脉。由于其高海拔，该地区空气干燥、稀薄，太阳辐射较强，气温偏低，少降水。图1为2012年9月25日研究区所在ASTER影像，影像上显示有成片云团及少量积雪覆盖，蓝色圈出区域为实验云区，其面积约20×20 km2。

根据摄影测量技术，构建立体像对云顶参数解算模型，需先布设多组已知控制点进行模型的参数计算和检验。通过遥感软件ENVI，在ASTER 3N、3B像对中提取特征点像点坐标和其对应的空间三维坐标(物方坐标)，选定36对控制点为模型构建点，另选定16对加密点用于模型误差检验，布点如图1所示。

图1 研究区ASTER影像及控制点布点图

1.2 数据简介

ASTER(advanced spaceborne thermal emission and reflection radiometer)是Terra卫星上载有的5种对地观测仪器之一，全称“高级星载热发射反射辐射计”，包括从可见光到热红外共14个光谱通道，主要用来增进对地球表面或近地球表面和低层大气的了解，并着力解决土地利用与覆盖、自然灾害、短期天气变动、水文等方面的问题[14]。ASTER可见光波段3N、3B通过线性阵列相机同步同轨方式双角度独立成像，3N为星下点垂直成像波段，3B为观测天顶角27.6°的后视成像波段，两个相机对同一地区先后取景成像的时间间隔为55 s，满足双像高重叠度拍摄、具有一定基高比的立体成像必要条件。使用ASTER 3N、3B两波段的L1A级影像数据，空间分辨率15 m。

2 研究方法

通过两种方法获取云迹风风场：①云参数解算模型法：构建立体像对云顶参数解算模型，解算包括云顶风场、云三维参数在内的云参数；②云迹风反演法：提取双像同名像点对应的物方位置信息，结合双像成像时间间隔计算云迹风风场。两种方法的数据基础都是采用相关系数匹配法获得高精度、高密度云顶采样点[15]，将采样点作为已知数据代入模型解算云参数，或用于提取对应像点的位置信息。最后将两种方法的计算结果进行对比验证，并采用卫星影像云产品和探空数据对两种云迹风算法进行检验。

2.1 云参数解算模型

ASTER 3N和3B两波段成像影像均为推帚式线性阵列传感器扫描成像，即探测器每扫描一行，就构成一条中心投影的影像，如此随着传感器不断沿轨运动，从而构成连续影像。每个成像瞬间，影像摄影中心S、任意物点A和对应像点a三点共线[16]，根据摄影测量原理，3N影像共线方程为

(1)

式(1)中：

(2)

式中：(x,0)为第i行投影带像片坐标系中的像点；f为相机主距；(XA,YA,ZA)为对应地面点的空间坐标；前两者合称描述摄影中心与影像之间相关位置的内方位元素。外方位元素包含6个，(XSi,YSi,ZSi)为第i行的摄影中心S的空间坐标；航向倾角φ、旁向倾角ω及像片旋角κ是描述摄影光束空间姿态的角元素。

ASTER 3N、3B两波段为达到一定的图像重叠度，将3B相机以倾斜角θ=-27.6°(核线与Z轴夹角)进行成像扫描，结合式(1)得到3B影像的共线构像方程为

(3)

根据相关文献[17-18]，认为投影带行与行之间的外方位元素是线性相关的，任意行与中心行之间的一次线性关系方程式为

(4)

式(4)中：XS0、YS0、ZS0、φ0、ω0、κ0为中心行外方位元素；ΔXS、ΔYS、ΔZS、Δφ、Δω、Δκ为外方位元素变化率；XSi、YSi、ZSi、φi、ωi、κi为第i行的外方位元素；l0、li为中心行、i行的行号。

受北半球低纬度高空西风气流影响，云移动速度较快(最高能达到300 km/h以上)，对云参数解算精度影响较大，因此风须作为模型的一项重要参数(本文仅考虑水平风)。假设ASTER 3N上一云点的空间坐标为(XA,YA,ZA)，后视影像3B在3N成像ts(55 s)后对相同地面重新成像，设该云点以(Vx,Vy)的速度移动，则3B影像上同名云点空间坐标表示为(XA+Vxt,YA+Vyt,ZA)，代入式(3)，则考虑风速对云移动影响的3B影像共线构像方程为

(5)

方程的解算分为以下3步。

(1)后方交会解算模型的外方位元素及外方位元素变化率。对于3N和3B影像，均采用式(1)，通过代入地面控制点的物方和像点坐标、内方位元素等已知数据，解求3N和3B影像分别对应的外方位元素及外方位元素变化率。

(2)获取3N和3B影像摄影方向线和地面两个交点的距离。运用式(1)结合数字高程模型(digital elevation mode，DEM)，求得某一云点3N和3B影像摄影方向线和地面两个交点的水平距离(Sx,Sy)，其和云点移动的实际水平距离具有相关关系。

(3)前方交会解算未知云点的物方坐标。式(1)和式(5)代入已知参数和解求的参数，形成最终的云顶三维参数解算模型，代入匹配获取的双像同名云点像点坐标，求得云点对应的物方坐标(XA,YA,ZA)和云移动速度(Vx,Vy)。

2.2 云迹风反演

ASTER影像获取自极轨Terra-ASTER传感器，3N、3B两个通道的影像可组建立体像对实现云迹风反演。本文研究中所使用的是ASTER原始的L1A级数据，反演步骤如下。

2.2.1 ASTER L1A数据几何校正

ASTER影像具有较高的对地定位精度，可无控制点通过内置定位数据进行几何校正，定位后的影像在沿轨方向和垂直轨道方向的定位精度估计值分别为47 m和54 m[19]。

遥感软件ENVI中有专门针对ASTER L1A数据进行几何校正的工具。其中精度相对较低的方法为：根据ASTER头文件中包含的地理信息，通过Map-Georeference ASTER-Georeference Data对影像进行几何校正，该方法仅需用户自定义校正点数进行校正，处理速度快捷；使用IGM文件高精度校正：通过Map-Georeference from Input Geometry-Georeference from IGM进行校正，该方法处理速度较慢，但精度较高。

2.2.2 云迹风反演原理

对于经过几何校正预处理的两幅影像云图，云迹风的反演步骤一般为：示踪云的选取、示踪云的追踪(图像匹配)[20]、示踪云地理定位、云迹风风场(风速风向)推算。图2为云迹风反演原理示意图[1]，即在左影像3N上[图2(a)]确定示踪云A，在右影像3B搜索区C中[图2(b)]匹配同名目标云B，再求得A、B的X向移动距离ΔX和Y向移动距离ΔY，A、B直线移动距离为ΔS，结合两影像间隔的时间，即得到水平风矢。

图2 云迹风反演原理示意图

3 结果与分析

3.1 控制点对云参数解算模型的检验

为评价云参数解算模型的可用性，模型建立以后，利用选取的控制点对模型进行误差检验。将包括用于解算和加密检验在内的共52对控制点代入模型，得到控制点物方坐标，并同控制点真实坐标对比，求得模型计算控制点X、Y、Z方向的绝对误差，如图3所示。

由图3看出，三个误差分布没有明显的规律及相关性，但全都控制在-30～35 m，即2个像元之内，认为模型的稳定性较好。

图3 云参数解算模型解算控制点的X、Y、Z方向的绝对误差

3.2 云参数解算模型计算云顶风场结果

对云参数解算模型计算的云顶风场进行系统误差纠正后，分布示意图如图4所示，各色各等级风符号代表不同风速风向下的风矢，风速和风向的统计结果如图5所示。

图4 云参数解算模型计算的云顶风场分布示意图

图5 云参数解算模型计算的云顶风场的风速和风向统计结果

联合图4和图5得到，实验云区中云参数解算模型计算得到的风矢共2 885个，基本为西南偏北风向，实验区上半部分的风速明显高于下半部分。风速分布在8～24 m/s，其中12～20 m/s范围的风矢占90.23%，97.54%的风矢风向分布在0～25°。临近风矢间风速风向一致，同一云团的风速风向渐进变化。

3.3 云迹风反演结果

通过影像位置信息反演的云迹风风场分布示意图如图6所示，风速风向统计结果如图7所示。

从图6和图7看出，实验云区中反演的云迹风风场基本为西南偏北风向，符合青藏高原夏秋季热低压，且高空盛行西风气流的实际情况，实验区上半部分的风速明显高于下半部分。风速分布在8～24 m/s，其中12～20 m/s范围的风矢占92.13%，97.18%的风矢风向分布在0～25°，存在个别0°以下、30°以上异常风向的风矢。相邻风矢和同一云团的风矢间变化合理，与模型计算结果高度一致。

图6 云迹风反演的云顶风场分布示意图

图7 云迹风反演的云顶风场的风速和风向统计结果

3.4 结果检验

通过云参数解算模型计算的风场和云迹风反演获取的风场进行互检验，以及借助影像云产品与探空资料获得的风场规律佐证两个方法的可用性及可信性。

3.4.1 云参数解算模型计算风场与云迹风反演结果互检验

对比两种方法获取的风场结果，二者风速风向的分布范围、分布趋势、分布规律及变化规律都保持高度一致，初步证明了云参数解算模型计算风场与云迹风反演风场两种方法的可行性。

3.4.2 MISR云产品检验

MISR传感器与ASTER传感器一同搭载在Terra卫星上，二者对同一地区拍摄时间几近相同，二者有较高的可比性。MISR多角度影像二级产品提供了两种云数字产品，1999—2011年提供的云数字产品为MIL2TCST，提供分辨率为1.1 km的云高产品，包括不考虑风速影响的和考虑风速影响的两种云高，还提供了分辨率为70.4 km的风速产品，为6个角度影像计算风速中的最大风速；2012年，改进了云参数解算算法的MIL2TCSP产品，风速分辨率提高为17.6 km，为6个角度影像计算风速的平均风速，包含了南北方向和东西方向。对于云顶高度的验证，根据文献[21]进行，对于云移动风速的验证，采用了二代产品MIL2TCSP。实验云区大小约20×20 km2，仅占风产品的1～2个格点，实验结果与风产品风向一致，但二者风速不具有对比性。

3.4.3 探空资料检验

探空资料一般可作为实地验证的最佳对比资料，由于未能获得实验区同时段探空数据，且探空数据无法体现大面积风场情况，因此借助国家气象信息中心全国119个站无线电探空仪观测的1980—2006年共27年的月平均标量风速资料[22]作为参考，该资料是对风要素的气候平均处理，体现了中国12个月下不同等压面上高空平均风速常年的基本情况，具有稳定性和可信性。

云顶三维参数解算模型获取的实验区云顶海拔高度分布在6 723.4～9 148.4 m，处在300～500 hPa，结合上述探空资料，得到对应高度下同时段的高空月平均风速范围为12～23 m/s，而本文研究中两种云迹风风场计算方法分别有94.87%、96.24%的风速处在该范围，吻合度较高，证明了本文两种方法的可信性。

4 结论

云顶三维参数解算模型是遥感和摄影测量理论相结合计算云顶风场的新方法，而影像云迹风反演方法是近年来发展起来的监测高空风信息的一种有效手段，二者都可精细化监测高空风场信息。详细介绍了云顶三维参数解算模型和影像云迹风反演法获取风场的两种方法，对模型结果进行了控制点检验，得到误差保持在2个像元内，模型稳定性较好；对两种方法获取的风场进行互检验，并通过云产品、探空资料对两种方法的风场规律进行检验，其风速风向的分布范围、分布趋势、分布规律及变化规律都保持高度一致，验证效果较好，获取的风场质量好、密度高，清晰地反映出了高空的风场结构，证明了两种方法的可行性和可信性。不过，已有的风产品分辨率极低，为70.4 km和17.6 km，难以借助其验证本文计算结果，因此本文研究对风的验证无法达到实时精细化和定量化。未来考虑引用其他卫星多角度影像的云迹风计算结果进行验证分析。