CN102759731B - 基于星载激光测高仪回波的海洋表面风、浪特性反演方法 - Google Patents

Abstract

本发明涉及基于星载激光测高仪回波的海洋表面风、浪特性反演方法。本发明首先对星载激光测高仪回波波形进行预处理和特征提取，再建立提取的回波参数与海洋表面的风、浪场之间的数学关系，即可根据回波参数来反演计算海洋表面波浪斜率、波高和海洋表面风速。由于激光光束具有远小于微波波束的发散角的特点，因此通过星载激光测高仪回波脉宽反演海面风速，与目前普遍应用的由微波雷达测高仪的雷达散射截面反演风速方法相比，具有更高的脚点定位精度和空间分辨率。

Description

基于星载激光测高仪回波的海洋表面风、浪特性反演方法

技术领域

本发明涉及一种海洋表面风、浪特性反演方法，尤其是涉及一种基于星载激光测高仪回波的海洋表面风、浪特性反演方法。

背景技术

海洋表面波浪斜率是重要的海态参数，同时波浪斜率还与海面上空风场强相关。目前对海洋表面波浪斜率的遥感测量主要由机载或星载微波雷达进行：根据雷达反射回波计算目标散射截面，并由此反演海面平均斜率(MSS,mean square slope)，进而通过风浪模型推算海平面上方风速。星载激光测高仪是一种安置于低轨卫星平台，通过发射激光脉冲、接收激光脉冲经地物反射后的微弱回波，由两者之间的时间差以及卫星平台的姿态、位置计算地物高程的新型卫星载荷。对地观测激光测高仪如果具有全波形记录能力，则可以进一步演变成同时具备测距和反演目标表面特性的激光雷达遥感系统。从上世纪70年代开始用于海洋表面测量的星载微波雷达高度计兴起并逐渐成熟，成为全球范围内海洋表面波浪形态观测的主要手段。相对于微波雷达测高仪，激光测高仪具有更小的光束发散角和工作波长、更精确的空间定位能力和更高的测量精度，并具备使用波形参数来反演被测目标信息的功能，是对地高精度测量与遥感的重要发展方向。世界上首台具有全波形记录能力的对地观测星载激光测高仪是美国2003年发射的ICESat卫星上搭载的地球科学激光测高系统--GLAS系统。该系统在轨运行的3年中，记录下了包括陆地、海洋、海冰和极区冰盖等不同地物的海量回波波形。由于GLAS系统主要用作对极区（主要是南北极）冰盖的测量，从回波形态中分析海洋表面特性、反演海洋表面风速的工作成果少见报道。

发明内容

本发明主要是解决现有技术所存在的技术问题；提供了一种基于星载激光测高仪回波的海洋表面风、浪特性反演方法。

本发明的上述技术问题主要是通过下述技术方案得以解决的：

一种基于星载激光测高仪回波的海洋表面风、浪特性反演方法，其特征在于，包括以下步骤：

步骤1，设定采样点数以及窗口长度，利用GLAS测高仪系统获取海洋回波波形信号后，计算信号标准差δ_s；并计算回波白噪声标准差δ_n；根据得到的δ_s和δ_n获取高斯滤波器宽度σ；

步骤2，由所得高斯滤波器宽度σ通过连续型高斯滤波器基于高斯滤波函数进行采样并归一化处理，得出离散高斯滤波器参数h_k；

步骤3，根据步骤3中得到的数字滤波器参数h_k获取滤波后信号；

步骤4，对步骤3中滤波后信号做参数提取，得到单高斯回波的峰值A_i，脉宽σ_i和重心点T_i位置；

步骤5，根据步骤4中得到的脉宽σ_i基于测高仪系统参数测高仪轨道高度z，光束发散角θ_T，发射脉冲宽度σ_f，接收系统展宽σ_h和真空光速c，获取海平面上方风速w。

在上述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，所述的步骤1中，回波白噪声标准差δ_n的具体获取方法如下：

a.计算接收离散波形S_k所有点的纵坐标幅度均值S_mean；

b.计算从波形尾端小于S_mean的第一个点起始到波形末尾点集合的标准差作为δ_n；

3.根据权利要求2所述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，其特征在于，所述的步骤3中，滤波后信号的获取是根据对数字滤波器参数h_k与回波离散信号S_k作卷积运算得到。

在上述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，所述的步骤4中，参数提取的具体方法如下：

步骤4.1.对步骤3中滤波后信号求二阶差分，找出全部拐点；

步骤4.2.针对步骤4.1中的全部拐点获取满足条件的n个拐点对：满足条件：波形中第2i-1个是正向直线外点，且第2i个是正向直线内点的拐点对，其中i=1,2,…n。

