机载合成孔径雷达高度计高程参数贝叶斯估计

杨磊; 周弘昊; 黄博; 廖仙华; 夏亚波

doi:10.11999/JEIT220322

机载合成孔径雷达高度计高程参数贝叶斯估计

doi: 10.11999/JEIT220322

1.
中国民航大学天津市智能信号与图像处理重点实验室天津 300300
2.
中国工程物理研究院电子工程研究所绵阳 621999

基金项目: 国家自然科学基金(61601470)，天津市自然科学基金(16JCYBJC41200)

详细信息

作者简介:
杨磊：男，副教授，研究方向为高分辨SAR成像及机器学习理论应用

周弘昊：男，硕士生，研究方向为机载雷达高度表系统及信号处理

黄博：男，博士生，研究方向为雷达高度表系统及信号处理

廖仙华：男，硕士生，研究方向为高分辨SAR成像及统计采样技术应用

夏亚波：男，硕士生，研究方向为高分辨SAR成像及统计采样技术应用

通讯作者:
黄博　vick123y@163.com

中图分类号: TN953
计量
- 文章访问数: 632
- HTML全文浏览量: 262
- PDF下载量: 87
- 被引次数: 1
出版历程
- 收稿日期: 2022-03-25
- 修回日期: 2022-05-24
- 网络出版日期: 2022-06-01
- 刊出日期: 2023-04-10

Elevation Estimation for Airborne Synthetic Aperture Radar Altimetry Based on Parameterized Bayesian Learning

1.
Tianjin Key Laboratory for Advanced Signal Processing, Civil Aviation University of China, Tianjin 300300, China
2.
Institute of Electronic Engineering, Chinese Academy of Engineering Physics, Mianyang 621999, China

Funds: The National Natural Science Foundation of China (61601470), The Natural Science Foundation of Tianjin (16JCYBJC41200)

摘要

摘要: 机载合成孔径雷达高度计(SARA)由于具有高航向分辨率，因此受到广泛关注。然而，现有的SARA地面高程重跟踪方法多基于最小二乘算子，高程参数估计精度和算法抑噪性能均存在上限，容易造成高程参数估计结果过拟合，对复杂高程变化适应能力有限。为此，该文提出一种基于参数化贝叶斯统计学习方法的机载SARA重跟踪算法(PR-Bayes)。通过引入目标场景地形先验概率模型，并结合模型驱动机器学习方法，可实现对目标高程信息重跟踪可信估计，从而有效避免估计参数过拟合问题。该算法基于布朗模型(BM)对SARA回波进行复杂模型参数反演，并设计哈密顿蒙特卡洛(HMC)统计采样器，实现对目标场景地形高度的参数估计。基于该文所提算法，分别通过点目标模拟和DEM半实物模拟对该算法进行有效性验证及高程参数估计精度验证，并通过实测数据验证该算法的实用性。
- 合成孔径雷达高度计 /
- 贝叶斯学习 /
- 重跟踪 /
- 布朗模型 /
- 哈密顿蒙特卡洛采样
Abstract: Airborne Synthetic Aperture Radar Altimeter (SARA) is capable of exploiting the high-resolution in along-track, which has been attracted wide concerns. However, the existing re-tracking methods are mostly based on the least square operator. The performance of estimation accuracy and noise suppression of the operator are limited due to neglect of noise factors and accordance over-fitting problem. In this paper, Parameterized Retracking Bayes (PR-Bayes) algorithm is proposed under the framework of Bayesian machine learning. By introducing a prior probability model of the terrain scene, and combining with model-driven machine learning method, the elevation information of re-tracking with reliable estimation of the target can be achieved. The problem about over-fitting can be alleviated effectively. In this algorithm, Brown Model (BM) is used to recover complicated model parameters of SARA echo. Then, Hamilton Monte Carlo (HMC) statistical sampler is designed to estimate the terrain height of the scene with a high accuracy and reliable confidence. The accuracy and validity of this algorithm are verified by point target simulation and semi-physical simulation based on DEM respectively, and the practicability is proved by the airborne raw SARA data.
- Synthetic Aperture Radar Altimeter(SARA) /
- Bayesian learning /
- Retracking /
- Brown Model (BM) /
- Hamilton Monte Carlo(HMC) sampling

HTML全文

1. 引言

机载雷达高度计对地形测绘和飞行安全均具有重要意义。传统雷达高度计主要利用有限脉冲测量地面高程参数^[1]。而新型的合成孔径雷达高度计(Synthetic Aperture Radar Altimeter, SARA)沿航向引入合成孔径思想，充分考虑因运动生成的多普勒带宽，可有效提高沿航向高程分辨率和测高精度^[1,2]。

SARA通过天线指向天底点正视的方式，获得沿航线向天底点1维的绝对高度信息^[2,3]。在技术上SARA采用了延时多普勒补偿技术，根据机载SARA沿方位向运动产生的多普勒频率，在航线形成诸多观测点对目标地形高度探测做出贡献，有效提高观测次数^[1,3]。与SAR处理方式类似，回波信号在SARA中通过去斜、A/D采集、快速傅里叶变换(Fast Fourier Transform, FFT)等处理，并通过方位多普勒信息进行延迟校正获得延时多普勒域的2维像(Delay Doppler Map, DDM)，故SARA又被称为延时多普勒测高计^[3,4]。

雷达回波信号含有回波噪声和仪器热噪声等诸多影响，故必须对回波进行一定的反演处理，才能从含有噪声的曲线中得到观测参数的估计值，这一过程即是对回波的跟踪处理^[4,5]。高度计对回波的跟踪处理有高度跟踪和重跟踪两个阶段。高度跟踪主要考虑回波信号捕获的稳健性并获得一个较为粗略的高度值。为了获得更高精度的高程参数，需在粗跟踪的基础上对回波进一步进行重跟踪处理。回波重跟踪的实质是采用一定的算法调整回波模型中的参数，使接收的回波信号与回波模型达到最佳匹配，从而反演出各项有用的高程特征参数。因此重跟踪的关键有两点，回波模型的建立以及适当的参数估计方法，即重跟踪算法^[4]。

国内外基于回波模型的重跟踪估计算法主要可分为两大类，一类结合经典最小二乘拟合算法实现对模型的参数估计^[4]，然而复杂环境下高维数据的高精度参数估计需求则易引发过拟合问题。针对该问题相关学者提出了第2类算法，在SARA模型基础上引入先验信息有效避免过拟合使得精度得到提升。一般正则先验灵活度不高，故选择引入能表征信号不确定性的概率先验，如Halimi等人^[5]通过引入统计先验的方式，利用最大后验算法进行参数求解，然而该方法只可提供点估计值，无法给出高阶统计信息，因此实用性差。

