2. 海岛(礁)测绘技术国家测绘局重点实验室, 青岛 266510;
3. 福州大学环境与资源学院, 福州 350108;
4. 中国测绘科学研究院, 北京 100039;
5. 台湾交通大学土木工程学系, 台湾新竹 300
2. Key Laboratory of Surveying and Mapping Technology on Island and Reef of SBSM, Qingdao 266510, China;
3. College of Environment and Resources, Fuzhou University, Fuzhou 350108, China;
4. Chinese Academy of Surveying and Mapping, Beijing 100039, China;
5. Department of Civil Engineering, National Chiao Tung University, Hsinchu 300, China
卫星测高技术能够大范围、高精度、周期性地探测海洋上的多种自然现象及其变化, 其观测数据被广泛地应用于地球重力场模型、平均海面、海潮模型、海底地形和洋流等方面的研究[1, 2].在现代科技条件下, 卫星测高仪在深海区域测得海面高的精度可以达到±2~3cm.然而, 在近岸海域, 沿海陆地地形、岛屿、潮汐、地球物理因素和仪器硬件响应等的影响造成雷达测高脉冲的反射波形不规则, 无法精确地求得卫星与星下点反射面之间的距离[3].因此, 为了得到精确的距离观测值, 需要重新确定波形前缘中点的位置, 并通过与波形预设门比较, 计算出波形重定距离改正值, 对地球物理数据记录(Geophysical Data Record, GDR)中给出的卫星到海面的距离进行改正, 此过程称为波形重定.对卫星测高波形进行波形重定的主要目的是通过对波形前缘中点的重新确定来减少短波随机噪声, 从而改善GDR中瞬时海面高的精度[3].在近岸海域, 由于对流层延迟、电离层延迟、地球物理改正、轨道和仪器校正等改正技术较为成熟, 因此利用波形重定方法来改善卫星测高数据质量成为十分有效的手段.
EnviSat是欧洲空间局(European Space Agency, ESA)于2002年3月发射的太阳同步卫星, 轨道高度为800km, 轨道倾角为98°, 重复周期为35天(http://envisat.esa.int/).发布的EnviSat测高数据中SGDR (Sensor Geophysical Data Record)为传感器数据, 它除了包含GDR数据外, 还包含由离散傅里叶变换和快速傅里叶变换算法分别处理得到的18 Hz的波形数据.GDR数据中不仅包含常规数据, 还包含了EnviSat利用Ice-1、Ice-2、Sea-ice和Ocean1等重定算法改正后的海面高.Ice-1和Ice-2算法主要是针对冰面进行波形重定; Sea-ice算法是针对海冰面回波进行波形重定; Ocean1算法是针对海洋面回波进行波形重定[4].针对近岸海域, 本文研究中只采用了EnviSat利用星载的Ocean1算法重定的海面高作为比较研究对象.
为了改善近岸海域的卫星测高数据质量, 本文在借鉴已有波形重定经验的基础上, 对地中海区域近海岸的EnviSat18Hz的测高波形数据进行了波形重定方法研究, 设计了卫星测高波形重定软件, 对EnviSat进行波形重定, 并与ESA所提供的EnviSat利用星载Ocean1算法重定后的海面高以及EGM2008大地水准面起伏进行了比较.
2 波形重定方法波形重定方法主要用来确定出波形实际的前缘中点与设计门之间的偏移量, 然后根据公式计算出距离改正值.目前主要的波形重定方法有拟合算法(如β参数法)、重心偏移法(OCOG)、阈值(Threshold)方法等[3, 5].本文研究中除了利用上述方法外, 还采用了Guo等[5]提出的改进Threshold算法, 并且在研究中对该方法进行了进一步的优化.
2.1 β参数算法β参数算法[6]是美国国家航空航天局(NASA)的Martin等[7]于1983提出的, 该方法基于Brown平均脉冲反射模型, 采用适当的参数函数对测高仪波形进行拟合.β参数算法需要有良好的初始值, 利用最小二乘平差或极大似然估计方法, 通过迭代计算求得β参数, 其中5-β参数法主要用于求解单一反射面反射的复杂波形, 如图 1所示.但是如果波形呈现尖锥形状时, 使用β参数算法将会使计算迭代不收敛而无法得到计算结果.本文研究中采用了5-β参数法的线性形式[3, 6, 8~12], 即
(1) |
式中,
y(t)为t时的采样功率; β 1为返回波形的热噪声; β 2为回波功率强度; β 3为前缘中点, 由接收到能量至最大振幅的一半, 可与预设门差值计算距离改正量; β 4为波形前缘斜率, 并提供了星下点有效波高信息; β 5为波形后缘斜率; P(x)为误差函数; Q为线性函数, 它模拟波形后缘逐渐衰减的回波.
利用最小二乘法求解非线性方程(1)确定5个参数的过程中, 未知参数的初始值和波形采样功率加权的确定对计算结果影响显著.在本文研究中, 未知参数的初始值和权的确定采用了Anzenhofer等[9]在处理ERS-1波形资料时所使用的方法.
2.2 OCOG (Off-Center of Gravity)算法Wingham等[13]于1986年提出了OCOG算法, 如图 2所示.OCOG算法是基于统计规律得到的简单波形重定算法, 其基本思想是找到每个返回波形的重心, 以数值方式统计出波形振幅、宽度与重心位置.该方法计算简单方便, 但与反射表面的物理意义无关.为了降低热噪声对波形的影响, Deng [3]在2004年对OCOG算法进行了改进.OCOG算法的公式为
(2) |
(3) |
(4) |
(5) |
式中, N为阀门的总个数; n为波形开始和结束时应剔除的偏差波形个数; Pi(t)为波形的第i个阀门的功率值; A为振幅; W为宽度; COG (Center of Gravity)为波形的重心; LEP (leading edge position)为前缘中点.
2.3 Threshold算法Threshold算法是Davis在1997年提出的, 利用OCOG计算的矩形, 根据振幅、最大波形采样等给出门槛值, 在与前缘陡峭部分相交门槛的几个临近采样值之间进行线性内插, 以确定重定点[14].它是一种统计算法, 具有OCOG算法的优点, 并且比OCOG算法所确定的跟踪门的位置更精确. Threshold算法[11, 14]的公式为
(6) |
(7) |
(8) |
式中, Pi、Pk分别为第i个、第k个阀门的功率值; PN为波形前5个阀门功率值的平均值; Th为门槛值; Gk为第k个大于Tl的阀门; Gr为前缘中点.
2.4 优化的Threshold算法在顾及波形物理机制基础上, Guo等[5]提出了改进的Threshold波形重定算法, 其主要思想:首先进行波形分析, 找出波形中所包含的子波形; 然后利用改进的波形公式进行重定.所谓子波形, 是指由于测高仪发射的微波脉冲受到陆地、海岸等地形影响, 在卫星观测波形中形成的具有多个波形前缘的波形, 这些多个波形前缘称为子波形[15].其判断为设第i个阀门对应的回波功率为Pi, 令回波间隔功率差的均值为
(9) |
整个波形回波间隔功率差的标准差为
(10) |
相邻回波功率差为
(11) |
整个波形回波单差的标准差为
(12) |
图 3为子波形的判断流程图.在本文研究中, 判断子波形时采用的开始阀门i为10, 结束阀门n为118, 子波形的判断标准为ε1=0.2S, ε2=0.2S1, 门槛值为0.3.
利用上述方法, 对整个波形进行搜索, 确定一个波形中所含的所有子波形, 然后利用公式(13)确定出每个子波形的功率门槛值:
(13) |
式中, A由公式(2)确定; Tn为子波形中第二个阀门的波形功率值.
对于每个子波形, 其波形前缘中点位置可以通过公式(2)、(6)、(8)和(13)计算得到.利用此方法, 一个波形可能确定多个前缘, 计算得到多个重定后的距离改正, 但在这些改正中只有一个是准确的, 本文研究中采用了最小重定距离改正的方法确定重定后的距离.
2.5 波形重定距离改正计算经波形重定后, 可以确定出波形的实际前缘中点的位置; 根据波形设计阀门的位置以及光速, 可以计算出重定距离改正dr.
(14) |
式中, ΔGa为一个阀门时间间隔(EnviSat中当脉冲宽度为320 MHz时, 阀门时间间隔为3.125ns); c为光速(即299792458 m/s); Gr为波形重定后的前缘中点的位置; G0为波形设计阀门的位置.对于不同的测高卫星, 其波形设计阀门的位置不同, 而对于同一颗测高卫星当采用不同的数据记录格式或者同一记录格式下采用的阀门起止范围不同时, 其设计阀门的位置也不相同.如EnviSat的设计阀门位置, 当阀门起止范围是0~127时, 其设计阀门为45;当阀门起止范围是1~128时, 其设计阀门为46[4].
3 波形重定软件设计基于上述波形重定方法, 本文利用FORTRAN和Matlab语言进行了波形重定软件设计.如图 4所示, 该软件通过对测高波形进行重定, 以及引入外部判断准则, 最后得到改正后的海面高(Sea Surface Height, SSH)、海平面异常值(Sea Level Anomaly, SLA)以及IMP值(Improvement Percentage, 见公式(15)).重定软件中各组成部分的功能如下:
(i) SGDR数据的读取与编辑:使用默认的标志读取并编辑SGDR数据, 输出18 Hz的波形数据并绘制波形图;
(ii)波形重定:使用所有的重定方法对波形进行重定.首先需要利用OCOG算法对波形进行重定, 然后将其结果作为其他重定算法的先验值; 对于5-β算法, 在重定时要引入波形采样功率的加权值及迭代溢出准则; 而改进的Threshold算法需进行子波形判断, 然后利用OCOG和Threshold算法计算出各波形的波形前缘中点, 并通过判断选取最优的前缘中点.
(iii)重定结果评判:引入EGM2008模型的大地水准面高数据作为重定结果的评判准则, 然后输出重定后的SSH、SLA以及IMP值.
4 EnviSat测高波形重定地中海是世界上最大的陆间海, 面积约250多万km2, 平均深度1450 m, 最深处5092 m, 海域中有南欧三大半岛及西西里岛、撒丁岛、科西嘉岛等岛屿[16].本研究选择地中海海域为实验区, 给出EnviSat几种波形重定方法的结果.
EnviSat有多个Pass通过地中海, 在本文研究中选取了其中Cycle59的Pass801作为研究对象(如图 5所示), 该Pass穿过科西嘉岛, 经过撒丁岛东部海域, 最近距离撒丁岛15.1km左右.由于40°N~43.3°N之间的GDR数据质量太差而无法使用, 所以研究区域选为38°N~40°N和43.3°N~44.3°N.本研究采用了5种波形重定方法对EnviSat的18Hz (128bins)测高波形进行了重定, 并将重定结果与EGM2008[17]确定的大地水准面起伏进行比较, 通过计算IMP值判断重定结果的质量.
利用IMP值来判断重定后海面高质量的方法是Hwang等[8]于2006年提出的, 其计算公式为
(15) |
式中, δraw为未重定的海面高与大地水准面高差值的标准差; δretracked为重定后的海面高与大地水准面高差值的标准差.若IMP为负值, 表明加入重定距离改正后的海面高比未重定前的海面高精度低.
表 1和表 2给出了在两个研究区域内各波形重定方法的波形重定成功率和IMP值, 其中EnviSat-ocean1的各项数据是由EnviSat的GDR数据中所提供的利用星载Ocean1重定算法[4]重定的海面高(18Hz)计算得到.除5-β参数法在两个区域的波形重定成功率分别为86.97%和78.27%外, 其他方法的波形重定成功率均为100%.利用不同重定方法计算得到的IMP值不同, 其中OCOG算法的重定结果最差, 主要是由于OCOG算法使用了全部的波形数据, 而当波形受到近海岸地形、地球物理环境、硬件和偏离星下点等的影响时, 波形中含有较大的噪声, 使得利用OCOG算法得到的海面高精度较低, 所以OCOG算法主要用来为其他算法求定初始值.利用本文优化的Threshold算法和Guo等的原始算法计算得到的IMP值较其他算法计算得到的IMP大, 表明该算法的重定效果最好.
在本文研究中, 采用了以下两个准则判断波形重定方法的好坏: (1)波形重定成功率, 即可重定波形个数与波形总个数之比; (2) IMP值和重定后海面高的平滑度, 所谓海面高的平滑度是指相邻点的海面高之间的差值要小, 不能产生大的突变, 保持良好的光顺性.图 6和图 7是利用本文提到的5种重定算法计算得到的海面高与由EGM2008模型计算得到的大地水准面高以及EnviSat的SGDR数据中提供的利用星载Ocean1算法重定的海面高的对比图.比较图 6和图 7可知, 利用Threshold算法和改进的Threshold重定后的海面高较利用星载Ocean1算法重定的海面高平滑, 并且重定后的海面高与EGM2008模型计算得到的大地水准面高更相似, 且两种改进的threshold方法重定后的海面高平滑等十分接近.综合表 1和表 2可知, 近海岸波形重定中优化的Threshold算法是目前较好的一种重定算法.对比表 1和表 2可知, 两个试验区域利用相同的波形重定方法进行波形重定后, 海面高的改善程度不同.在38°N~40°N范围内, Pass801的地面轨迹距岛屿和海岸较远, 最近处为40.8km, 除38.0°N~38.6°N的水深达到-390m~-1300m外, 其余区域的水深都低于-1700m, 测高波形为海洋波形和似海洋波形, 其质量较好, 因此GDR中未重定的海面高的精度较高, 所以重定后海面高的改善效果不明显; 在43.3°N~44.3°N范围内, Pass801的地面轨迹距岛屿和海岸较近, 最近处为6.6km左右, 水深最浅处为-208m左右, 大部分区域的水深高于-1700m, 离最近海岸的距离均小于40km, 测高波形为似海洋波形和阶梯状波形, 其质量相对较差, 因此GDR中未重定的海面高的精度相对较低, 所以重定后海面高的改善效果明显.因此, 在近岸区域进行波形重定是十分必要的, 可以显著提高卫星测高数据的精度及其可用数据量.
Vignudelli et al. [18]曾利用TOPEX/Poseidon (T/ P)的GDR数据获得科西嘉海峡的SLA, 并与利用MSS CLS01和LEGOS计算得到的平均海平面进行比较, 得到较好的结果.Fenoglio-Marc等[19]利用T/P的GDR和R-GDR数据以及ERS-2的GDR数据对地中海区域的SLA进行了研究, 并将结果与潮汐测量数据进行了比对.在他们的研究中仅使用了测高卫星的GDR数据没有进行波形重定的研究.与他们的研究相比, 采用波形重定后的GDR数据得到的SSHs的精度显著提高, SLA偏差在±15cm之间, 并且无需在一定限制条件对GDR数据进行预处理.
5 结论本文探讨了近岸海域卫星测高波形重定方法, 以地中海为例利用5种波形重定方法对EnviSat测高波形进行了重定研究.由于OCOG算法使用了全部的波形数据, 而当波形受到近海岸地形、地球物理因素、硬件和偏离星下点等的影响时, 波形中含有较大的噪声, 使得利用OCOG算法得到的海面高精度较低, 所以OCOG算法主要用来为其他算法求定初始值.除OCOG算法外, 其他波形重定方法可以显著提高EnviSat在近岸海域GDR数据的精度, 其中优化的Threshold方法的重定结果是最好的.这说明该方法较适合于近岸测高波形的重定, 但该方法对于不同的测高卫星需采用不同的子波形的判断标准.当IMP值越大, 重定后的海面高越平滑, 那么这种波形重定方法的重定效果就越好.通过研究表明, 波形重定方法可以有效改善近岸海域卫星测高GDR数据的质量, 增加近海岸海域测高数据的利用率, 使得测高数据可以用近岸区域大地测量、大地水准面和重力异常确定、地球物理、物理海洋等方面的研究.
致谢感谢德国Technische University Darmstadt提供的EnviSat的SGDR数据.
[1] | Fu L L, Cazenave A. Satellite altimetry and earth sciences:a handbook of techniques and applications. San Diego: Academic Press, 2001 . |
[2] | 王广运, 王海英, 许国昌. 卫星测高原理. 北京: 科学出版社, 1995 . Wang G Y, Wang H Y, Xu G C. The Theory of Satellite Altimetry (in Chinese). Beijing: Science Press, 1995 . |
[3] | Deng X L.Improvement of geodetic parameter estimation in coastal regions from satellite radar altimetry[Ph.D.thesis].Perth:Curtin University of Technology, 2004 http://www.oalib.com/references/15818946 |
[4] | European Space Agency (ESA).EnviSat RA2/MWR Product Handbook.ESA Issue 2.1, 2006 |
[5] | Guo J Y, Hwang CW, Chang X T, et al. Improved threshold retracker for satellite altimeter waveform retracking over coastal sea. Progress in Natural Science , 2006, 16(7): 732-738. DOI:10.1080/10020070612330061 |
[6] | Martin T V, Zwally H J, Brenner A C, et al. Analysis and tracking of continental ice sheet radar altimeter waveforms. J Geophys.Res. , 1983, 88(C3): 1608-1616. DOI:10.1029/JC088iC03p01608 |
[7] | Brown G S. The average impulse response of a rough surface and its applications. IEE-E Trans Ant Pop , 1977, AP-25(1): 67-74. |
[8] | Hwang C, Guo J, Deng X, et al. Coastal gravity anomalies from retracked Geosat/GM altimetry:improvement, limitation and the role of airborne gravity data. J Geod , 2006, 80: 204-216. DOI:10.1007/s00190-006-0052-x |
[9] | Anzehofer M, Shum C K, Rentsh M.Coastal altimetry and applications.Rep 464, Dept Geod Sci and Surveying, Ohio State University, Columbus, 1999 |
[10] | Hwang C, Hsu H, Deng X. Marine gravity anomaly from satellite altimetry:a comparison of methods over shallow waters. Berlin:Springer, IrUernational Association of Geodesy SymPosia , 2003, 126: 59-66. DOI:10.1007/978-3-642-18861-9 |
[11] | 高永刚, 郭金运, 岳建平. 卫星测高在陆地湖泊水位变化监测中的应用. 测绘科学 , 2008, 33(6): 73–75. Gao Y G, Guo J Y, Yue J P. Lake level variations measured with satellite altimetry (in Chinese). Science of Surveying and Mapping (in Chinese) , 2008, 33(6): 73-75. |
[12] | 褚永海, 李建成, 张燕, 等. EnviSat测高数据波形重跟踪分析研究. 大地测量与地球动力学 , 2005, 25(1): 76–80. Chu Y H, Li J C, Zhang Y, et al. Analysis and investigation of waveform retracking data of EnviSat. Journal of Geodesy and Geodynamics (in Chinese) , 2005, 25(1): 76-80. |
[13] | Wingham D J, Rapley C G, Griffiths H. New techniques in satellite altimeter tracking systems. Proc IGARSS'86 Symp, Zurich , 1986: 1339-1344. |
[14] | Davis C H. A robust threshold retracking algorithm for measuring ice-sheet surface elevation change from satellite radar altimeter. IEEE Trans.Geosci.Remote.Sens. , 1997, 35(4): 974-979. DOI:10.1109/36.602540 |
[15] | Guo J Y, Gao Y G, Hwang C W, et al. A multisubwaveform parametric retracker of the radar satellite altimetric waveform and recovery of gravity anomalies over coastal oceans. Sci.China Ser D-Earth Sci. , 2009. DOI:10.1007/s11430-009-0171-3 |
[16] | 地中海.http://baike.baidu.com/view/15817.htm, [2008.12.8] Mediterranean Sea.http://baike.baidu.com/view/15817.htm, [2008.12.8] |
[17] | Pavlis N K, Holmes S A, Kenyon S C, et al.An earth gravitational model to degree 2160:EGM2008.The 2008General Assemly of the European Geoscience Union, Vienna.Apr 13-18, 2008 |
[18] | Vignudelli S, Cipollini P, Roblou L, et al. Improved satellite altimetry in coastal systems:Case study of the Corsica Channel (Mediterranean Sea). Geophys.Res.Lett. , 2005, 32: L07608. |
[19] | Fenoglio-Mare L, Vignudelli S, Humbert A, et al.An assessment of satellite altimetry in proximity of the Mediterranean coastline 3rd EnviSat Symposium Proceedings, ESA Publications Division, 2007.SP-636 |