步骤4.3.在第i组拐点对区间内中找出最大值作为预估计高斯振幅A_{i_pre}，两拐点均值为预估计高斯重心位置T_{i_pre}，找出区间内信号值等于0.8倍预估高斯振幅两侧点位置，其两点间距为预估计高斯脉宽σ_{i_pre}；

步骤4.4.以估计的高斯参数为初始波形，基于梯度下降法或梯度下降-高斯牛顿法的LM算法对滤波后信号进行高斯拟合，当达到设定拟合收敛条件时，拟合结束得出各个高斯分量的拟合参数A_i，T_i和σ_i。

在上述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，所述的步骤1中，高斯滤波器宽度σ基于公式（式15）：

$\mathrm{\&sigma;}=<k\&CenterDot;\frac{{\mathrm{\&delta;}}_n}{{\mathrm{\&delta;}}_s}>$

其中k为常正整数，〈〉为最近取整函数，δ_n为权利要求2中计算的回波白噪声标准差，δ_s为信号标准差，可以通过计算回波信号当前点与窗口长度为j的前后2j+1个回波信号点的标准差获得，通常情况下长度2j+1约为5%左右信号总采样点长度。

在上述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，所述高斯滤波函数为：

$h\left(t\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}}\mathrm{\&sigma;}}\exp \left(\frac{{-t}^2}{{2\mathrm{\&sigma;}}^2}\right)$

其中，σ为滤波器宽度，exp()为指数函数，t为横坐标时间轴。

在上述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，所述的步骤4.4中，定义收敛条件为ε≤δn，基于梯度下降法或梯度下降-高斯牛顿法的LM算法对滤波后信号进行高斯拟合，当拟合满足条件：

时，拟合结束，其中S_k为离散激光高度计测高仪回波信号，L为回波信号的采样点数，w_k为式(12)的离散采样点，其采样间隔与回波信号相同。其中，LM算法是通用算法，输入初始参数，预拟合函数类型和拟合结束条件后，其会运算自收敛或无法收敛；这里输入初始参数是4.1-4.3步骤中计算的A_{i_pre}，T_{i_pre}和σ_{i_pre}；拟合收敛条件为ε≤δn；预拟合函数类型为12式，拟合结果是得出12式中的参数A_i，T_i和σ_i。

在上述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，还包括一个步骤6：即由回波波形反演得到风速数据之后，利用海面均方斜率公式s²=0.003+0.00512w能够获取海面平均斜率；根据风浪波高公式SWH=4σ_ξ=0.064w²能够获取波浪高度。

因此，本发明具有如下优点：因此，本发明主要有以下优点：1）激光测高仪光束发散角通常为亚毫弧量级，远小于雷达测高仪微波波束发散角的特点，具有更高的脚点定位精度和空间分辨率。2）由于使用回波的脉冲宽度反演风速，与回波整体能量无关，不受大气衰减变化的影响，具有较强的通用性和抗干扰性。3）与目前通用的雷达测高仪使用雷达散射截面反演风速的原理和途径完全不同，可以与雷达测高仪的测量数据进行校验或补充。

附图说明

图1是本发明中天顶方向入射的激光测高仪探测海洋表面回波的工作示意图。

图2是本发明应用GLAS系统发射的激光脉冲时间波形，其中，横轴为相对时间轴，单位为ns；纵轴为幅值，单位为V，援引GLA01 i_tx_wf V。

图3是本发明中GLAS海洋回波处理结果，其中，实线表示海洋回波的原始波形GLA01 i_rng_wf V，虚线为经过自适应高斯滤波和高斯拟合后波形，具体为图6中第一组海洋回波数据。

图4是本发明GLAS公布的海洋回波高斯滤波处理后的波形结果，与图3处理结果对比有很好的一致性，两者拟合具体参数在图6中给出。

图5是2003年2月21日零时NCEP公布的全球表面上方12m位置的风速信息，NCEP提供全球数值天气分析数据集，数据集在格林尼治时间00:00、06:00、12:00和18:00时更新，其中，包含地表大气压、风速、可降水量和以标准大气压层给出的温度、相对湿度和位势高度数据。其全球坐标纬度为北向为正，南向为负；经度为从0度经线向东为正到359度，360度经线与0度重合；其中风速数据援引NCEP第280行，281行的经向风和纬向风的均方根值。