针对以上问题，本文提出一种基于贝叶斯统计机器学习方法的机载SARA回波参数化重跟踪(Parameterized Retracking Bayes, PR-Bayes)算法，使用布朗模型(Brown Model, BM)作为回波重跟踪拟合形式，为解决参数估计过拟合问题，对目标地形引入统计先验信息，并在压缩感知(Compressed Sensing, CS)的框架下，建立接收回波与目标高程参数之间的数学模型。然后，在先验建模的基础上，针对先验与似然函数不共轭的问题进行分层建模，构建目标参数与似然函数之间的联系，进而推导高程参数后验分布解析公式并进行求解。因所得后验分布较为复杂，传统基于马尔科夫链的蒙特卡洛(Markov Chain Monte Carlo, MCMC)类采样算法在此高维情况下易出现随机游走^[4,6]，因此本文引入基于哈密顿动能方程的哈密顿蒙特卡洛(Hamilton Monte Carlo, HMC)算法，通过哈密顿方程有效改善原马尔科夫链参数更新的随机性^[7,8]。对参数后验分布解析公式进行统计采样，得到该参数估计值的同时并获得此估计值如置信区间、置信度等高阶统计信息，实现复杂地形高程参数可信估计^[9,10]。本文针对所提PR-Bayes算法，进行点目标仿真、半实物模拟以及实测数据相关实验，并通过与传统最小二乘(Least Squares, LS)算法进行数据对比，验证了所提算法的有效性以及优越性。

2. 高度计回波信号模型

2.1 回波信号模型

如图1(a)所示，为获得天底点的1维绝对高度信息，SARA在传统高度计的基础上引入合成孔径思想，通过多普勒锐化处理(Doppler Beam Sharpening, DBS)，有效提高传统高度计雷达功率利用率不足问题。如图1(b)所示，根据多普勒频率对同心圆形式的脉冲足迹添加多普勒线进行细分，细分后的脉冲足迹均代表一回波信号^[2,3]。但非天底点回波信号不可直接用于高度测量，需对这些信号进行相关延迟校正。经延迟校正后的多个回波信号有效提高了观测次数，在雷达功率利用率提升的基础上进一步提高了观测精度。

图 1 合成孔径雷达高度计(SARA)工作原理示意图

下载: 全尺寸图片幻灯片

图1(a)为SARA原理几何模型，建立3维坐标系，SARA以一恒定速度 ${V_{\text{s}}}$ 沿高度为 ${h_0}$ ，方向为 $X$ 轴的航线飞行，在点G处对正下方坐标为 $\left( {{X_0},{Y_0},{Z_0}} \right)$ 的点目标T进行高度测量，其中 ${X_0}$ = ${V_{\text{s}}}{t_{\text{s}}}$ ，故SARA欲获得的高度信息 $h$ = $\left| { {{\boldsymbol{GT}}} } \right|$ = $\left| {{h_0} - {Z_0}} \right|$ ，经时延后收到如式(1)所示的回波信号

$\begin{split} {s_{\text{r}}}\left( {t,{t_{\text{s}}}} \right) =& {A_0}\exp \left\{ {\text{j}}\Bigr[ 2\pi {f_0}\left( {t - {{2R\left( {{t_{\text{s}}}} \right)}/ {\text{c}}}} \right) \right. \\ & \left. \left. + \pi {K_{\text{r}}}{{\left( {t - {{2R\left( {{t_{\text{s}}}} \right)}/ {\text{c}}}} \right)}^2} \right] \right\} \end{split}$

(1)

其中， $t$ 和 ${t_{\text{s}}}$ 别为快时间和慢时间， ${A_0}$ 为散射系数， $R$ = $\sqrt {h_{{\text{ref}}}^2 + {{\left( {{V_{\text{s}}}{t_{\text{s}}}} \right)}^2}}$ 为当前时刻SARA与观测目标之间的直线距离，选择场景中心距离 ${h_{{\text{ref}}}}$ 作为参考高度， ${f_0}$ 为中心频率， ${\text{c}}$ 为光速， ${T_p}$ 为脉冲宽度， ${K_{\text{r}}}$ 为调频率。对接收到的回波数据进行解线性频调以及RVP校正等处理^[4]，沿航向对数据通过FFT进行DBS处理可得

$\begin{split} s\left( {t,{f_{\rm{d}}}\left( m \right)} \right) = &{A_0}\exp \left\{ { - {\text{j}}{{4\pi {K_{\text{r}}}\left( {R - {h_{{\text{ref}}}}} \right)t}/ {\text{c}}}} \right\}\\ & \cdot {{\rm{sinc}}} \left\{ {{T_{\text{s}}}\left[ {{{2{X_0}{V_{\text{s}}}}/ {\left( {\lambda {h_{{\text{ref}}}}} \right)}} - {f_{\rm{d}}}\left( m \right)} \right]} \right\} \end{split}$

(2)

其中， ${f_{\rm{d}}}\left( m \right)$ 为飞机沿航线向运动时产生的第 $m$ 个多普勒频率。此时可通过多普勒频率对回波延迟进行校正补偿。设 ${f_{\rm{d}}}\left( m \right)$ 时刻SARA与目标点T 的距离为 $R(m)$ ，则延迟补偿距离量 $\Delta R$ = $\left| {h - R\left( m \right)} \right|$ ，对于此校正量采用相位补偿的方式可避免插值等复杂操作^[4]，补偿相位如式(3)所示

$\varPhi \left( {{f_{\rm{d}}}\left( m \right),t} \right) = {\text{exp}}\left\{ {{\text{j}}\frac{{4\pi {K_{\text{r}}}}}{{\text{c}}}\frac{{{\lambda ^2}{h_{{\text{ref}}}}}}{{8V_{\text{s}}^2}}f_{\rm{d}}^2\left( m \right)t} \right\}$

(3)

为此，对式(2)乘以延迟补偿相位式(3)再进行距离压缩可得多视前的延迟多普勒回波信号，如式(4)所示

$\begin{split} s\left( {f,{f_{\rm{d}}}\left( m \right)} \right) =& {A_0}{{\rm{sinc}}} \left\{ {{T_p}\left[ {f + {{2{K_{\text{r}}}\left( {h - {h_{{\text{ref}}}}} \right)} / {\text{c}}}} \right]} \right\}\\ & \cdot {{\rm{sinc}}} \left\{ {{T_{\text{s}}}\left[ {{{2{X_0}{V_{\text{s}}}}/ {\left( {\lambda {h_{{\text{ref}}}}} \right)}} - {f_{\rm{d}}}\left( m \right)} \right]} \right\} \end{split}$

(4)

其中， $f$ 为距离压缩域频率，回波经上述处理后获得DDM，进行多视处理后可获得该脉冲足迹内的点目标响应函数，如图2所示。其中，图2(a)为延迟补偿前的DDM，图2(b)为补偿后得到的DDM，对图2(b)进行多视处理得到如图2(c)所示的多视回波，其中横坐标为延迟距离门单元数。在粗跟踪高度 ${h_{{\text{ref}}}}$ 的基础上，高度计便可根据该回波函数进一步进行重跟踪处理，反演获得该点更高精度的高度数据。

图 2 2维延时多普勒像(DDM)及点目标响应函数

下载: 全尺寸图片幻灯片

2.2 回波参数模型

