2019, Vol. 44引用本文 |
| 利用IceBridge测高数据计算海冰干舷的方法研究 |  [PDF全文] |
海冰是全球气候系统的重要组成部分，它在调节极地地区大气和海洋相互作用方面发挥着重要作用，获取准确的海冰参数，对于研究气候变化相关科学问题至关重要[1, 2]。“冰桥计划(Operation IceBridge, OIB)”是在地球两极进行的最大航空遥感科学观测工程，填补已停止运行的冰、云和陆地高程卫星(ice, cloud and land elevation satellite, ICESat)与2018年9月发射的ICESat-2之间数据间隙，延续激光测高仪时间序列。OIB已经获得了独特的时间序列，可用于分析南北极海冰的年际变化[3, 4]。利用机载或星载高度计提取海冰干舷和厚度的首要工作是获取瞬时海面高[5-7]，因此，获取准确的海面高对于计算海冰干舷、厚度至关重要[8]。
利用OIB获取海面高的方法主要有目视选取法、最低点/阈值法、表面反射率法等。目视选取法是利用机载地形扫描测绘仪1B级测高数据(airborne topographic mapper level-1B elevation and return strength, ATM L1B)和同步的数字测图系统相机1B级影像(digital mapping system level-1B imagery, DMS L1B)目视选取冰间水道(开阔水面、薄冰等)内较低高程点计算局部海面高[9]，这种方法精度较高，但人工选取工作量较大，其主要用于对比验证；Wang等[9]利用高程最低阈值0.2%计算局部海面高，该方法需要DMS影像辅助目视识别是否含有冰间水道；Kwok等[10]将表面反射率作为激光脚点是否在冰间水道内的指标，把高反射率海冰表面与低反射率冰间水道内激光脚点分离开来，实验中发现由于无雪覆盖较厚的冰面与少量雪面上激光具有相近的表面反射率，仅利用该方法很难分离出冰间水道内高程较低的激光脚点。
为了解决海冰干舷计算中目视法、阈值法、表面反射率法的局限性。本文提出了利用表面反射率法和高程最低点结合的方法。首先根据实验确定的冰间水道内激光脚点表面反射率截止值(Rcutoff)，将部分脚点分离出来; 然后再根据开阔水面、薄冰、较厚的冰面高程差异，在沿轨1 km分段取高程最低的10个点来计算局部海面高，以实现程序化快速、准确地计算海冰干舷的目的。
1 海冰干舷计算的方法流程利用表面反射率法和高程最低点结合的方法计算海冰干舷的流程如图 1所示，为达到仅利用ATM L1B测高数据程序化快速准确计算海冰干舷的目的，首先需要确定用于提取冰间水道内ATM L1B脚点的Rcutoff。利用Rcutoff提取出冰间水道内ATM L1B脚点，再使用高程最低点法获取局部海面高，用于计算海冰干舷。
 |
| 图 1 海冰干舷计算流程图 Fig.1 Flow Chart of Sea Ice Freeboard Retrieval |
1.1 确定表面反射率截止值
表面反射率(R)是ATM接收激光信号强度与发射激光信号强度之比。确定Rcutoff是实现本文提出的利用ATM L1B数据快速准确计算海冰干舷的基础。由于南极和北极海冰表面状况的差异，以及2017年OIB对相关仪器进行了更新。因此有必要对新的表面反射率截止值进行研究[11]。为了能够更加准确地获取Rcutoff，利用DMS L1B影像目视选取冰间水道中ATM L1B样本点，并分析样本点表面反射率的概率密度分布。图 2(a)是DMS L1B影像，可以清晰地识别出冰间水道和海冰，图 2(b)和图 2(c)分别是同步获取的ATM L1B点的表面反射率和高程分布，ATM L1B点颜色从蓝到红，代表的表面反射率和高程值逐渐增大，从中可以看出表面反射率分布和高程分布之间存在较好的一致性，并且冰间水道内脚点与海冰上脚点的表面反射率和高程均存在明显差异，说明了利用Rcutoff提取冰间水道内的脚点，再使用高程最低点法计算海面高的方法的可行性。
 |