图6是使用GLAS海洋回波计算波形参数结果与GLAS官方处理公布结果对比。高斯拟合后的参数包括峰值，重心位置和脉冲宽度，其中预处理后提取出的脉冲宽度用于海平面上方风速反演计算。表中UTC时间为以2000年1月1日12时为基准0点的相对时间，实际为波形编号；重心位置为高斯重心与最后一个回波点的时间差值，因此为负数。使用该方法处理回波波形误差均小于1%。表1中所有海洋回波均为2003年2月21日00:13分到00:40分数据，与图5中NCEP气象数据时间接近，便于减小误差；GLAS官方数据是援引GLA05 d_para2[01][02][03]中分别给出Std Fit Peak 1Amplitude Sigma(V)，Std Fit Peak 1 Location Sigma(ns)，Std fit Peak1 Width Sigma(ns)，即波形数据处理结束后计算出的回波信号高斯波形的振幅，重心位置和均方根脉宽。

图7是反演海洋表面上方风速数据与NCEP公布风速对比。表1中海洋回波通过计算处理后的脉冲宽度使用10式反演海平面上方风速w，其中GLAS系统参数为轨道高度z=600km，光束发散角θ_T=110urad，采样间隔Δt=1ns，对应系统展宽σ_h=1/12^0.5ns，发射脉宽σ_f使用同样方法做滤波拟合提取高斯参数；通过对比结果反演风速准确度很好。

图8是海洋表面斜率和波高计算数据。通过表2反演风速结果，根据Cox-Munk海面均方斜率公式s²=0.003+0.00512w计算海面平均斜率；根据Phillips风浪波高公式SWH=4σ_ξ＝0.064w²计算波浪高度(Significant waveheight)，其中w为风速，σ_ξ为RMS波高。

具体实施方式

下面通过实施例，并结合附图，对本发明的技术方案作进一步具体的说明。

实施例：

首先介绍一下本发明所需要的理论基础：

1.回波脉宽与海面风速、波浪斜率之间关系的理论推导

根据Gardner的理论，发射激光脉冲经过一次菲涅尔衍射，入射到地球表面经过目标反射，再经过一次菲涅尔衍射，到达接收望远镜视场，望远镜接收到回波功率由(1)式表示，式中ρ为反射目标横截面坐标矢量，β(ρ)光斑内目标反射率，ξ(ρ)目标表面高程轮廓，A_R为接收望远镜面积，z为测高仪轨道高度，c为真空光速，Σ目标光斑照射区域，a(ρ,z)和f(t)分别为激光脉冲在空域和时域分布函数。

$P\left(t\right)=\frac{A_R{T_a}^2}{a^2}{\&Integral;\&Integral;}_{\mathrm{\&Sigma;}}\mathrm{\&beta;}\left(\mathrm{\&rho;}\right)\&CenterDot;{\left|a(\mathrm{\&rho;},z)\right|}^2\&CenterDot;{\left|f\right[t-\frac{2z}{c}-\frac{{\mathrm{\&rho;}}^2}{\mathrm{cz}}+\frac{2\mathrm{\&xi;}\left(\mathrm{\&rho;}\right)}{c}\left]\right|}^2{d}^2\mathrm{\&rho;}---\left(1\right)$

通常情况下，激光发射脉冲在空域和时域分布都可近似为高斯函数，分别由(2)式和(3)式表示。式中Q为发射激光脉冲的能量，θ_T为初始光束发散角，σ_f为发射脉冲RMS宽度：

${\left|a(\mathrm{\&rho;},z)\right|}^2=\frac{Q}{2\mathrm{\&pi;}{\left(z\tan {\mathrm{\&theta;}}_T\right)}^2}\exp ({-\mathrm{\&rho;}}^2/2z^2{\tan }^2{\mathrm{\&theta;}}_T)---\left(2\right)$

${\left|f\left(t\right)\right|}^2=\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&sigma;}}_f}\exp (-\frac{t^2}{2{\mathrm{\&sigma;}}_f^2})---\left(3\right)$

按照Kodis理论，在近天顶入射的条件下，海洋表面反射率由表面随机分布镜面点的后向散射决定，散射功率正比于被照亮的镜面点个数n。每个镜面点散射可以看作是一个正切球，球半径等于镜面点表面的长短轴半径r_a和r_b的几何平均。在光学领域散射截面σ=πδ〈|r_ar_b|〉，归一化到Ω=4π的球面度，并代入Barrick的结论n〈|r_ar_b|〉=(1+ρ²/z²)²p[ξx,ξy|ξ(ρ)]后，海洋表面反射率β由(4)式表示：

$\mathrm{\&beta;}\left(\mathrm{\&rho;}\right)=\frac{\mathrm{\&delta;n}<\left|r_a{r}_b\right|>}{4}=\mathrm{\&delta;}{(1+\frac{{\mathrm{\&rho;}}^2}{z^2})}^{2}p[{\mathrm{\&xi;}}_x,{\mathrm{\&xi;}}_y|\mathrm{\&xi;}\left(\mathrm{\&rho;}\right)]/4---\left(4\right)$