回波模型的建立是重跟踪算法的先决步骤，2.1节对SARA回波信号进行推导得到表达式(4)，回波模型便是对此进行相关建模。式(4)对方位向多普勒频率和高度向信息分别进行了描述，如式(4)所示。高度向和一般成像SAR的点目标响应相同，而方位向上多了一个sinc函数^[4]。故此，根据文献[4,5]回波建模，基于天线方向图相乘的思路，将方位向sinc函数作为波束锐化后的子波束方向图，将此波束方向图与物理天线方向图相乘建立SARA的布朗回波模型，获得一个关于延迟参数 $\tau$ 的函数。

现SARA回波模型大多用于粗糙海表面的回波信号，经调研，Brown^[11,12]提出布朗模型的建立是基于海表面或者具有复杂异构特征的陆地，本文尝试将布朗模型适当修改并用于部分陆地地形。传统雷达高度计回波模型为3项卷积模型，其中包含强散射体高度的概率密度函数(Probability Density Function, PDF)、雷达的点目标响应(Point Target Response, PTR)和平坦表面冲激响应(Flat Surface Impulse Response, FSIR)。而SARA在此基础上引入多普勒频率，故本文采用2维布朗模型^[5]

$s(t,{f_{\rm{d}}}\left( m \right)) = {\text{FSIR}}(t,{f_{\rm{d}}}\left( m \right)) * {\text{PDF}}\left( t \right) * {\text{PTR}}(t,{f_{\rm{d}}}\left( m \right))$

(5)

其中， $t$ = ${N_r}T$ ， ${N_r}$ 为回波距离门数，一般取 ${N_r}$ =128， $T$ 为时间分辨率。根据文献[11,12]所提，FSIR当满足一定条件时，不仅适用于海表面，同样也适用于具有异构特征的复杂地形，其表达式为

$\begin{split} {\text{FSIR}}(t,m) =& \frac{{{P_u}}}{\pi }\exp \left[ - \frac{{4{\text{c}}}}{{\gamma {h_{{\text{ref}}}}}}(t - \tau T)\right]U(t - \tau T) \\ & \times \left[{\phi _{t,m + 1}} - {\phi _{t,m}}\right]\\[-10pt] \end{split}$

(6)

其中， $\tau$ 为高度延迟距离门数， $\gamma$ 为天线波束宽度， $U(\cdot)$ 为Heaviside函数^[5]， ${\phi _{t,m}} = {{\rm{Re}}} [\arctan ({{X_m}}/ {\sqrt {R - X_m^2} })]$ ， ${X_m}$ 为多普勒频率 ${f_{\rm{d}}}\left( m \right)$ 所对应坐标。根据文献[4,5]，除绝对高程参数 $\tau$ ，当针对海洋SARA回波时，FSIR和PDF分别描述了代表了海洋风速 ${P_u}$ 和海面有效波高 ${\text{SWH}}$ 。模型参数估计方式与坐标下降算法类似，因本文重点关注绝对高程参数估计，对 $\tau$ 进行估计时， ${P_u}$ 和 ${\text{SWH}}$ 均被固定不予估计，对此下文不再赘述。

PDF通过对波高的统计性质来描述海面上同时出现高低不等的海浪波高所形成的粗糙海表面^[13,14]。而大多实际地球表面(如海洋，陆地)都可认为是粗糙的，故当高度计观察表面为均匀粗糙表面时，通常可用散射单元表面高度PDF表示粗糙度。类比于此，根据Warner^[15]推导近似表征地面回波特征的平均散射系数近似表达 ${\sigma _{\text{s}}}{\text{ = }}{{n{{\cos }^n}\left( \theta \right)} \mathord{\left/ {\vphantom {{n{{\cos }^n}\left( \theta \right)} 4}} \right. } 4}$ ，将其作为均方根粗糙度来定义适用于粗糙地面的PDF^[13,14]

${\text{PDF}}\left( t \right) = {{\exp \left( { - {{{t^2}} \mathord{\left/ {\vphantom {{{t^2}} {2\sigma _{\text{s}}^2}}} \right. } {2\sigma _{\text{s}}^2}}} \right)} \mathord{\left/ {\vphantom {{\exp \left( { - {{{t^2}} \mathord{\left/ {\vphantom {{{t^2}} {2\sigma _{\text{s}}^2}}} \right. } {2\sigma _{\text{s}}^2}}} \right)} {\sqrt {2\pi \sigma _{\text{s}}^2} }}} \right. } {\sqrt {2\pi \sigma _{\text{s}}^2} }}$

(7)

其中， $\theta$ 为入射角，n为拟合系数，通常设定为一常数。SARA 2维PTR表达式为

${\text{PTR}}(t,{f_{\rm{d}}}\left( m \right)) = {{\rm{sinc}}} {\left( {{t \mathord{\left/ {\vphantom {t T}} \right. } T}} \right)^2}{{\rm{sinc}}} {\left( {{{{f_{\rm{d}}}\left( m \right)} \mathord{\left/ {\vphantom {{{f_d}\left( m \right)} F}} \right. } F}} \right)^2}$

(8)

其中， $F$ 为频率分辨率。与SARA回波处理类似，式(5)表达的DDM也进行多视处理^[3]，得到

${{\bf{BM}}} (\tau ) = \sum\limits_m^M {s\left( {t,{f_{\rm{d}}}\left( m \right)} \right)}$

(9)

此时，可用布朗模型对式(4)回波功率函数的延迟距离门数 $\tau$ 进行估计，获得如图2(c)中半功率点所处的距离门位置，从而获得较高精度的高度值。

2.3 高程参数重跟踪模型

引入布朗模型之后，根据上文测高原理，可将SARA测高过程以矩阵形式表示，建立式(10)信号模型

${\boldsymbol{Y}}={\boldsymbol{A}}\left[{{\boldsymbol{I}}}_{M}\cdot {\bf{BM}}\left(\tau \right)\odot {{\boldsymbol{F}}}_{{\rm{d}}}\right]+{{\boldsymbol{C}}}_{N}$

(10)

其中， ${\boldsymbol{Y}} \in {\mathbb{C}^{K \times N}}$ 为SARA接收数据域回波信号，建模时需考虑降采样处理，设 $K \le M$ , $M$ , $N$ 分别为方位向和距离向采样点数， ${{\boldsymbol{C}}_N} \in {\mathbb{C}^{K \times N}}$ 为探测中的杂波或加性噪声， ${\boldsymbol{A}} \in {\mathbb{C}^{K \times M}}$ 为方位向傅里叶字典^[16,17]，式(10)中括号内部为DDM 2维像，其中 ${{\boldsymbol{I}}_M}$ 为 $M \times {\text{1}}$ 的单位向量， ${\bf{BM}}\left( \cdot \right)$ 为式(9)表示的布朗模型关于延迟距离门数 $\tau$ 的函数，故可通过 ${{\boldsymbol{I}}}_{M}\cdot{\bf{BM}}\left(\tau \right)$ 的形式来表示多视前的DDM，大小为 $M \times N$ 。 ${{\boldsymbol{F}}_{\rm{d}}} \in {\mathbb{C}^{M \times N}}$ 表示与多普勒频率相关的延迟补偿相位信息矩阵， $\odot$ 表示哈达玛(Hadamard)积，DDM经 ${{\boldsymbol{F}}_{\rm{d}}}$ 处理后可还原成延迟补偿前的延时多普勒像。故根据前文所述， ${\boldsymbol{A}}$ , ${{\boldsymbol{F}}_{\rm{d}}}$ 构建为

