专利名称：基于星载激光测高仪回波的海洋表面风、浪特性反演方法
技术领域：
本发明涉及一种海洋表面风、浪特性反演方法，尤其是涉及一种基于星载激光测高仪回波的海洋表面风、浪特性反演方法。
背景技术：
海洋表面波浪斜率是重要的海态参数，同时波浪斜率还与海面上空风场强相关。目前对海洋表面波浪斜率的遥感测量主要由机载或星载微波雷达进行根据雷达反射回波计算目标散射截面，并由此反演海面平均斜率(MSS, mean square slope),进而通过风浪模型推算海平面上方风速。星载激光测高仪是一种安置于低轨卫星平台，通过发射激光脉冲、接收激光脉冲经地物反射后的微弱回波，由两者之间的时间差以及卫星平台的姿态、位置计算地物高程的新型卫星载荷。对地观测激光测高仪如果具有全波形记录能力，则可以进 一步演变成同时具备测距和反演目标表面特性的激光雷达遥感系统。从上世纪70年代开始用于海洋表面测量的星载微波雷达高度计兴起并逐渐成熟，成为全球范围内海洋表面波浪形态观测的主要手段。相对于微波雷达测高仪，激光测高仪具有更小的光束发散角和工作波长、更精确的空间定位能力和更高的测量精度，并具备使用波形参数来反演被测目标信息的功能，是对地高精度测量与遥感的重要发展方向。世界上首台具有全波形记录能力的对地观测星载激光测高仪是美国2003年发射的ICESat卫星上搭载的地球科学激光测高系统一GLAS系统。该系统在轨运行的3年中，记录下了包括陆地、海洋、海冰和极区冰盖等不同地物的海量回波波形。由于GLAS系统主要用作对极区(主要是南北极)冰盖的测量，从回波形态中分析海洋表面特性、反演海洋表面风速的工作成果少见报道。
发明内容
本发明主要是解决现有技术所存在的技术问题；提供了一种基于星载激光测高仪回波的海洋表面风、浪特性反演方法。本发明的上述技术问题主要是通过下述技术方案得以解决的一种基于星载激光测高仪回波的海洋表面风、浪特性反演方法，其特征在于，包括以下步骤步骤1，设定采样点数以及窗口长度，利用GLAS测高仪系统获取海洋回波波形信号后，计算信号标准差Ss;并计算回波白噪声标准差δη;根据得到的63和6 获取高斯滤波器宽度σ ;步骤2，由所得高斯滤波器宽度σ通过连续型高斯滤波器基于高斯滤波函数进行采样并归一化处理，得出离散高斯滤波器参数hk ；步骤3，根据步骤3中得到的数字滤波器参数hk获取滤波后信号；步骤4，对步骤3中滤波后信号做参数提取，得到单高斯回波的峰值Ai，脉宽σ 1和重心点Ti位置；步骤5，根据步骤4中得到的脉宽σ i基于测高仪系统参数测高仪轨道高度Z，光束发散角θτ，发射脉冲宽度Of，接收系统展宽Oh和真空光速C，获取海平面上方风速W。在上述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，所述的步骤I中，回波白噪声标准差δη的具体获取方法如下a.计算接收离散波形Sk所有点的纵坐标幅度均值Smean ；b.计算从波形尾端小于Snrean的第一个 点起始到波形末尾点集合的标准差作为δ ；3.根据权利要求2所述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，其特征在于，所述的步骤3中，滤波后信号的获取是根据对数字滤波器参数hk与回波离散信号Sk作卷积运算得到。在上述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，所述的步骤4中，参数提取的具体方法如下步骤4. I.对步骤3中滤波后信号求二阶差分，找出全部拐点；步骤4. 2.针对步骤4. I中的全部拐点获取满足条件的η个拐点对满足条件波形中第2i-l个是正向直线外点，且第2i个是正向直线内点的拐点对，其中i=l，2，…η。步骤4. 3.在第i组拐点对区间内中找出最大值作为预估计高斯振幅Aum，两拐点均值为预估计高斯重心位置Tum,找出区间内信号值等于O. 8倍预估高斯振幅两侧点位置，其两点间距为预估计高斯脉宽oipM;步骤4. 4.以估计的高斯参数为初始波形，基于梯度下降法或梯度下降-高斯牛顿法的LM算法对滤波后信号进行高斯拟合，当达到设定拟合收敛条件时，拟合结束得出各个高斯分量的拟合参数AyTi和0i。在上述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，所述的步骤I中，高斯滤波器宽度σ基于公式(式15):O = Ik- — }
\ δ、/其中k为常正整数，O为最近取整函数，δ η为权利要求2中计算的回波白噪声标准差，δ s为信号标准差，可以通过计算回波信号当前点与窗口长度为j的前后2j+l个回波信号点的标准差获得，通常情况下长度2j+l约为5%左右信号总采样点长度。在上述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，所述高斯滤波函数为
I — 12Ii(I) ='..1....广.......exp(——Γ)
ν2πσ 2σ-其中，σ为滤波器宽度，exp O为指数函数，t为横坐标时间轴。在上述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，所述的步骤4.4中，定义收敛条件为ε ( δη，基于梯度下降法或梯度下降-高斯牛顿法的LM算法对
滤波后信号进行高斯拟合，当拟合满足条件JyL (sLZwI)Ig5时，拟合结束，其中Sk为
y I L
离散激光高度计测高仪回波信号，L为回波信号的采样点数，Wk为式(12)的离散采样点，其采样间隔与回波信号相同。其中，LM算法是通用算法，输入初始参数，预拟合函数类型和拟合结束条件后，其会运算自收敛或无法收敛；这里输入初始参数是4. 1-4. 3步骤中计算的Ai,, Ti,和oipM;拟合收敛条件为ε ( δη;预拟合函数类型为12式，拟合结果是得出12式中的参数Ai,凡和σ it)在上述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，还包括一个步骤6 :即由回波波形反演得到风速数据之后，利用海面均方斜率公式S2=O. 003+0. 00512w能够获取海面平均斜率；根据风浪波高公式SWH=4 σ ξ=0. 064w2能够获取波浪高度。因此，本发明具有如下优点因此，本发明主要有以下优点1)激光测高仪光束发散角通常为亚毫弧量级，远小于雷达测高仪微波波束发散角的特点，具有更高的脚点定位精度和空间分辨率。2)由于使用回波的脉冲宽度反演风速，与回波整体能量无关，不受大气衰减变化的影响，具有较强的通用性和抗干扰性。3)与目前通用的雷达测高仪使用雷达散射截面反演风速的原理和途径完全不同，可以与雷达测高仪的测量数据进行校验或补充。
图I是本发明中天顶方向入射的激光测高仪探测海洋表面回波的工作示意图。图2是本发明应用GLAS系统发射的激光脉冲时间波形,其中,横轴为相对时间轴，单位为ns ;纵轴为幅值，单位为V，援引GLAOl i_tx_wf V。图3是本发明中GLAS海洋回波处理结果，其中，实线表示海洋回波的原始波形GLAOl i_rng_wf V，虚线为经过自适应高斯滤波和高斯拟合后波形，具体为图6中第一组海洋回波数据。图4是本发明GLAS公布的海洋回波高斯滤波处理后的波形结果，与图3处理结果对比有很好的一致性，两者拟合具体参数在图6中给出。图5是2003年2月21日零时NCEP公布的全球表面上方12m位置的风速信息，NCEP提供全球数值天气分析数据集，数据集在格林尼治时间00:00、06:00、12:00和18:00时更新，其中，包含地表大气压、风速、可降水量和以标准大气压层给出的温度、相对湿度和位势高度数据。其全球坐标纬度为北向为正，南向为负；经度为从O度经线向东为正到359度，360度经线与O度重合；其中风速数据援引NCEP第280行，281行的经向风和纬向风的均方根值。图6是使用GLAS海洋回波计算波形参数结果与GLAS官方处理公布结果对比。高斯拟合后的参数包括峰值，重心位置和脉冲宽度，其中预处理后提取出的脉冲宽度用于海平面上方风速反演计算。表中UTC时间为以2000年I月I日12时为基准O点的相对时间，实际为波形编号；重心位置为高斯重心与最后一个回波点的时间差值，因此为负数。使用该方法处理回波波形误差均小于1%。表I中所有海洋回波均为2003年2月21日00:13分到00:40分数据，与图5中NCEP气象数据时间接近，便于减小误差；GLAS官方数据是援引 GLA05 d_para2
中分别给出 Std Fit Peak IAmplitude Sigma (V), Std FitPeak I Location Sigma (ns), Std fit Peakl Width Sigma (ns),即波形数据处理结束后计算出的回波信号高斯波形的振幅，重心位置和均方根脉宽。图7是反演海洋表面上方风速数据与NCEP公布风速对比。表I中海洋回波通过计算处理后的脉冲宽度使用10式反演海平面上方风速W，其中GLAS系统参数为轨道高度z=600km,光束发散角Θ T=110urad,采样间隔Δ t=lns,对应系统展宽σ h=l/12°_5ns,发射脉宽σ f使用同样方法做滤波拟合提取高斯参数；通过对比结果反演风速准确度很好。
图8是海洋表面斜率和波高计算数据。通过表2反演风速结果，根据Cox-Munk海面均方斜率公式S2=O. 003+0. 00512w计算海面平均斜率；根据Phillips风浪波高公式SWH=4 ο ξ = O. 064w2 计算波浪高度(Significant wave height),其中 w 为风速，σ ξ 为 RMS波闻。
具体实施例方式下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明。实施例首先介绍一下本发明所需要的理论基础 I.回波脉宽与海面风速、波浪斜率之间关系的理论推导根据Gardner的理论，发射激光脉冲经过一次菲涅尔衍射，入射到地球表面经过目标反射，再经过一次菲涅尔衍射，到达接收望远镜视场，望远镜接收到回波功率由(I)式表示，式中P为反射目标横截面坐标矢量，β (P)光斑内目标反射率，ξ (P)目标表面高程轮廓，Ak为接收望远镜面积，ζ为测高仪轨道高度，c为真空光速，Σ目标光斑照射区域，a(P，ζ)和f(t)分别为激光脉冲在空域和时域分布函数。P{t)=^ JJβ{ρ)■ \a(p, ■ [/Ir-^-￡1+d2p(I)
Z~I c cz c通常情况下，激光发射脉冲在空域和时域分布都可近似为高斯函数，分别由(2)式和(3)式表示。式中Q为发射激光脉冲的能量，θτ为初始光束发散角，Of为发射脉冲RMS宽度
,VIa(p,r)l' =—-^expf^p2 Hz1 tan2 θτ)(2)
2π(ζ &ηθ Y
丨爾丨2 =-/=—6ΧΡ("7~τ)(3)^2χσ 2σ按照Kodis理论，在近天顶入射的条件下，海洋表面反射率由表面随机分布镜面点的后向散射决定，散射功率正比于被照亮的镜面点个数η。每个镜面点散射可以看作是一个正切球，球半径等于镜面点表面的长短轴半径1^和rb的几何平均。在光学领域散射截面O=Ji δ〈 |rarb|〉，归一化到 Ω =4 π 的球面度,并代入 Barrick 的结论 η〈 | rarb |〉= (1+P 2/ζ2)2ρ[ξχ, ξγ ξ (P)]后，海洋表面反射率β由(4)式表示β{ρ) = Η\γ \) = δ(1 + ￡1)-ρ[ξ% j \ξ{ρ)]/4 (4)
4(4)式中δ是海水与大气之间镜面反射率，与入射激光波长λ有关；P [ ξ χ，ξ y I ξ (ρ)]为海面斜率在波高为ξ时的条件概率密度，ξ χ和I,分别表示光斑照明区域海水表面横截面χ和y方向波浪斜率。根据Cox-Munk经验公式，在忽略波浪毛细波和地球自转等影响的理想情况下，波高Ρ( ξ )和波浪表面轮廓ρ( ξ χ, ξ y)满足高斯分布，且仅与风速有关。由于波浪同一点处的高度和斜率是不相关的，则波高和波浪斜率的联合概率密度函数可由(5)式表示
权利要求
1.一种基于星载激光测高仪回波的海洋表面风、浪特性反演方法，其特征在于，包括以下步骤 步骤1，设定采样点数以及窗口长度，利用GLAS测高仪系统获取海洋回波波形信号后，计算信号标准差Ss;并计算回波白噪声标准差δη;根据得到的63和6 获取高斯滤波器宽度σ ; 步骤2，由所得高斯滤波器宽度σ通过连续型高斯滤波器基于高斯滤波函数进行采样并归一化处理，得出离散高斯滤波器参数hk ； 步骤3，根据步骤3中得到的数字滤波器参数hk获取滤波后信号； 步骤4，对步骤3中滤波后信号做参数提取，得到单高斯回波的峰值Ai，脉宽σ i和重心点Ti位置； 步骤5，根据步骤4中得到的脉宽σ i基于测高仪系统参数测高仪轨道高度z，光束发散角θτ，发射脉冲宽度Of，接收系统展宽和真空光速C，获取海平面上方风速W。
2.根据权利要求I所述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，其特征在于，所述的步骤I中，回波白噪声标准差Sn的具体获取方法如下 a.计算接收离散波形Sk所有点的纵坐标幅度均值S_n； b.计算从波形尾端小于S_n的第一个点起始到波形末尾点集合的标准差作为δη;
3.根据权利要求2所述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，其特征在于，所述的步骤3中，滤波后信号的获取是根据对数字滤波器参数hk与回波离散信号Sk作卷积运算得到。
4.根据权利要求I所述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，其特征在于，所述的步骤4中，参数提取的具体方法如下 步骤4. I.对步骤3中滤波后信号求二阶差分，找出全部拐点； 步骤4. 2.针对步骤4. I中的全部拐点获取满足条件的η个拐点对满足条件波形中第2i-l个是正向直线外点，且第2i个是正向直线内点的拐点对，其中i=l，2，…η。
步骤4. 3.在第i组拐点对区间内中找出最大值作为预估计高斯振幅Aum，两拐点均值为预估计高斯重心位置Tum,找出区间内信号值等于O. 8倍预估高斯振幅两侧点位置，其两点间距为预估计高斯脉宽Oi,; 步骤4. 4.以估计的高斯参数为初始波形，基于梯度下降法或梯度下降-高斯牛顿法的LM算法对滤波后信号进行高斯拟合，当达到设定拟合收敛条件时，拟合结束得出各个高斯分量的拟合参数AyTi和0i。
5.根据权利要求I所述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，其特征在于，所述的步骤I中，高斯滤波器宽度σ基于公式
6.根据权利要求I所述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，其特征在于，所述高斯滤波函数为
7.根据权利要求4所述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，其特征在于，所述的步骤4.4中，定义收敛条件为ε ( δη，基于梯度下降法或梯度下降-高斯牛顿法的LM算法对滤波后信号进行高斯拟合，当拟合满足条件
8.根据权利要求I所述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，其特征在于，还包括一个步骤6 :即由回波波形反演得到风速数据之后，利用海面均方斜率公式S2=O. 003+0. 00512w能够获取海面平均斜率；根据风浪波高公式SWH=4 σ ξ=0. 064w2能够获取波浪高度。
全文摘要
本发明涉及基于星载激光测高仪回波的海洋表面风、浪特性反演方法。本发明首先对星载激光测高仪回波波形进行预处理和特征提取，再建立提取的回波参数与海洋表面的风、浪场之间的数学关系，即可根据回波参数来反演计算海洋表面波浪斜率、波高和海洋表面风速。由于激光光束具有远小于微波波束的发散角的特点，因此通过星载激光测高仪回波脉宽反演海面风速，与目前普遍应用的由微波雷达测高仪的雷达散射截面反演风速方法相比，具有更高的脚点定位精度和空间分辨率。
文档编号G01S7/48GK102759731SQ20121025996
公开日2012年10月31日 申请日期2012年7月25日 优先权日2012年7月25日
发明者周辉, 李松, 葛兵, 郑国兴, 马跃 申请人:武汉大学