(4)式中δ是海水与大气之间镜面反射率，与入射激光波长λ有关；p[ξ_x,ξ_y|ξ(ρ)]为海面斜率在波高为ξ时的条件概率密度，ξ_x和ξ_y分别表示光斑照明区域海水表面横截面x和y方向波浪斜率。根据Cox-Munk经验公式，在忽略波浪毛细波和地球自转等影响的理想情况下，波高p(ξ)和波浪表面轮廓p(ξ_x,ξ_y)满足高斯分布，且仅与风速有关。由于波浪同一点处的高度和斜率是不相关的，则波高和波浪斜率的联合概率密度函数可由(5)式表示：

$p[{\mathrm{\&xi;}}_x,{\mathrm{\&xi;}}_y|\mathrm{\&xi;}\left(\mathrm{\&rho;}\right)]=p\left(\mathrm{\&xi;}\right)p({\mathrm{\&xi;}}_x,{\mathrm{\&xi;}}_y)=\frac{1}{\sqrt{2\mathrm{\&pi;}}{\mathrm{\&sigma;}}_{\mathrm{\&xi;}}}\exp \left(\frac{-{\mathrm{\&xi;}}^2}{2{\mathrm{\&sigma;}}_{\mathrm{\&xi;}}^2}\right)\&CenterDot;\frac{1}{{\mathrm{\&pi;s}}^2}\exp (-\frac{{\mathrm{\&xi;}}_x^2+{\mathrm{\&xi;}}_y^2}{s^2})---\left(5\right)$

式(5)中σ_ξ=0.016w²为Phillips的风浪波高谱，s²=0.003+0.00512w为Cox-Munk海面均方斜率，w为在海平面上方12.5米处的平均风速。将式(2)-(5)代入(1)式得到测高仪探测器输出的回波信号强度的时间分布函数S(t)由(6)式表示为：

$S\left(t\right)=\frac{\mathrm{\&eta;}}{\mathrm{h\&upsi;}}\&CenterDot;\frac{\mathrm{\&delta;}{A}_R{T_a}^2}{4z^2}{\&Integral;\&Integral;}_{\mathrm{\&Sigma;}}{\left|a(\mathrm{\&rho;},z)\right|}^2{(1+\frac{{\mathrm{\&rho;}}^2}{z^2})}^2\&Integral;p({\mathrm{\&xi;}}_x,{\mathrm{\&xi;}}_y,\mathrm{\&xi;})\mathrm{d\&xi;}{\left|f\right[t-\frac{2z}{c}-\frac{{\mathrm{\&rho;}}^2}{\mathrm{cz}}+\frac{2\mathrm{\&xi;}}{c}\left]\right|}^2*h\left(t\right)d^2\mathrm{\&rho;}---\left(6\right)$

(6)式中h为普朗克常数，v为光波频率，η为接收系统光子效率，(6)式中卷积g(t)=|f(t)|²*h(t)式中，h(t)表示发射脉冲，g(t)是接收系统的脉冲响应；对(1)式回波功率的不同阶次时间积分分别可以得出平均回波光子数N，回波重心T和回波均方脉宽σ，其中回波重心T和回波脉宽σ分别由(7)式和(8)式表示：

$T_i=\frac{\underset{-\&infin;}{\overset{+\&infin;}{\text{\&Integral;}}}\mathrm{tS}\left(t\right)\mathrm{dt}}{\underset{-\&infin;}{\overset{+\&infin;}{\&Integral;}}S\left(t\right)\mathrm{dt}}=\frac{2z}{c}+\frac{{\mathrm{\&sigma;}}_r^2}{\mathrm{cz}}-\frac{2<\mathrm{\&xi;}>}{c}---\left(7\right)$

${\mathrm{\&sigma;}}_i^2=\frac{\underset{-\&infin;}{\overset{+\&infin;}{\&Integral;}}{(t-T_i)}^{2}S\left(t\right)\mathrm{dt}}{\underset{-\&infin;}{\overset{+\&infin;}{\&Integral;}}S\left(t\right)\mathrm{dt}}={\mathrm{\&sigma;}}_h^2+{\mathrm{\&sigma;}}_f^2+\frac{4}{c^2}({\mathrm{\&sigma;}}_{\mathrm{\&xi;}}^2-{<\mathrm{\&xi;}>}^2)+{\left(\mathrm{cz}\right)}^{-2}\&Integral;{({\mathrm{\&rho;}}^2-{\mathrm{\&sigma;}}_r^2)}^2{b}_2(\mathrm{\&rho;},z)d^2\mathrm{\&rho;}---\left(8\right)$