$\left. \begin{aligned} & {\boldsymbol{A}} = \left[ {{\boldsymbol{a}}\left( {{f_{\rm d}}\left( 1 \right)} \right),{\boldsymbol{a}}\left( {{f_{\rm d}}\left( 2 \right)} \right), \cdots ,{\boldsymbol{a}}\left( {{f_{\rm d}}\left( M \right)} \right)} \right] \\ & {\boldsymbol{a}}\left( {{f_{\rm d}}\left( m \right)} \right) = {\left\{ {\exp \left[ { - {\text{j}}2\pi {f_{\rm d}}\left( m \right){t_1}} \right],\exp \left[ { - {\text{j}}2\pi {f_{\rm d}}\left( m \right){t_2}} \right], \cdots ,\exp \left[ { - {\text{j}}2\pi {f_{\rm d}}\left( m \right){t_K}} \right]} \right\}^{\text{T}}} \end{aligned} \right\}$

(11)

$\left.\begin{aligned} & {{\boldsymbol{F}}}_{\rm d}={\left[b\left({f}_{\rm d}\left(1\right)\right),b\left({f}_{\rm d}\left(2\right)\right),\cdots ,b\left({f}_{\rm d}\left(M\right)\right)\right]}^{\text{T}}\\ & {\boldsymbol{b}}\left({f}_{\rm d}\left(m\right)\right)=\text{exp}\left[\text{j}\frac{4\pi {K}_{\text{r}}}{\text{c}}\cdot\frac{{\lambda }^{2}{h}_{\text{ref}}^{2}}{8{V}_{\text{s}}^{2}}{f}_{\rm d}^{2}\left(m\right){t}_{1},\text{j}\frac{4\pi {K}_{\text{r}}}{\text{c}}\cdot \frac{{\lambda }^{2}{h}_{\text{ref}}^{2}}{8{V}_{\text{s}}^{2}}{f}_{\rm d}^{2}\left(m\right){t}_{2},\cdots ,\text{j}\frac{4\pi {K}_{\text{r}}}{\text{c}}\cdot \frac{{\lambda }^{2}{h}_{\text{ref}}^{2}}{8{V}_{\text{s}}^{2}}{f}_{\rm d}^{2}\left(m\right){t}_{N}\right]\end{aligned}\right\}$

(12)

其中， ${\left[ \cdot \right]^{\text{T}}}$ 代表转置运算。对接收回波以及重跟踪回波模型对目标高程参数 $\tau$ 建立数学模型后，本文通过模型驱动的机器学习方式，面向式(10)求解高程参数 $\tau$ 问题，提出一种参数化重跟踪贝叶斯算法PR-Bayes。

3. 高度计回波统计建模

当式(10)所示重跟踪参数模型建立完后，为避免复杂地形高程参数估计出现过拟合现象，需对模型中各参数选取合适的先验统计建模，推导参数后验概率密度函数解析式并获得闭合解，实现高精度可信估计。

3.1 高程参数先验建模

相比于正则先验的灵活性不高，概率先验可根据目标特性自适应建立随机概率模型，充分考虑目标特性的“不确定性”。为了获得高精度高程信息，需对式(10)中的高程参数 $\tau$ 选取合适的先验分布进行建模。因 $\tau$ 代表天底点1维绝对高度数据，其对于空域的接收回波具有明显的稀疏特性，因此本文对高程参数 $\tau$ 的先验分布约束为稀疏特性良好的拉普拉斯分布^[17,18]，其表达式为

$f\left( {\tau |\lambda } \right) = {\text{Laplace}}\left( {\tau |\lambda } \right) \propto {\text{exp}}\left( { - {{\left| \tau \right|} \mathord{\left/ {\vphantom {{\left| \tau \right|} \lambda }} \right. } \lambda }} \right)$

(13)

其中， ${\text{Laplace}}(\cdot)$ 为拉普拉斯分布， $\tau$ 与 $\lambda$ 均为正值， $\lambda$ 为此拉普拉斯分布的尺度参数。

3.2 回波数据似然函数建模

完成先验建模后，根据贝叶斯推论，还需对式(10)进行似然函数建模。为完成该建模，需先对噪声选取合适的概率分布，考虑到噪声功率的随机性，因高斯噪声方差 ${\sigma ^2}$ 具有随机性，故对其建模为逆伽马分布^[16]

$\begin{split} f\left( {{\sigma ^2}|a,b} \right) =& \prod\limits_{k = 1}^K {\prod\limits_{n = 1}^N {\mathcal{I}\mathcal{G}\left( {\sigma _{kn}^2|a,b} \right)} } \\ =& \prod\limits_{k = 1}^K {\prod\limits_{n = 1}^N {\left[ {\left( {{{{b^a}} / {\varGamma \left( a \right)}}} \right)\sigma _{kn}^{ - 2\left( {a{\text{ + }}1} \right)} \exp \left( { - {b/ {\sigma _{kn}^2}}} \right)} \right]} } \end{split}$

(14)

其中， $\mathcal{I}\mathcal{G}(\cdot )$ 代表逆伽马分布， $a$ 和 $b$ 分别为逆伽马分布的形状参数和尺度参数。此时便可对回波数据 ${\boldsymbol{Y}}$ 进行似然建模。通常，距离压缩后不同距离单元之间独立同分布，且服从圆对称复高斯分布^[16]。由式(10)可得回波数据似然函数服从建模为

$\begin{split} f({\boldsymbol{Y}}\mid \tau ,{\sigma ^2}) =& \prod\limits_{k = 1}^K {\prod\limits_{n = 1}^N {\mathcal{C}\mathcal{N}} \left( {{{\boldsymbol{Y}}_{kn}}\mid \tau ,\sigma _{kn}^2} \right)} \\ =& \prod\limits_{k = 1}^K \prod\limits_{n = 1}^N {{\left( {2\pi \sigma _{kn}^2} \right)}^{ - \frac{1}{2}}}\\ & \cdot \exp \left\{ { - {{{{\left| {{{\left[ {{\boldsymbol{Y}} - {\boldsymbol{AX}}} \right]}_{kn}}} \right|}^2}} {\left( {2\sigma _{kn}^2} \right)}}} \right\} \end{split}$

(15)

其中， ${\boldsymbol{X}}$ = ${{\boldsymbol{I}}_{\boldsymbol{M}}}{\bf{BM}}\left( \tau \right) \odot {{\boldsymbol{F}}_{\rm{d}}}$ 。

3.3 分层贝叶斯建模

如式(13)、式(15)所示的先验与似然函数模型，二者在贝叶斯推论下非共轭故无法获得闭合解，针对此问题，需对式(13)的先验进行分层建模以实现二者共轭。根据文献[16,18]可知，可将拉普拉斯分布分两层建模为