| 图 2 沿轨DMS L1B影像和ATM L1B数据高程、表面反射率分布 Fig.2 DMS L1B Image and the Spatial Distribution of ATM L1B Data Elevation and Apparent Reflectivity Along the Fight Track |
图 2(d)是图 2(a)影像区域目视选取的冰间水道内ATM L1B点，用于统计分析确定冰间水道内ATM L1B点的Rcutoff。图 3是目视选取的冰间水道内89 218个ATM L1B样本点的概率密度分布图，参照Kowk等[10]方法，取冰间水道内样本ATM L1B点概率密度曲线峰值对应的表面反射率R=0.33作为Rcutoff，用于后续ATM L1B数据处理中替代目视选取冰间水道中ATM L1B点，以实现程序化快速处理获取海冰干舷的要求。
 |
| 图 3 冰间水道内ATM L1B点的表面反射率概率密度分布 Fig.3 Probability Density Distribution of ATM L1B Samples from the ATM L1B Points over Leads, Vertical Line Shows Rcutoff at 0.33 |
1.2 冰间水道探测方法
利用DMS L1B影像是冰间水道识别最准确、最直接的方法，但是面对海量数据，需要巨大的人力，且效率低下。为了能够快速准确地处理ATM L1B数据，利用上述实验得到的Rcutoff来判断ATM沿轨1 km分段是否含有冰间水道，如果该分段上脚点的表面反射率在0~Rcutoff之间数目超过10个，则认为该1 km分段含有冰间水道。另外，如果某个1 km分段或部分分段与相邻分段获取的海面高相差较大，则可能该分段并没有冰间水道，这时需要利用DMS L1B影像加以验证。
1.3 局部海面高的计算方法实验中发现海冰表面状态复杂，少量雪面或光滑海冰上与冰间水道内脚点表面反射率相近，仅仅利用Rcutoff很难将冰间水道内较低ATM L1B点单独提取出来。较低的表面反射率和高程值是判断脚点是否位于冰间水道的有效指标[10]。参照Wang等[9]利用目视法在1幅DMS影像(沿轨约600 m)上选取海面高连接点数目3~9个的原则，利用Rcutoff判断每个1 km分段是否含有冰间水道，如果有冰间水道，提取反射率在0~0.33之间的脚点，取高程最低10个ATM L1B点，按1倍中误差剔除后，剩余ATM L1B点高程均值作为局部海面高；实验段A、B、C、D中没有冰间水道的分段，根据每个实验段中含有冰间水道的1 km分段与沿轨距离(实验段第1个采样点为起点)之间线性拟合方程插值得到。
1.4 海冰干舷的计算方法海冰干舷通过计算冰面高程与邻近冰间水道高程之差得到[9, 12, 13]。ATM高程数据是相对于WGS-84参考椭球高，而实际海面高受潮汐、洋流、大气荷载的影响[11]。某一点的瞬时海面高hssh[14]：
| $ {h_{{\rm{ssh}}}} = {h_{{\rm{geoid}}}} + {h_{{\rm{tides}}}} + {h_{{\rm{pressure}}}} + {h_{{\rm{dynamic}}}} $ | (1) |
式中，hgeoid为大地水准面高; htides为潮汐影响; hpressure为大气荷载影响; hdynamic为海面动态地形。hgeoid、htides、hpressure这3个参数分别来源于地球重力场模型EGM2008、海潮模型TPXO6.2、ECMWF Interim reanalysis数据。对于海冰干舷计算1 km的分段尺度上，脚点间隔3~4 m，目前数值模拟值尚不能满足所需的cm级精度和空间分辨率要求[9]。因此，在每个1 km分段ATM L1B高程与该分段的hssh求差即可得到海冰干舷[7, 9, 15]。
沿轨1 km分段海冰干舷计算原理，如图 4。
| $ F = {h_{\rm{s}}} - {h_{{\rm{ssh}}}} $ | (2) |
 |