其中， $b_2(\mathrm{\&rho;},z)={\left|a_i(\mathrm{\&rho;},z)\right|}^2{(1+\frac{{\mathrm{\&rho;}}^2}{z^2})}^{2}p({\mathrm{\&xi;}}_x,{\mathrm{\&xi;}}_y)/\&Integral;{\left|a_i(\mathrm{\&rho;},z)\right|}^2{(1+\frac{{\mathrm{\&rho;}}^2}{z^2})}^{2}p({\mathrm{\&xi;}}_x,{\mathrm{\&xi;}}_y)d^2\mathrm{\&rho;},$

对于波高满足高斯分布的波浪而言<ξ>=0。

由于(7)式和(8)式中对ρ的积分仍然有波浪斜率分布p(ξ_x,ξ_y)，需要将其转换为与ρ有关的函数。假设海面斜率分布为一维形式，如图1所示海面镜面点只有当该点斜率等于天顶角正切值情况下产生的后向散射才能被望远镜接收，即tanθ=ρ/z。对ρ积分得到(9)式：

${\mathrm{\&sigma;}}_r^2=\frac{{\&Integral;}_0^{\&infin;}{(1+\frac{{\mathrm{\&rho;}}^2}{z^2})}^2\exp \left(\frac{{-\mathrm{\&rho;}}^2}{2z^2{\tan }^2{\mathrm{\&theta;}}_T}\right)\exp \left(\frac{-{\mathrm{\&rho;}}^2}{z^2{s}^2}\right){\mathrm{\&rho;}}^3\mathrm{d\&rho;}{\&Integral;}_0^{2\mathrm{\&pi;}}\mathrm{d\&theta;}}{{\&Integral;}_0^{\&infin;}{(1+\frac{{\mathrm{\&rho;}}^2}{z^2})}^2\exp \left(\frac{{-\mathrm{\&rho;}}^2}{2z^2{\tan }^2{\mathrm{\&theta;}}_T}\right)\exp \left(\frac{-{\mathrm{\&rho;}}^2}{z^2{s}^2}\right)\mathrm{\&rho;d\&rho;}{\&Integral;}_0^{2\mathrm{\&pi;}}\mathrm{d\&theta;}}=\frac{{2z}^2}{({\tan }^{-2}{\mathrm{\&theta;}}_T+2s^{-2})}---\left(9\right)$

将(9)式代入(8)式，并代入Phillips风浪波高谱和Cox-Munk海面均方斜率表达式，可得激光测高仪探测到的回波脉宽表达式为：

${\mathrm{\&sigma;}}_i^2={\mathrm{\&sigma;}}_h^2+{\mathrm{\&sigma;}}_f^2+\frac{4}{c^2}{\mathrm{\&sigma;}}_{\mathrm{\&xi;}}^2+\frac{{\mathrm{\&sigma;}}_r^4}{c^2{z}^2}={\mathrm{\&sigma;}}_h^2+{\mathrm{\&sigma;}}_f^2+\frac{4}{c^2}{\left(0.016w^2\right)}^2+\frac{4z^2}{c^2}{({\tan }^{-2}{\mathrm{\&theta;}}_T+2/(0.003+0.00512w))}^{-2}---\left(10\right)$

(10)式显示：测高仪回波脉宽仅与风速w及若干系统参数相关，可以利用回波波形的脉宽来计算风速；由回波波形反演得到风速数据之后，利用Cox-Munk海面均方斜率公式s²=0.003+0.00512w还可计算海面平均斜率(Mean square slope)；利用Phillips风浪波高公式SWH=4σ_ξ=0.064w²计算波浪高度(Significant wave height)。

2.回波信号噪声抑制理论

激光测高仪通过接收和处理激光发射脉冲经大气传输在目标表面反射后的回波波形数据，得到测距结果，并可反演目标表面特性。一般假设测高仪初始发射的激光脉冲的时间和空间分布均满足高斯分布，在不考虑大气散射效应对脉冲影响（仅考虑大气衰减）情况下，激光发射脉冲经过一次菲涅尔衍射入射到地表目标；经由目标表面的反射、再经过一次菲涅尔衍射，进入望远镜视场。由于星载激光测高仪在地面的光斑直径往往在几十米的量级，因此回波信号通常为多模式复杂曲线，由混入噪声的若干高斯分量组成；而海洋回波通常情况下则为近似单高斯型回波。

一维归一化高斯滤波函数表示为式(11)，其中σ为高斯函数宽度。

$h\left(t\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}}\mathrm{\&sigma;}}\exp \left(\frac{{-t}^2}{2{\mathrm{\&sigma;}}^2}\right)---\left(11\right)$

由m个高斯分量组成的有效回波w(t)由式(12)表示，其中A_i，σ_i和T_i分别为第i个高斯分量的幅度，宽度和位置，实际激光测高仪回波s(t)混入了白噪声n(t)，即s(t)=w(t)+n(t)。