$f\left( {\tau |\alpha } \right) = \mathcal{C}\mathcal{N}\left( {\tau |0,\alpha } \right) = {\left( {2\pi \alpha } \right)^{ - \frac{1}{2}}}\exp \left( { - {{{\tau ^2}} \mathord{\left/ {\vphantom {{{\tau ^2}} {2{\alpha ^2}}}} \right. } {2{\alpha ^2}}}} \right)$

(16)

$\begin{split} f\left( {\alpha |\lambda } \right) = & {\text{Gamma}}\left( {\alpha |\eta ,\lambda } \right) = \left( {{{{\lambda ^\eta }} \mathord{\left/ {\vphantom {{{\lambda ^\eta }} {\varGamma \left( \eta \right)}}} \right. } {\varGamma \left( \eta \right)}}} \right){\alpha ^{ - \eta - 1}}\\ & \cdot {\text{exp}}\left( { - \lambda \alpha } \right) \end{split}$

(17)

其中，将第1层建模为关于高程参数 $\tau$ 的复高斯分布， $\alpha$ 为超参数，由于高斯函数固有的平滑特性，式(13)不属于稀疏分布范畴。为促进先验的稀疏特性并简化后验分布的推导，本文选取与高斯分布成对共轭的伽马分布作为超参数 $\alpha$ 的先验，并在分层模型的第2层将超参数 $\alpha$ 建模成如式(17)的形式，其中 $\lambda$ 为尺度参数。此时可根据式(16)和式(17)计算边缘分布，如当 $\eta$ 的值很小时， $f\left( {\tau |\lambda } \right) = \displaystyle\int {f\left( {\tau |\alpha } \right)} \cdot f\left( {\alpha |\lambda } \right){\text{d}}\alpha$ ，该分布服从拉普拉斯分布^[16,18]。上述贝叶斯分层中参数与超参数的关系可用有向无环图DAG进行表示，如图3所示。

图 3 贝叶斯分层模型DAG

下载: 全尺寸图片幻灯片

综上所述，提出的贝叶斯模型包括参数 $\tau ,{\sigma ^2}$ 以及超参数 $\lambda$ ，可将 $\tau$ 和其他参数的先验联合表示其后验概率密度，故该贝叶斯模型的后验分布可以从式(18)层次结构计算得出

$\begin{split} & f(\tau ,\alpha ,\lambda ,{\sigma ^2}|{\boldsymbol{Y}}) \propto\\ & \frac{{f({\boldsymbol{Y}}|\tau ,{\sigma ^2})f(\tau |\alpha )f(\alpha |\lambda )f({\sigma ^2})}}{{\displaystyle\int {\displaystyle\int {\displaystyle\int {\displaystyle\int {f({\boldsymbol{Y}}|\tau ,{\sigma ^2})f(\tau |\alpha )f(\alpha |\lambda )f({\sigma ^2}){\rm{d}}\tau {\rm{d}}\alpha {\rm{d}}\lambda {\rm{d}}{\sigma ^2}} } } } }} \end{split}$

(18)

其中，等号右侧分子部分可由式(14)—式(17)计算得出，而分母则为对边缘概率密度函数进行多重积分获得。然而该分母的多重积分在实际中无法获得闭合解，通常考虑变分贝叶斯(Variational Bayes, VB)或者MCMC统计采样两种解决方案，VB算法通常用于计算典型模型，复杂模型计算困难，问题依赖性强^[19,20]。而MCMC算法的普适性更强，故本文采用后者对该类型后验分布进行解析计算获得高精度参数估计。

4. 高程参数贝叶斯统计学习

针对复杂后验分布计算求解，本文提出的PR-Bayes算法引入新型HMC采样，与常规MCMC算法依赖马尔科夫链进行递归采样不同，HMC算法基于物理学哈密顿动态方程，计算更为简便，通过状态转移的方式对参数后验概率密度函数进行参数更新以获得解析式闭合解。

4.1 参数后验分布计算

由式(15)、式(16)和式(17)，对高程参数 $\tau$ 进行边缘后验分布推导，如式(19)所示

$\begin{split} f\left(\tau |{\boldsymbol{Y}},\alpha ,\lambda ,{\sigma }^{2}\right)\propto & f\left({\boldsymbol{Y}}|\lambda ,{\sigma }^{2}\right)\cdot f\left(\tau |\alpha \right)\cdot f\left(\alpha |\lambda \right)\\ \propto& \displaystyle \prod _{k\text{=1}}^{K}{\displaystyle \prod _{n=1}^{N}{\left(\alpha {\sigma }_{kn}^{2}\right)}^{-\frac{\text{1}}{2}}}{\alpha }^{c-1}\\ & \cdot \mathrm{exp}\left\{-{\left|{\left({\boldsymbol{Y}}-{\boldsymbol{AX}}\right)}_{kn}\right|}^{2}/2{\sigma }_{kn}^{2}\right. \\ & \left. -\left({\tau }^{\text{2}}/\left(2{\alpha }^{\text{2}}\right)+\lambda \alpha \right)\right\}\\[-10pt] \end{split}$

(19)

在获得主要目标参数 $\tau$ 的后验分布后，可由式(14)和式(15)对噪声参数 ${\sigma ^2}$ 进行后验推导，噪声参数 ${\sigma ^2}$ 服从

$\begin{split} f\left({\sigma }^{2}|{\boldsymbol{Y}},\tau ,\lambda \right)\propto& f\left({\boldsymbol{Y}}|\tau ,{\sigma }^{2}\right)\cdot f\left({\sigma }^{2}\right)\\ \propto &\displaystyle \prod _{k\text{=1}}^{K}{\displaystyle \prod _{n=1}^{N}\mathcal{I}\mathcal{G}}\left({\sigma }_{kn}^{2}|a+\frac{\text{1}}{2},\right.\\ & b+{\left|{\left({\boldsymbol{Y}}-{\boldsymbol{AX}}\right)}_{kn}\right|}^{2}/2) \end{split}$

(20)

可看出噪声参数 ${\sigma ^2}$ 服从常规逆伽马分布。在完成上述后验推导后，由式(16)、式(17)对超参数 $\alpha$ 进行推导

$\begin{split} f\left(\alpha |\tau ,\lambda \right)\propto &f\left(\tau |\alpha \right)\cdot f\left(\alpha \right)\\ \propto & {\alpha }^{\eta -\frac{3}{2}}\cdot \mathrm{exp}\left\{-\left(\lambda \alpha +{\tau }^{2}/\left(2{\alpha }^{2}\right)\right)\right\} \end{split}$

(21)

高程参数 $\tau$ 及其超参数 $\alpha$ 的边缘后验分布并不满足于常规分布函数，对于此类非常规分布，便可使用上文所提的HMC算法对其进行求解^[7,10,20]，接下来详细叙述HMC算法求解目标参数的计算过程。

4.2 哈密顿蒙特卡洛抽样算法