| 图 4 雪-海冰模型以及参数的定义 Fig.4 Snow-Sea Ice Model and Definition of Parameters |
式中，F表示海冰总干舷；为海冰表面雪厚(Fs)和冰干舷(Ff)之和；hs表示雪面上高程；hssh表示瞬时海面高；Ti表示海冰厚度。
计算出每个1 km分段海冰干舷之后，将每个分段分别汇总到对应的实验段A、B、C、D即可得到每个实验段沿轨海冰干舷。
2 海面高和海冰干舷分析 2.1 ATM L1B与DMS L1B数据描述实验所用的ATM L1B数据和DMS L1B影像，均来源于美国冰雪数据研究中心(https://nsidc.org/icebridge/portal/)。
ATM是一种锥形扫描激光测距系统，工作波长为532 nm，脉冲重复频率为5 kHz，扫描速率为20 Hz, 最低雷达扫描角度15°(T2扫描仪)、23°(T3扫描仪)或30°(T4扫描仪)[9, 16, 17]。为了提供准确的高程估计，该仪器包括用于精确定位飞机的差分GPS系统，以及用于精确姿态确定的惯性导航系统(加速度计和陀螺仪)[10, 18]。利用从发射到接收激光脉冲的时间差、飞机的位置和姿态信息来测量表面高度[19]。
ATM L1B数据提供沿轨条带扫描高程、发射/接收激光强度、姿态等测量数据，参考框架是ITRF-2005，参考椭球是WGS-84椭球[11]。标称飞行高度(500 m)，跨轨扫描条带宽度约250 m，实际扫描宽度随飞行高度变化而改变，激光脚点大小约1 m，条带中心附近激光脚点间距约3~4 m，由于ATM锥形扫描，在条带边缘激光脚点间距亚m级，沿着飞行路径存在极其光滑的表面(开阔水面、无雪覆盖光滑海冰)，激光的镜面反射会导致部分测量点丢失[10]。Martin等[20]估算的ATM系统水平准确度74 cm，垂直准确度6.6 cm，垂直精度3 cm。
DMS是一种数字摄影测量系统，获取中低空高分辨率的彩色或全色波段影像，DMS L1B影像是DMS原始影像经地理位置和正射校正后的影像，其空间分辨率随飞行高度而变化，空间分辨率范围为0.015 m~2.5 m[21]，主要用于识别冰间水道，通过提取冰间水道内高程较低脚点计算局部海面高[9]。
2.2 实验区域概况实验区域位于南极威德尔海域(如图 5)，选取OIB_20171114航线中4个ATM L1B实验段A、B、C、D，利用表面反射率和高程最低点相结合的手段计算局部海面高，实验段的基本信息如表 1所示。
 |
| 图 5 OIB_20171114海冰航线和实验段A、B、C、D的位置 Fig.5 Operation IceBridge Sea Ice Flight in Weddell Sea of Antarctica on November 14, 2017 and Four Selected Study Sections of A, B, C and D |
| 表 1 实验段A、B、C、D基本信息 Tab.1 Basic Information of Sections A, B, C and D |
 |
2.3 局部海面高的精度分析
利用本文提出的实验方法，得到4个实验段的局部海面高，如表 2所示。实验段A、B、C、D局部海面高中误差的均值分别为：0.01 m、0.03 m、0.01 m、0.01 m，计算的局部海面高具有较高的精度，可以满足厘米级海冰干舷计算的要求，说明本文提出的方法可行且精度较高。
| 表 2 实验段A、B、C、D海面高和干舷数据统计 Tab.2 Statistics ofSea Ice Freeboard and Sea Surface Height about Sections of A, B, C and D |
 |
图 6中圆点为含有冰间水道1 km分段利用本文提出的方法计算的局部海面高，星号点为没有冰间水道的分段根据拟合方程插值的局部海面高，虚线为沿轨拟合线。决定系数R2在统计学中，常用来评估模型的拟合优度。
 |
| 图 6 实验段A、B、C、D沿轨每公里局部海面高与沿轨距离散点图 Fig.6 Scatter Plots of Local Sea Surface Heights over 1 km Section with Leads Against the Along-Track Distance |
为实验段A、B、C、D局部海面高与沿轨距离散点图，沿轨线性拟合决定系数R2分别为：0.99、0.95、0.94、0.89，R2均接近或大于0.9，且残差均匀分布在拟合线两侧，说明线性模型和参数估计值是可靠的。因此，对于没有冰间水道的分段，在无法直接获取的情况下，通过线性拟合方程插值获得的海面高是合理且可靠的。
2.4 沿轨海冰干舷的验证与分析实验段A、B、C、D沿轨海冰干舷如图 7，计算方法如§1.4所述，先分别计算各实验段中每个1 km分段的海冰干舷，最终汇总得到每个实验段沿轨海冰干舷。
 |
| 图 7 实验段A、B、C、D沿轨海冰干舷分布 Fig.7 SeaIce Freeboard Distribution along-Track from Sections of A, B, C and D |
综合图 7和表 2，可以得到实验段A、B、C、D海冰干舷均值分别为：0.67 m、0.56 m、0.47 m、0.48 m，这4个实验段海冰干舷均值呈现出自西向东(由A到D)干舷较小的趋势，这与Kwok等[2]计算的威德尔海该航迹2011/2014/2016该航线海冰干舷的空间分布具有一致性。实验段A、B、C、D海冰干舷中误差分别为：0.34 m、0.29 m、0.26 m、0.28 m，与Wang等[7]利用DMS影像和R＜0.25在南极别林斯高晋海计算的2011年4个实验段海冰干舷中误差均值0.34 m、0.32 m、0.24 m、0.34 m结果精度相近。
从图 7和表 2干舷中误差参数可以看到，沿轨海冰干舷呈现不断起伏的状态，一方面这可能与计算结果不确定度有关；另一方面，是由于南北极海冰表面本身并不水平，海冰受洋流、大洋循环、潮汐、海风等的作用下，海冰之间破裂形成冰间水道和挤压形成的冰脊并存，ATM激光雷达沿轨条带扫描，扫描脚点仅为1 m，在雪面上脚点间隔3~4 m，其中条带边缘脚点间隔更为密集(亚m级)，因此ATM测量数据能够呈现出m级范围内海冰表面起伏状况，这在一定程度上影响干舷中误差大小。
3 结束语激光雷达是监测南北极海冰、冰盖和冰川变化的重要传感器。“冰桥计划”自2009年启动以来，在南北极地区已获得近10年机载激光雷达观测时间序列，延续了ICESat建立的激光雷达测高数据时间序列，同时，为研究区域尺度海冰干舷、厚度的空间分布和年际变化积累了大量高精度、高分辨率的测量数据。为了准确、程序化处理海量ATM测高数据，根据激光脚点在海冰和冰间水道内表面反射率和高程分布特征，提取冰间水道内部分高程较低的激光脚点用于提供海冰干舷计算所需的沿轨海面高程参考。
利用实验确定的Rcutoff能有效识别冰间水道中ATM L1B点，根据ATM L1B点在不同类型表面(开阔水面、薄冰、海冰)高程差异，由Rcutoff提取的脚点中取高程最低点计算局部海面高，其结果具有较高的精度，能够满足cm级海冰干舷计算要求。
通过结合ATM L1B测高数据表面反射率截止值和高程最低点法，能够达到程序化快速、准确计算海冰干舷的目标，后续将利用该方法处理大量OIB南极威德尔海2009~2017年海冰航线ATM L1B数据，计算海冰干舷，来分析干舷的厚度梯度和年际变化。
| [1] |
Zhou L, Xu S, Liu J, et al. On the Retrieval of Sea Ice Thickness and Snow Depth Using Concurrent Laser Altimetry and L-Band Remote Sensing Data[J]. The Cryosphere, 2018, 12: 993-1 012. DOI:10.5194/tc-12-993-2018 |
| [2] |
Lemke P. Observations: Changes in Snow, Ice and Frozen Ground[R]. Climate Change the Physical Science Basis, IPCC, New York, USA, 2007
|
| [3] |
Kwok R, Kacimi S. Three Years of Sea Ice Freeboard, Snow Depth, and Ice Thickness of the Weddell Sea from Operation IceBridge and CryoSat-2[J]. The Cryosphere, 2018, 12: 2 789-2 801. DOI:10.5194/tc-12-2789-2018 |
| [4] |
冯准准, 程晓, 康婧, 等. 美国NASA冰桥(IceBridge)科学计划:进展与展望[J]. 遥感学报, 2013, 17(2): 410-422. |
| [5] |
Onana V D P, Kurtz N T, Farrell S L, et al. A Sea-Ice Lead Detection Algorithm for Use with High-Resolution Airborne Visible Imagery[J]. IEEE Transactions on Geoscience and Remote Sensing, 2013, 51(1): 38-56. DOI:10.1109/TGRS.2012.2202666 |
| [6] |
Yi D, Harbeck J P, Manizade S S, et al. Arctic Sea Ice Freeboard Retrieval with Waveform Characteristics for NASA's Airborne Topographic Mapper (ATM) and Land, Vegetation, and Ice Sensor (LVIS)[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 53(3): 1 403-1 410. |
| [7] |
Wang Xianwei, Guan Fang, Liu Junyuan, et al. An Improved Approach of Total Freeboard Retrieval with IceBridge Airborne Topographic Mapper (ATM) Elevation and Digital Mapping System (DMS) Images[J]. Remote Sensing of Environment, 2016, 184: 582-594. DOI:10.1016/j.rse.2016.08.002 |
| [8] |
Forsberg R, Skourup H. Arctic Ocean Gravity, Geoid and Sea-Ice Freeboard Heights from ICESat and GRACE[J]. Geophysical Research Letters, 2005, 32(21): L21502. DOI:10.1029/2005GL023711 |
| [9] |
Wang X, Xie H, Ke Y, et al. A Method to Automatically Determine Sea Level for Referencing Snow Freeboards and Computing Sea Ice Thicknesses from NASA IceBridge Airborne LiDAR[J]. Remote Sensing of Environment, 2013, 131: 160-172. DOI:10.1016/j.rse.2012.12.022 |
| [10] |
Kwok R, Cunningham G F, Manizade S S, et al. Arctic Sea Ice Freeboard from IceBridge Acquisitions in 2009: Estimates and Comparisons with ICESat[J]. Journal of Geophysical Research, 2012, 117(C02018): 1-14. |
| [11] |
Kurtz N T, Farrell S L, Studinger M, et al. Sea Ice Thickness, Freeboard, and Snow Depth Products from Operation IceBridge Airborne Data[J]. The Cryosphere, 2013, 7(4): 1 035-1 056. DOI:10.5194/tc-7-1035-2013 |
| [12] |
Xie H, Ackley S F, Yi D, et al. Sea-Ice Thickness Distribution of the Bellingshausen Sea from Surface Measurements and ICESat Altimetry[J]. Deep Sea Research Part Ⅱ: Topical Studies in Oceanography, 2011, 58(9-10): 1 039-1 051. DOI:10.1016/j.dsr2.2010.10.038 |
| [13] |
Kwok R, Maksym T. Snow Depth of the Weddell and Bellingshausen Sea Ice Covers from IceBridge Surveys in 2010 and 2011: An Examination[J]. Journal of Geophysical Research. Oceans, 2014, 119(7): 4 141-4 167. DOI:10.1002/2014JC009943 |
| [14] |
Coleman R. Satellite Altimetry and Earth Sciences: A Handbook of Techniques and Applications[J]. Eos Transactions American Geophysical Union, 2001, 82(34): 376. |
| [15] |
Zwally J H, Donghui Yi, Kwok R, et al. ICESat Measurements of Sea Ice Freeboard and Estimates of Sea Ice Thickness in the Weddell Sea[J]. Journal of Geophysical Research, 2008, 113(C12013): 1-19. |
| [16] |
Studinger M. IceBridge ATM L1B Elevation and Return Strength, Version 2[C]. NASA National Snow and Ice Data Center Distributed Active Archive Center, Boulder, Colorado, USA, 2018 http://nsidc.org/data/docs/daac/icebridge/ilatm1b/
|
| [17] |
Krabill W B, Abdalati W, Frederick E B, et al. Aircraft Laser Altimetry Measurement of Elevation Changes of the Greenland Ice Sheet: Technique and Accuracy Assessment[J]. Journal of Geodynamics, 2002, 34(3-4): 357-376. DOI:10.1016/S0264-3707(02)00040-6 |
| [18] |
Yi D, Kurtz N, Harbeck J, et al. Comparing Coincident Elevation and Freeboard From IceBridge and Five Different CryoSat-2 Retrackers[J]. IEEE Transactions on Geoscience and Remote Sensing, 2019, 57(2): 1 219-1 229. DOI:10.1109/TGRS.2018.2865257 |
| [19] |
Wang X, Holland D M. A Method to Calculate Elevation-Change Rate of Jakobshavn Isbrae Using Operation IceBridge Airborne Topographic Mapper Data[J]. IEEE Geoscience and Remote Sensing Letters, 2018, 15(7): 981-985. DOI:10.1109/LGRS.2018.2828417 |
| [20] |
Martin C F, Krabill W B, Manizade S S, et al. Airborne Topographic Mapper Calibration Procedures and Accuracy Assessment[R]. Maryland, U.S: National Aeronautics and Space Administration, 2012
|
| [21] |
Dominguez R. IceBridge DMS L1B Geolocated and Orthorectified Images[C]. NASA National Snow and Ice Data Center Distributed Active Archive Center, Boulder, Colorado, USA, 2018
|