$w\left(t\right)=\underset{1}{\overset{m}{\mathrm{\&Sigma;}}}{A}_i\exp \left[\frac{-{(t-T_i)}^2}{{{2\mathrm{\&sigma;}}_i}^2}\right]---\left(12\right)$

高斯滤波过程等价于信号s(t)与滤波函数g(t)的卷积，由式(13)表示：

$s\left(t\right)*h\left(t\right)=n\left(t\right)*h\left(t\right)+\underset{1}{\overset{n}{\mathrm{\&Sigma;}}}\frac{A_i{\mathrm{\&sigma;}}_i}{\sqrt{({\mathrm{\&sigma;}}^2+{{\mathrm{\&sigma;}}_i}^2)}}\exp \left[\frac{-{(t-T_i)}^2}{2({\mathrm{\&sigma;}}^2+{{\mathrm{\&sigma;}}_i}^2)}\right]---\left(13\right)$

通过式(13)可以得出，与高斯滤波函数卷积后各个高斯分量宽度σ_i’=(σ²+σ_i ²)^0.5，峰值A_i’=σA_i/(σ²+σ_i ²)^0.5，即滤波后不仅降低高斯信号峰值，也展宽了高斯信号脉宽。这里使用的测高仪回波是经光电和模数转换后数字化仪输出的离散信号，减去上一节方法估计的噪声均值，得到含有零均值白噪声的信号。离散高斯滤波器参数序列h_k由式(3)采样并归一化获得，与信号等采样间隔；GLAS数字化仪采样间隔bin为1ns，2ns和4ns3种情况，绝大多数为1ns；满足涵盖99%能量的3σ原则采样点数N≥6σ+1。假设白噪声方差为δ_n ²，有效信号方差为δ_f ²，回波信号整体方差为δ_s ²=δ_f ²+δ_n ²；滤波后方差δ₀ ²由式(14)表示。

${\mathrm{\&delta;}}_0^2=\frac{{\mathrm{\&delta;}}_s^2}{2\sqrt{\mathrm{\&pi;}}\mathrm{\&sigma;}}\&GreaterEqual;\frac{1.91{\mathrm{\&delta;}}_s^2}{N}={\mathrm{\&delta;}}_s^2\&CenterDot;{\mathrm{\&Sigma;}}_{k=1}^N{h}_k^2---\left(14\right)$

通过式(14)可以得出宽度σ越大，滤波参数序列N越长，对噪声的抑制效果越好；但式(5)表明这将导致有效高斯分量变形，结合上一节的分析，在背景区域没有有效信号即δ_f ²=0，回波信号方差δ_s ²=δ_n ²，只与噪声方差有关；而在有效信号区域δ_f ²>>δ_n ²，δ_s ²≈δ_f ²回波信号方差主要由有效信号决定，如果在不同信号方差区域使用滤波宽度σ随之变化的方法，在背景区域滤波宽度σ较大，在有效信号区域σ很小，则可以抑制大部分噪声并最大限度的保留有效信号的完整性。基于上述原则设计高斯滤波自适应宽度σ由式(15)决定：

$\mathrm{\&sigma;}=<k\&CenterDot;\frac{{\mathrm{\&delta;}}_n}{{\mathrm{\&delta;}}_s}>=<k\&CenterDot;\sqrt{\frac{{\mathrm{\&delta;}}_n^2}{{\mathrm{\&delta;}}_f^2+{\mathrm{\&delta;}}_n^2}}>---\left(15\right)$

其中k为常正整数，〈〉为最近取整函数，由式(15)所给出高斯滤波器宽度σ在背景区域σ=k；在有效信号区域σ=〈kδ_n/δ_s〉，而由于此时δ_s ²≈δ_f ²>>δ_n ²，滤波器宽度σ将随有效信号的尖锐程度变化，有效信号越尖锐σ越小，甚至趋近于0；而在有效信号较弱的过度区域滤波器宽度σ将介于上面两种情况之间；这不仅保证了在背景区域可以使用相同的最大k bin宽度的高斯滤波器滤除噪声，也使得在有效信号区域滤除噪声的同时最小限度的破坏有效信号波形；而由于使用数字滤波器，滤波器采样点数N=6σ+1，因此对σ采取最近取整，且规定σ最小取值为1bin。使用式(15)计算自适应宽度σ中需要有已知以下两点：1)信号标准差δ_s，δ_s可以通过计算回波信号当前点与窗口长度为j的前后2j+1个回波信号点的标准差获得，通常情况下长度2j+1约为5%左右信号总采样点长度。2)回波的白噪声标准差δ_n计算方法如流程A所示：