经典MCMC算法构造马尔科夫链进行迭代收敛实现对复杂分布的采样，然而高维特征下该方式易出现随机游走影响模型参数的采样性能。HMC算法通过哈密顿方程对该迭代的随机性进行约束，从而大大提高收敛效率。哈密顿动力学通过一个物体在某个时间点 $t$ 的位置 $x$ 及其动量 $p$ 来描述运动，即物体具有一定势能 $U\left( x \right)$ 和动能 $K\left( p \right)$ ，通过二者的标量和来定义物体的总能量 $H\left( {x,p} \right) = U\left( x \right) + K\left( p \right)$ 。此时可基于能量守恒定理，将此物体动能与势能实现相互转换并构造哈密顿动态分程^[7]，实现如式(22)分布采样

$f\left( x \right) \propto \exp \left[ { - U\left( x \right)} \right]$

(22)

HMC算法根据式(22)利用梯度信息进行蛙跳(leap frog)步骤更新变量。蛙跳法根据一定次序逐渐更新物体的位置和动量，先计算极小一段时间 $\delta {\text{/2}}$ 内的动量变化，再计算经过 $\delta$ 时间运动后的物体位置，最后再计算完成另外 $\delta {\text{/2}}$ 时间间隔的动量计算。这样便可计算出相同时间点上的 $p$ 和 $x$ ^[7]。这种方法为

$\left. \begin{split} & {p_i}\left( {t + \frac{\delta }{2}} \right) = {p_i}\left( t \right) - \frac{\delta }{2}\frac{{\partial U}}{{\partial {x_i}\left( t \right)}}\\ & {x_i}\left( {t + \delta } \right) = {x_i}\left( t \right) + \delta \frac{{\partial K}}{{\partial {p_i}\left( {t + {\delta / 2}} \right)}}\\ & {p_i}\left( {t + \delta } \right) = {p_i}\left( {t + \frac{\delta }{2}} \right) - \frac{\delta }{2}\frac{{\partial U}}{{\partial {x_i}\left( {t + \delta } \right)}} \end{split}\right\}$

(23)

针对本文高程回波模型问题，定义势能函数为

$U\left( \tau \right) = \sum\limits_{k = 1}^K {\sum\limits_{n = 1}^N {\left( {\frac{1}{{2\sigma _{kn}^2}} {{\left| {{{\left( {{\boldsymbol{Y}} - {\boldsymbol{AX}}} \right)}_{kn}}} \right|}^2} + \left( {{{{\tau ^{\text{2}}}}/ {\left( {2{\alpha ^{\text{2}}}} \right)}} + \lambda \alpha } \right)} \right)} }$

(24)

根据HMC采样原理，利用上述蛙跳步骤式(24)对变量 $\tau$ 进行更新，通过不断更新迭代直至达到收敛，通过采集变量 $\tau$ 的收敛样本并求其期望，便可获得目标的高度参数 $\tau$ 的估计值 ${\hat \tau _{{\text{HMC}}}}$ ，为使HMC采样正确收敛至目标函数，采样过程中需使用“接收-拒绝”步骤进行样本筛选。接收概率表示为^[17,18]

$\min \left\{ {1,\exp \left( {H\left( {\tau ,p} \right) - H\left( {\tau *,p*} \right)} \right)} \right\}$

(25)

其中， $\tau *,p*$ 代表收敛样本，综上。根据回波模型后验采样求解算法，可总结出PR-Bayes算法参数估计流程。

4.3 PR-Bayes算法流程

SARA参数估计流程如图4所示。在粗跟踪基础上对回波数据进行dechirp等处理得到多视回波。参数初始化完成后，输入经处理得到的多视结果和粗跟踪高度值，基于布朗模型采用PR-Bayes算法进行重跟踪参数建模，并根据式(19)—式(21)分别对参数 $\tau$ , ${\sigma ^2}$ , $\alpha$ 应用HMC算法进行采样。当获得收敛样本后，可将此样本用于参数和变量的估计和计算。此后，对采样样本计算期望，以降低误差影响，并最终得到目标参数 $\tau$ 的采样估计值，获得相应目标点的高度重跟踪结果，实现对目标地形高程参数求解。

图 4 本文所提PR-Bayes算法运算流程

下载: 全尺寸图片幻灯片

5. 实验验证

上面已对所提PR-Bayes算法原理进行阐述分析，本节通过数值精确度对该算法进行验证并设计对比实验，验证本文所提PR-Bayes算法的有效性及优越性。本节共包括点目标仿真、数字高程模型(Digital Elevation Model, DEM)半实物模拟以及实测数据验证3部分实验。

5.1 点目标仿真模拟实验

为了验证高维情况下，引入先验可改善过拟合现象，提高参数估计精度，本部分通过设立地面点目标，仿真SARA回波的形式，分别使用PR-Bayes和LS算法对仿真回波进行重跟踪处理，并通过得到的高程参数进行数值验证。在进行验证前，需使得布朗模型和回波有着良好的匹配度，同时为了更好地对算法性能进行数值验证，还需考虑实际探测过程中的噪声影响。本部分仿真设立无噪声、信噪比分别为0 dB和–10 dB 3种噪声模式。点目标仿真实验中仿真参数和实验结果分别如表1和图5所示。

表 1 点目标仿真雷达参数

雷达参数	数值	雷达参数	数值	雷达参数	数值
雷达载频	9.6 GHz	采样频率	300 MHz	波长	3.13 cm
脉冲宽度	4 μs	天线孔径	0.4 m	飞行高度	1000 m
信号带宽	200 MHz	脉冲重复频率	2000 Hz	载机速度	60 m/s

下载: 导出CSV

| 显示表格

图 5 点目标估计结果

下载: 全尺寸图片幻灯片

图5将得到的估计结果通过布朗模型表示并与仿真回波进行拟合度对比，可根据拟合度初步判断算法优劣。其中，可看到无噪声环境下仿真回波与布朗模型具有良好的匹配度。图5(a)为LS算法估计结果，图5(b)是PR-Bayes算法结果。从中可以看出，相较于LS算法，PR-Bayes算法有着更好的拟合度，且随着噪声的增加，LS算法拟合度逐渐下降，而PR-Bayes算法仍有较高的拟合度。表2总结了不同信噪比下得到的真值和相应算法参数估计结果，并采用式(26)进行精度提升计算，其中 ${\tau _0}$ 为真值。从结果可得到，噪声环境下，LS估计结果逐渐失真，而本文所提PR-Bayes算法仍有较高的估计精度

表 2 点目标仿真实验结果

噪声(dB)	真值	LS估计值	PR-Bayes估计值	精度提升(%)
无噪声	50.40	50.04	50.04	0.55
0	50.40	51.02	50.08	50.25
-10	50.40	51.49	49.79	44.23

下载: 导出CSV

| 显示表格

$\Delta \tau = \frac{{\left| {\left| {{\tau _{{\text{HMC}}}} - {\tau _0}} \right| - \left| {{{\hat \tau }_{{\text{LS}}}} - {\tau _0}} \right|} \right|}}{{\left| {{{\hat \tau }_{{\text{LS}}}} - {\tau _0}} \right|}} \times 100\%$

(26)

为进一步验证算法稳定性，本文设置蒙特卡洛实验进行大量数据测试，随机生成1000组回波进行参数估计，将 $\tau$ 的估计偏差作为评估指标，以此衡量两种算法的估计精度及稳定性，理论上大量随机实验结果服从高斯分布，实验结果如图6所示。两种算法估计偏差概率大致服从高斯分布，因PR-Bayes算法估计偏差集中在0附近且几乎无波动，而LS算法估计偏差波动较大，故PR-Bayes算法估计偏差明显优于LS算法。

图 6 蒙特卡洛模拟实验结果

下载: 全尺寸图片幻灯片

5.2 半实物模拟仿真实验

进行点目标仿真数据验证算法精度提升后，进一步选取DEM半实物仿真数据验证参数估计性能。选取部分DEM作为真值读取连续变化的高度信息，并基于此按照同心圆法获取仿真模拟回波数据，进而开展实验验证。针对不同场景该算法的适用性，选取两种地形进行模拟仿真实验，分别为具有缓变地形特征的平原地形和具有突变地形特征的山区地形。

实验结果如图7所示，其中上下两行分别为平原地形和山区地形的实验结果。图7(a)、图7(c)分别为两种地形的DEM，矩形框内区域为高度估计区域。图7(b)、图7(d)分别为重跟踪的结果对比曲线图，其中分别包含了DEM实际高度数据和LS, PR-Bayes两种算法的估计结果，可看出两种算法地形跟踪趋势均符合DEM真实高度变化。为定量分析算法精度，计算两种算法估计结果相对于真值的标准偏差(STandard Deviation, STD)值。因DEM数据高度真值已知，为便于计算对估计结果进行插值，STD定义为^[3]

图 7 DEM半实物模拟仿真实验结果

下载: 全尺寸图片幻灯片

${\text{STD = }}\sqrt {\frac{{\text{1}}}{I}\sum\limits_i^I {{{\left( {\tau - {{\hat \tau }_i}} \right)}^2}} }$

(27)

其中， $I$ 为估计点总数， $\tau$ 为该点真实高度值， $\hat \tau$ 为算法估计值，计算结果如表3所示。由表3可明确得出，PR-Bayes算法的STD值明显小于LS算法，同时PR-Bayes估计算法对于缓变地形和突变地形均具有良好的适用性，且有效改善了LS算法出现过拟合等问题，可实现对复杂地形的高度跟踪。

表 3 半实物仿真数据重跟踪结果定量分析

指标(m)	LS	PR-Bayes
平原STD	2.7852	2.0014
山区STD	1.7811	1.5231

下载: 导出CSV

| 显示表格

5.3 机载实测数据实验

接下来，通过机载合成孔径雷达高度计实测数据验证算法实际性能。本文所用实测数据由中国工程物理研究院电子工程研究所开发的机载合成孔径雷达高度计系统所得，系统参数如表4所示。由于无法事先获取地面场景高程真值，这里结合实际情况设定高程参考值，用以衡量测高精度。

表 4 机载实测数据参数

雷达参数	数值	雷达参数	数值	雷达参数	数值
雷达载频	9.6 GHz	采样频率	125 MHz	波长	3.13 cm
脉冲宽度	5 μs	天线孔径	0.4 m	飞行高度	2600 m
信号带宽	100 MHz	脉冲重复频率	2000 Hz	载机速度	60 m/s

下载: 导出CSV

| 显示表格

相关实验结果如图8所示，其中图8(a)为实测数据DDM 2维象，选取其中部分数据作为对象，分别使用两种算法，估计效果如图8(b)所示。从跟踪趋势可以看出，LS算法与PR-Bayes算法高度基本符合高程变化趋势，而本文所提算法估计结果，与实测地形高度变化趋势更为贴合。同时LS估计结果在个别点出现坏值，PR-Bayes算法估计结果更为平滑，有效改善了传统LS算法的过拟合问题。因无法实际获得探测地形高度的真实值，本文选取10个连续回波的高度估计平均值作为参考值^[21]，分别与两种算法跟踪结果进行STD值计算，结果如表5所示，其中PR-Bayes算法的STD值优于LS算法，实际证明了本文雷达高度计的参数化贝叶斯估计性能的优越性。

图 8 实测数据实验结果

下载: 全尺寸图片幻灯片

表 5 实测数据重跟踪结果定量分析

指标(m)	LS	PR-Bayes
STD	14.74	12.22

下载: 导出CSV

| 显示表格

6. 结论

本文基于机载SARA测高原理，针对传统LS重跟踪算法估计精度存在上限以及易出现过拟合等问题，提出了具有高精度的PR-Bayes可信算法。对布朗模型做出相关改进并作为重跟踪参数模型，基于贝叶斯学习框架建立回波和布朗模型高程参数之间的数学模型，设立合理有效的观测目标先验信息，改善了高维特征下过拟合问题，有效抑制了噪声影响，提高了测高精度。最后，通过点目标仿真实验、半实物模拟仿真实验和实测数据实验对所提PR-Bayes进行评估，验证了本文所提重跟踪算法优越性，可为SARA高精度高度测量提供方法论基础，为复杂地形跟踪提供技术基础。

图 1 合成孔径雷达高度计(SARA)工作原理示意图

下载: 全尺寸图片幻灯片

图 2 2维延时多普勒像(DDM)及点目标响应函数

下载: 全尺寸图片幻灯片

图 3 贝叶斯分层模型DAG

下载: 全尺寸图片幻灯片

图 4 本文所提PR-Bayes算法运算流程

下载: 全尺寸图片幻灯片

图 5 点目标估计结果

下载: 全尺寸图片幻灯片

图 6 蒙特卡洛模拟实验结果

下载: 全尺寸图片幻灯片

图 7 DEM半实物模拟仿真实验结果

下载: 全尺寸图片幻灯片

图 8 实测数据实验结果

下载: 全尺寸图片幻灯片

表 1 点目标仿真雷达参数

雷达参数	数值	雷达参数	数值	雷达参数	数值
雷达载频	9.6 GHz	采样频率	300 MHz	波长	3.13 cm
脉冲宽度	4 μs	天线孔径	0.4 m	飞行高度	1000 m
信号带宽	200 MHz	脉冲重复频率	2000 Hz	载机速度	60 m/s

下载: 导出CSV

表 2 点目标仿真实验结果

噪声(dB)	真值	LS估计值	PR-Bayes估计值	精度提升(%)
无噪声	50.40	50.04	50.04	0.55
0	50.40	51.02	50.08	50.25
-10	50.40	51.49	49.79	44.23

下载: 导出CSV

表 3 半实物仿真数据重跟踪结果定量分析

指标(m)	LS	PR-Bayes
平原STD	2.7852	2.0014
山区STD	1.7811	1.5231

下载: 导出CSV

表 4 机载实测数据参数