a.计算接收离散波形S_k所有点的纵坐标幅度均值S_mean。

b.计算从波形尾端小于S_mean的第一个点起始（从此点开始以后点对应波形幅度都小于幅度均值S_mean）到波形末尾点集合的标准差作为δ_n。

3.高斯拟合方法

经过滤波后的回波信号需要经过高斯拟合计算出回波中含有一个或多个高斯波形的峰值A_i，脉宽σ_i和重心点T_i位置，用于反演地表目标信息。高斯拟合主要经过流程B所示：

a.对滤波后信号求二阶差分，找出全部拐点；

b.求出满足条件：波形中第2i-1个是正向直线外点，且第2i个是正向直线内点的拐点对；

c.在第i组拐点对区间内中找出最大值作为预估计高斯振幅A_{i_pre}，两拐点均值为预估计高斯重心位置T_{i_pre}，找出区间内信号值等于0.8倍预估高斯振幅两侧点位置，其两点间距为预估计高斯脉宽σ_{i_pre}；

d.以估计的高斯参数为初始波形，基于最小二乘LM算法对滤波后信号进行高斯拟合，当达到设定拟合收敛条件时，拟合结束得出各个高斯分量的拟合参数A_i，T_i和σ_i。

LM算法拟合的目标是使拟合数据与实际数据之间的余量尽量接近背景白噪声，拟合数据与实际数据的余量的标准偏差ε为(16)式表示，S_k为离散激光高度计测高仪回波信号，L为回波信号的采样点数，w_k为式(12)的离散采样点，其采样间隔与回波信号相同。

$\mathrm{\&epsiv;}=\sqrt{{\mathrm{\&Sigma;}}_{k=1}^L\frac{{(S_k-w_k)}^2}{L}}---\left(16\right)$

通常采用的判据为余量的标准偏差ε应小于2倍的噪声标准差，即ε<2δ_n。如果某一轮迭代退出后，ε不满足该判据，则可能是一些较小的有效回波分量被当成噪声滤去了，应将这些分量重新计入初始模型，进行新一轮的LM算法迭代。

4.实际反演计算流程

由于激光测高仪海洋回波通常情况下为近似单高斯型回波，在噪声抑制和高斯参数提取中m=1，即只有一个峰值A_i，脉宽σ_i和重心点T_i位置。具体步骤如下流程：

a.GLAS测高仪海洋回波采样点数L=200，因此取j=5逐点计算信号标准差δ_s；根据流程A计算回波白噪声标准差δ_s；得出δ_s和δ_s后，根据式(7)计算高斯滤波器宽度σ。

b.由所得高斯滤波器宽度σ代入式(11)，即连续型高斯滤波器，由于要处理的波形为离散数字信号，滤波器采样点数N=6σ+1，滤波器采样间隔与回波信号采样间隔相等，对式(11)采样并归一化处理，得出h_k。

c.所得数字滤波器参数h_k与回波离散信号S_k作卷积运算，得出滤波后信号。

d.按照流程B对滤波后信号做参数提取，其中收敛条件为ε≤δn_n，当拟合满足(16)式条件时，拟合结束得出单高斯回波的峰值A_i，脉宽σ_i和重心点T_i位置。

e.得出的脉宽σ_i根据式(10)，并代入测高仪系统参数测高仪轨道高度z，光束发散角θ_T，发射脉冲宽度σ_f，接收系统展宽σ_h和真空光速c，计算海平面上方风速w。

f.由回波波形反演得到风速数据之后，利用Cox-Munk海面均方斜率公式s²=0.003+0.00512w还可计算海面平均斜率(Mean square slope)；利用Phillips风浪波高公式SWH=4σ_ξ＝0.064w²计算波浪高度(Significantwave height)。

本文中所描述的具体实施例仅仅是对本发明精神作举例说明。本发明所属技术领域的技术人员可以对所描述的具体实施例做各种各样的修改或补充或采用类似的方式替代，但并不会偏离本发明的精神或者超越所附权利要求书所定义的范围。

Claims (7)

1.一种基于星载激光测高仪回波的海洋表面风、浪特性反演方法，其

特征在于，包括以下步骤：

步骤1，设定采样点数以及窗口长度，利用GLAS测高仪系统获取海洋回波波形信号后，计算信号标准差δ_s；并计算回波白噪声标准差δ_n；根据得到的δ_s和δ_n获取高斯滤波器宽度σ；所述高斯滤波函数为：

$h\left(t\right)=\frac{1}{\sqrt{2\mathrm{\&pi;}}\mathrm{\&sigma;}}\exp \left(\frac{{-t}^2}{2{\mathrm{\&sigma;}}^2}\right)$