雷达参数	数值	雷达参数	数值	雷达参数	数值
雷达载频	9.6 GHz	采样频率	125 MHz	波长	3.13 cm
脉冲宽度	5 μs	天线孔径	0.4 m	飞行高度	2600 m
信号带宽	100 MHz	脉冲重复频率	2000 Hz	载机速度	60 m/s

下载: 导出CSV

表 5 实测数据重跟踪结果定量分析

指标(m)	LS	PR-Bayes
STD	14.74	12.22

下载: 导出CSV

参考文献(21)

[1]	刘鹏, 许可, 王磊, 等. 合成孔径雷达高度计与传统高度计精度比对分析与机载试验验证[J]. 电子与信息学报, 2016, 38(10): 2495–2501. doi: 10.11999/JEIT151354 LIU Peng, XU Ke, WANG Lei, et al. Precision comparison and airborne experiment validation between SAR altimeter and conventional altimeter[J]. Journal of Electronics &Information Technology, 2016, 38(10): 2495–2501. doi: 10.11999/JEIT151354
[2]	QUILFEN Y and CHAPRON B. On denoising satellite altimeter measurements for high-resolution geophysical signal analysis[J]. Advances in Space Research, 2021, 68(2): 875–891. doi: 10.1016/j.asr.2020.01.005
[3]	SHI Lingwei, XU Ke, LIU Peng, et al. Height precision of SAR altimeter and conventional radar altimeter based on flight experimental data[J]. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2016, 9(6): 2676–2686. doi: 10.1109/JSTARS.2016.2550030
[4]	杨双宝, 翟振和, 许可, 等. 合成孔径雷达高度计数据处理方法[J]. 遥感技术与应用, 2017, 32(6): 1083–1092. doi: 10.11873/j.issn.1004-0323.2017.6.1083 YANG Shuangbao, ZHAI Zhenhe, XU Ke, et al. The ground process segment of SAR altimeter[J]. Remote Sensing Technology and Application, 2017, 32(6): 1083–1092. doi: 10.11873/j.issn.1004-0323.2017.6.1083
[5]	HALIMI A, MAILHES C, TOURNERET J Y, et al. A semi-analytical model for delay/Doppler altimetry and its estimation algorithm[J]. IEEE Transactions on Geoscience and Remote Sensing, 2014, 52(7): 4248–4258. doi: 10.1109/TGRS.2013.2280595
[6]	NEMETH C and FEARNHEAD P. Stochastic gradient Markov chain Monte Carlo[J]. Journal of the American Statistical Association, 2021, 116(533): 433–450. doi: 10.1080/01621459.2020.1847120
[7]	CHAARI L, TOURNERET J Y, CHAUX C, et al. A Hamiltonian Monte Carlo method for non-smooth energy sampling[J]. IEEE Transactions on Signal Processing, 2016, 64(21): 5585–5594. doi: 10.1109/TSP.2016.2585120
[8]	VAJARGAH K F, BENIS S G, and GOLSHAN H M. Detection of the quality of vital signals by the Monte Carlo Markov Chain (MCMC) method and noise deleting[J]. Health Information Science and Systems, 2021, 9(1): 26. doi: 10.1007/s13755-021-00157-5
[9]	BRESSON G, CHATURVEDI A, RAHMAN M A, et al. Seemingly unrelated regression with measurement error: Estimation via Markov Chain Monte Carlo and mean field variational Bayes approximation[J]. The International Journal of Biostatistics, 2021, 17(1): 75–97. doi: 10.1515/ijb-2019-0120
[10]	BAISTHAKUR S and CHAKRABORTY A. Experimental verification for load rating of steel truss bridge using an improved Hamiltonian Monte Carlo-based Bayesian model updating[J]. Journal of Civil Structural Health Monitoring, 2021, 11(4): 1093–1112. doi: 10.1007/s13349-021-00495-8
[11]	BROWN G. The average impulse response of a rough surface and its applications[J]. IEEE Transactions on Antennas and Propagation, 1977, 25(1): 67–74. doi: 10.1109/TAP.1977.1141536
[12]	BROWN G S. A useful approximation for the flat surface impulse response[J]. IEEE Transactions on Antennas and Propagation, 1989, 37(6): 764–767. doi: 10.1109/8.29363
[13]	何华锋, 戴嘉琪, 贺友. 弹载雷达导引头测高回波模型仿真研究[J]. 电光与控制, 2018, 25(11): 11–14,29. doi: 10.3969/j.issn.1671-637X.2018.11.002 HE Huafeng, DAI Jiaqi, and HE You. Simulation of altimeter echo model for missile-borne radar seeker[J]. Electronics Optics &Control, 2018, 25(11): 11–14,29. doi: 10.3969/j.issn.1671-637X.2018.11.002
[14]	EGIDO A and SMITH W H F. Fully focused SAR altimetry: Theory and applications[J]. IEEE Transactions on Geoscience and Remote Sensing, 2017, 55(1): 392–406. doi: 10.1109/TGRS.2016.2607122
[15]	WARNER B D. Radar terrain return from the spherical earth at near vertical incidence[R]. AD260172, 1961.
[16]	杨磊, 夏亚波, 毛欣瑶, 等. 基于分层贝叶斯Lasso的稀疏ISAR成像算法[J]. 电子与信息学报, 2021, 43(3): 623–631. doi: 10.11999/JEIT200292 YANG Lei, XIA Yabo, MAO Xinyao, et al. Sparse ISAR imaging algorithm based on Bayesian-lasso[J]. Journal of Electronics &Information Technology, 2021, 43(3): 623–631. doi: 10.11999/JEIT200292
[17]	YANG Lei, ZHAO Lifan, BI Guoan, et al. SAR ground moving target imaging algorithm based on parametric and dynamic sparse Bayesian learning[J]. IEEE Transactions on Geoscience and Remote Sensing, 2016, 54(4): 2254–2267. doi: 10.1109/TGRS.2015.2498158
[18]	杨磊, 李埔丞, 李慧娟, 等. 稳健高效通用SAR图像稀疏特征增强算法[J]. 电子与信息学报, 2019, 41(12): 2826–2835. doi: 10.11999/JEIT190173 YANG Lei, LI Pucheng, LI Huijuan, et al. Robust and efficient sparse-feature enhancement for generalized SAR imagery[J]. Journal of Electronics &Information Technology, 2019, 41(12): 2826–2835. doi: 10.11999/JEIT190173
[19]	PROTIN F, JULES M, NGUYEN D T, et al. Unified modelling of epidemics by coupled dynamics via Monte-Carlo Markov Chain algorithms[J]. arXiv: 2106.13463, 2021.
[20]	YAO Yu and STEPHAN K E. Markov chain Monte Carlo methods for hierarchical clustering of dynamic causal models[J]. Human Brain Mapping, 2021, 42(10): 2973–2989. doi: 10.1002/hbm.25431
[21]	HALIMI A, MAILHES C, TOURNERET J Y, et al. Bayesian estimation of smooth altimetric parameters: Application to conventional and delay/Doppler altimetry[J]. IEEE Transactions on Geoscience and Remote Sensing, 2016, 54(4): 2207–2219. doi: 10.1109/TGRS.2015.2497583