其中，σ为滤波器宽度，exp()为指数函数，t为横坐标时间轴；

步骤2，由所得高斯滤波器宽度σ通过连续型高斯滤波器基于高斯滤波函数进行采样并归一化处理，得出离散高斯滤波器参数序列h_k；

步骤3，根据步骤2中得到的离散高斯滤波器参数序列h_k获取滤波后信号；

步骤4，对步骤3中滤波后信号做参数提取，得到单高斯回波的峰值A_i，脉宽σ_i和重心点T_i位置；

步骤5，根据步骤4中得到的脉宽σ_i，基于测高仪系统参数：测高仪轨道高度z，光束发散角θ_T，发射脉冲宽度σ_f，接收系统展宽σ_h和真空光速c，获取海平面上方风速w,计算海平面上方风速w公式如下：

${\mathrm{\&sigma;}}_i^2={\mathrm{\&sigma;}}_h^2+{\mathrm{\&sigma;}}_f^2+\frac{4}{c^2}{\mathrm{\&sigma;}}_{\mathrm{\&xi;}}^2+\frac{{\mathrm{\&sigma;}}_r^4}{c^2{z}^2}={\mathrm{\&sigma;}}_h^2+{\mathrm{\&sigma;}}_f^2+\frac{4}{c^2}{\left(0.016w^2\right)}^2+\frac{{4z}^2}{c^2}{({\tan }^{-2}{\mathrm{\&theta;}}_T+2/(0.003+0.00512w))}^{-2}$ 。

2.根据权利要求1所述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，其特征在于，所述的步骤1中，回波白噪声标准差δ_n的具体获取方法如下：

a.计算接收离散波形S_k所有点的纵坐标幅度均值S_mean；

b.计算从波形尾端小于S_mean的第一个点起始到波形末尾点集合的标准差作为δ_n。

3.根据权利要求2所述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，其特征在于，所述的步骤3中，滤波后信号的获取是根据对离散高斯滤波器参数序列h_k与回波离散信号S_k作卷积运算得到。

4.根据权利要求1所述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，其特征在于，所述的步骤4中，参数提取的具体方法如下：

步骤4.1.对步骤3中滤波后信号求二阶差分，找出全部拐点；

步骤4.2.针对步骤4.1中的全部拐点获取满足条件的n个拐点对：满足条件：波形中第2i-1个是正向直线外点，且第2i个是正向直线内点的拐点对，其中i=1,2,…n；

步骤4.3.在第i组拐点对区间内中找出最大值作为预估计高斯振幅A_{i_pre}，两拐点均值为预估计高斯重心位置T_{i_pre}，找出区间内信号值等于0.8倍预估高斯振幅两侧点位置，其两点间距为预估计高斯脉宽σ_{i_pre}；

步骤4.4.以估计的高斯参数为初始波形，基于梯度下降法或梯度下降-高斯牛顿法的LM算法对滤波后信号进行高斯拟合，当达到设定拟合收敛条件时，拟合结束得出各个高斯分量的拟合参数A_i，T_i和σ_i。

5.根据权利要求2所述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，其特征在于，所述的步骤1中，高斯滤波器宽度σ基于公式：

$\mathrm{\&sigma;}=<k\&CenterDot;\frac{{\mathrm{\&delta;}}_n}{{\mathrm{\&delta;}}_s}>$

其中k为常正整数，〈〉为最近取整函数，δ_n为权利要求2中计算的回波白噪声标准差，δ_s为信号标准差，可以通过计算回波信号当前点与窗口长度为j的前后2j+1个回波信号点的标准差获得，通常情况下长度2j+1约为5%左右信号总采样点长度。

6.根据权利要求4所述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，其特征在于，所述的步骤4.4中，定义收敛条件为ε≤δ_n，基于梯度下降法或梯度下降-高斯牛顿法的LM算法对滤波后信号进行高斯拟合，当拟合满足条件：

时，拟合结束，其中S_k为离散激光高度计测高仪回波信号，L为回波信号的采样点数，w_k为式十二的离散采样点，其采样间隔与回波信号相同；

$w\left(t\right)=\underset{1}{\overset{m}{\mathrm{\&Sigma;}}}{A}_i\exp \left[\frac{{-(t-T_i)}^2}{2{{\mathrm{\&sigma;}}_i}^2}\right]$ 式十二。

7.根据权利要求1所述的基于星载激光测高仪回波的海洋表面风、浪特性反演方法，其特征在于，还包括一个步骤6：即由回波波形反演得到风速数据之后，利用海面均方斜率公式s²=0.003+0.00512w能够获取海面平均斜率；根据风浪波高公式SWH=4σ_ξ=0.064w²能够获取波浪高度。

Images