一种线上贝叶斯压缩水下成像方法及系统、设备和介质

1.本发明属于声呐成像技术领域，特别是涉及一种线上贝叶斯压缩水下成像方法及系统、设备和介质。


背景技术：

2.水下成像系统能够通过声回波重建目标强度。声成像作为一种重构水下环境和目标的有效方法，被广泛应用于海洋勘探、测绘和环境监测。在考虑远场声传播和后向散射的情况下，波束形成是开发声成像系统最常用的方法，其模型和伴随算子可以通过快速傅里叶变换(fft)和快速傅里叶反变换(ifft)有效实现。
3.利用压缩感知方法实现波束形成能够有效提高成像质量。常用的求解方法包括正交匹配追踪(omp)、快速迭代收缩阈值算法(fista)、近似消息传递(amp)及稀疏贝叶斯学习(sbl)等。上述方法在处理多快拍波束形成问题时，一般采用“成块”处理的方式，各快拍之间采用相互独立假设，例如多快拍sbl(msbl)和块近端梯度法(bpg)。其明显缺陷在于，难以构建各快拍之间的结构关系。
4.另一种方式是建立动态系统来实现多快拍的求解问题。具体方法是，引入不随时间变换的状态转移矩阵构建各快拍之间满足的状态转移方程，然后通过基于kalman滤波的稀疏编码方法进行高精度重构。然而，利用这种方法进行水下图像重构的存在如下两个困难：一是不能通过物理方法显式地给出状态转移矩阵的形式，另一方面，基于稀疏编码的求解方法大多数是以迭代机制实现，在满足某些收敛条件时才能判定输出，不但计算量非常大，而且迭代次数具有不确定性，给在线成像系统的实现造成了困难。


技术实现要素：

5.本发明为了解决现有技术中的问题，提出一种线上贝叶斯压缩水下成像方法及系统、设备和介质。本发明通过建立基于波束形成的测量方程和状态转移方程，给出非时变的先验稀疏假设，该假设保证了待重构的源向量的稀疏一致性；通过期望最大(em)算法推导得到超参数的更新公式，包括稀疏超参数γ，噪声不确定度β和状态转移矩阵d；为降低计算量，采用锯齿滞后机制对参数进行近似估计，结合kalman滤波和rts平滑得到周期性的参数迭代计算公式。
6.本发明是通过以下技术方案实现的，本发明提出一种线上贝叶斯压缩水下成像方法，所述方法包括：
7.接收并采样存储水下目标的声回波，根据声呐系统的测量范围划分数据帧，将测量数据划分为长度为δ的数据块，并将待重构的状态更新数据划分为长度为δ
′
的数据块；按照规定模式完成超参数初始化，所述超参数初始化包含两种初始化方式，即，确定数值初始化和最优初始化；
8.将测量数据划分的数据块和感知矩阵输入线上贝叶斯压缩水下成像模型，将最新存储的长度为δ的测量数据块输入kalman滤波-rts平滑模块，通过kalman滤波器得到前向
递归的估计均值和协方差矩阵，并将结果馈入rts平滑器得到后向递归估计均值和协方差矩阵；
9.根据声呐系统参数建立测量方程和状态转移方程，通过em算法完成超参数估计；在实现时，启动锯齿滞后周期性更新机制，在第l个更新周期，根据kalman滤波-rts平滑模块输出的源向量状态统计信息，更新计算超参数；将更新得到的超参数和长度为δ
′
的状态更新数据块重新输入kalman滤波-rts平滑模块，得到长度为δ
′
的状态估计均值和协方差矩阵，其中均值向量输出为本周期的源向量估计结果；更新时间窗和对应的长度为δ的测量数据块，输入kalman滤波-rts平滑模块启动第l 1个周期的状态估计。
10.本发明还提出一种线上贝叶斯压缩水下成像系统，所述系统包括：
11.测量数据接收模块，用于接收并采样存储水下目标的声回波，根据声呐系统的测量范围划分数据帧，将测量数据划分为长度为δ的数据块，并将待重构的状态更新数据划分为长度为δ
′
的数据块；按照规定模式完成超参数初始化，所述超参数初始化包含两种初始化方式，即，确定数值初始化和最优初始化；
12.kalman滤波-rts平滑模块，用于将测量数据划分的数据块和感知矩阵输入线上贝叶斯压缩水下成像模型，将最新存储的长度为δ的测量数据块输入kalman滤波-rts平滑模块，通过kalman滤波器得到前向递归的估计均值和协方差矩阵，并将结果馈入rts平滑器得到后向递归估计均值和协方差矩阵；
13.图像重构模块，用于根据声呐系统参数建立测量方程和状态转移方程，通过em算法完成超参数估计；在实现时，启动锯齿滞后周期性更新机制，在第l个更新周期，根据kalman滤波-rts平滑模块输出的源向量状态统计信息，更新计算超参数；将更新得到的超参数和长度为δ
′
的状态更新数据块重新输入kalman滤波-rts平滑模块，得到长度为δ
′
的状态估计均值和协方差矩阵，其中均值向量输出为本周期的源向量估计结果；更新时间窗和对应的长度为δ的测量数据块，输入kalman滤波-rts平滑模块启动第l 1个周期的状态估计。
14.本发明提出一种电子设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现所述线上贝叶斯压缩水下成像方法的步骤。
15.本发明提出一种计算机可读存储介质，用于存储计算机指令，所述计算机指令被处理器执行时实现所述线上贝叶斯压缩水下成像方法的步骤。
16.本发明的有益效果是：
17.1、本发明所提出贝叶斯压缩成像方法以贝叶斯压缩感知和kalman滤波为基础，具有高精度图像重构能力，获得更高的成像质量；
18.2、本发明所提出的成像方法具有非迭代特征，即不需要通过反复迭代实现参数估计收敛，因此同时具备低计算量和计算量的可预测性，保证了线上算法的实现稳定性；同时，可证明采用的锯齿滞后更新的超参数具有收敛性；
19.3、本发明采用最大期望算法可得到状态转移矩阵的估计结果，从而充分利用了源向量各快拍之间的结构关系，以进一步提高成像精度。
附图说明
20.图1是锯齿滞后更新机制的时序原理图；
21.图2(a)是稀疏超参数γ收敛曲线；
22.图2(b)是噪声不确定度β收敛曲线；
23.图2(c)是相关系数ρ收敛曲线；
24.图3(a)是重构rmse与测量快拍数的性能曲线；
25.图3(b)是重构srr与测量快拍数的性能曲线；
26.图3(c)是运行时间与测量快拍数的性能曲线；
27.图3(d)是重构rmse与阵元数量的性能曲线；
28.图3(e)是重构srr与阵元数量的性能曲线；
29.图3(f)是运行时间与阵元数量的性能曲线；
30.图3(g)是重构rmse与snr的性能曲线；
31.图3(h)是重构srr与snr的性能曲线；
32.图3(i)是运行时间与snr的性能曲线；
33.图3(j)是重构rmse与稀疏度的性能曲线；
34.图3(k)是重构srr与稀疏度的性能曲线；
35.图3(l)是运行时间与稀疏度的性能曲线；
36.图4(a)是online-ksbl方法的成像结果；
37.图4(b)是msbl方法的成像结果；
38.图4(c)是dsc方法的成像结果；
39.图4(d)是tmsbl方法的成像结果。
40.图3(a)
‑‑
图3(l)中，带有方块的曲线为确定数值初始化online-ksblδ
′
＝3；带有三角的曲线为最优初始化online-ksblδ
′
＝3；带有圆点的曲线为msblδ
′
＝3。
具体实施方式
41.下面将结合本发明实施例中的附图对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
42.如图1所示，本实施例提供的线上贝叶斯压缩水下成像方法，采用锯齿滞后更新机制，该机制依赖于两个延迟参数，即状态更新数据块δ
′
和测量数据块δ。从序列图中可见，必有δ
′
≤δ，即符合锯齿状延迟。源向量的统计信息以块的形式估计，数据长度为δ
′
，同时完成一次超参数更新。在第l更新阶段，按顺序更新源向量xk，k∈[k
l
1,k
l
δ
′
]，在图1中以深灰色进行标识。其中起始快拍号为k
l
＝(l-1)δ
′
。同时，在该阶段用于超参数估计的测量数据块yk′
，k
′
∈[k
l
1,k
l
δ 1]，为了简化符号表示，采用k
l
′
＝k
l
δ 1。另外，图中浅灰色表示的测量数据是可以丢弃的，在以后的处理中也不会再使用。所述超参数，包括稀疏超参数、噪声不确定度和转移矩阵均以锯齿滞后周期性更新机制递归估计。
[0043]
本发明提出一种线上贝叶斯压缩水下成像方法，所述方法包括：
[0044]
接收并采样存储水下目标的声回波，根据声呐系统的测量范围划分数据帧，将测量数据划分为长度为δ的数据块，并将待重构的状态更新数据划分为长度为δ
′
的数据块；按照规定模式完成超参数初始化，所述超参数初始化包含两种初始化方式，即，确定数值初
始化和最优初始化；
[0045]
将测量数据划分的数据块和感知矩阵输入线上贝叶斯压缩水下成像模型，将最新存储的长度为δ的测量数据块输入kalman滤波-rts平滑模块，通过kalman滤波器得到前向递归的估计均值和协方差矩阵，并将结果馈入rts平滑器得到后向递归估计均值和协方差矩阵；
[0046]
根据声呐系统参数建立测量方程和状态转移方程，通过em算法完成超参数估计；在实现时，启动锯齿滞后周期性更新机制，在第l个更新周期，根据kalman滤波-rts平滑模块输出的源向量状态统计信息，更新计算超参数；将更新得到的超参数和长度为δ
′
的状态更新数据块重新输入kalman滤波-rts平滑模块，得到长度为δ
′
的状态估计均值和协方差矩阵，其中均值向量输出为本周期的源向量估计结果；更新时间窗和对应的长度为δ的测量数据块，输入kalman滤波-rts平滑模块启动第l 1个周期的状态估计。
[0047]
确定状态更新数据块长度δ
′
和测量数据块长度δ，并提供两种超参数初始化方案，即确定数值初始化和最优初始化，所述确定数值初始化规定所有超参数初始化为固定常数或常向量；所述最优初始化采用传统稀疏贝叶斯学习模型估计稀疏超参数向量γ0和噪声不确定度β0作为初值。
[0048]
利用kalman滤波器和rts平滑器得到源向量的统计信息。具体的：
[0049]
步骤1、在第l更新阶段，采用kalman滤波进行正向递归估计源向量的滤波均值和协方差p
t|t
，时间窗t∈[k
l
1,k
l
′
]。具体步骤如下：
[0050]
输入参数：测量数据块y
t
，t∈[k
l
1,k
l
′
]；感知矩阵a，转移矩阵d，稀疏协方差矩阵γ，噪声不确定度β，第l-1更新阶段的最终输出均值和协方差
[0051]
for t∈[k
l
1,k
l
′
]
[0052]
预测：
[0053]
p
t|t-1
＝dp
t-1|t-1dh
(i-dhd)γ
[0054]
增益：g
t
＝p
t|t-1ah
(ap
t|t-1ah
βi)-1
[0055]
滤波：
[0056]
p
t|t
＝p
t|t-1-gap
t|t-1
[0057]
end
[0058]
步骤2、采用rts平滑进行反向递归估计平滑均值和协方差t∈[k
l
1,k
l
′
]。具体步骤如下：
[0059]
输入参数：kalman滤波输出的状态统计信息和p
t|t
，t∈[k
l
1,k
l
′
]；转移矩阵d；
[0060]
初始化：
[0061]
for t＝k
l
′
,k
l
′‑
1,
…
,k
l
2
[0062]
增益：
[0063]
平滑：
[0064][0065]
if t≠k
l
′
[0066][0067]
end。
[0068]
以kalman滤波-rts平滑得到的源向量统计信息为基础，通过测量方程和状态转移方程，结合锯齿滞后周期性更新机制，更新估计超参数，具体为：
[0069]
步骤(1)：在已知感知矩阵和测量回波数据前提下，在时域上建立用于目标强度成像的线性方程，即测量方程；
[0070]
步骤(2)：对相邻时序的源向量建立状态转移方程，其中状态转移矩阵暂时是未知的，与水下成像中快拍结构关系的未知性一致；
[0071]
步骤(3)：根据期望最大算法，建立各超参数对应的q函数；
[0072]
步骤(4)：在上述q函数的基础上，给出更新公式和近似迭代公式。
[0073]
所述步骤(1)具体为：
[0074]
在时间序列t∈[1,2,
…
,t]建立如下测量方程：
[0075]yt
＝ax
t
n
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(1)
[0076]
其中，x
t
为t时刻源向量或状态变量，y
t
为t时刻源向量或状态变量的输出值，a∈cm×n为波束形成对应的感知矩阵，具有非时变特性，可以表示为：a＝[a(θ1),a(θ2),
…
,a(θn)]，其中预定义的波达角(doa)θn的转向矢量可以表示为：
[0077][0078]
其中，阵元数目为m，预定义doa角度采样数目为n，d为线性阵列的阵元间隔，λ为声信号的波长，n
t
满足复高斯分布，其中β为当前噪声不确定性，im为m
×
m的单位矩阵。
[0079]
所述步骤(2)具体为：
[0080]
在时间序列t∈[1,2,
…
,t]建立如下状态转移方程：
[0081]
x
t
＝dx
t-1
z
t
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(3)
[0082]
其中，z
t
满足复高斯分布状态转移矩阵d假设为对角阵，能够有效降低过拟合的发生，γ为时序源向量的协方差矩阵，其具有对角阵形式，即γ＝diag(γ)。
[0083]
所述步骤(3)具体为：
[0084]
由测量方程(1)得到似然概率密度：
[0085][0086]
假设源向量具有稀疏性，结合状态转移方程(3)得到转移概率密度：
[0087][0088]
因此得到联合概率密度：
[0089][0090]
其中，y
1:k
＝[y1,y2,
…
,yk]并且k∈[1,2,
…
,t]；当t＝1时,初始条件p(x1|x0；γ)＝p(x1；γ)；
[0091]
建立关于稀疏超参数γ的q函数：
[0092][0093]
为简化后续推导，期望表示符号简写为ek[
·
]，公式(7)中，第二项可以进一步化简为：
[0094][0095]
其中，p
1|k
和分别为x1相对于联合概率密度p(x
1:k
|y
1:k
；γ,β)的估计协方差矩阵和均值，其结果由rts平滑器估计得到；由此公式(7)可以写成：
[0096][0097]
其中，c
1|k
为t＝1时的中间矩阵，中间矩阵f
t|k
表示为：
[0098]ft|k
＝c
t|k
dc
t-1|kdh-c
t,t-1|kdh-dc
t-1,t|k
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(10)
[0099][0100][0101][0102]
建立关于噪声不确定性β的q函数：
[0103][0104]
假设状态转移矩阵d为对角阵，其对角元素d＝[d1,d2,
…
,dn]，满足|di|《1；i＝1,2,
…
,n；关于状态转移矩阵d的q函数写成：
[0105][0106]
为简化后续推导，假设γ-1
(i-dhd)的对角元素远小于1，即|di|2＜＜1，i＝1,2,
…
,n；于是q函数(15)可以重新写成：
[0107][0108]
其中，c
t|k
(i,i)表示矩阵c
t|k
的对角元素。
[0109]
所述步骤(4)具体为：
[0110]
锯齿波滞后更新机制，引入如下近似：
[0111][0112][0113][0114]
据此，在第l更新阶段，稀疏超参数γ的递归迭代公式为：
[0115][0116]
噪声不确定度β的递归迭代公式为：
[0117][0118]
状态转移矩阵d的对角阵元素di的递归迭代公式为：
[0119][0120]
本发明还提出一种线上贝叶斯压缩水下成像系统，所述系统包括：
[0121]
测量数据接收模块，用于接收并采样存储水下目标的声回波，根据声呐系统的测量范围划分数据帧，将测量数据划分为长度为δ的数据块，并将待重构的状态更新数据划分为长度为δ
′
的数据块；按照规定模式完成超参数初始化，所述超参数初始化包含两种初始化方式，即，确定数值初始化和最优初始化；
[0122]
kalman滤波-rts平滑模块，用于将测量数据划分的数据块和感知矩阵输入线上贝叶斯压缩水下成像模型，将最新存储的长度为δ的测量数据块输入kalman滤波-rts平滑模块，通过kalman滤波器得到前向递归的估计均值和协方差矩阵，并将结果馈入rts平滑器得到后向递归估计均值和协方差矩阵；
[0123]
图像重构模块，用于根据声呐系统参数建立测量方程和状态转移方程，通过em算法完成超参数估计；在实现时，启动锯齿滞后周期性更新机制，在第l个更新周期，根据kalman滤波-rts平滑模块输出的源向量状态统计信息，更新计算超参数；将更新得到的超参数和长度为δ
′
的状态更新数据块重新输入kalman滤波-rts平滑模块，得到长度为δ
′
的状态估计均值和协方差矩阵，其中均值向量输出为本周期的源向量估计结果；更新时间窗和对应的长度为δ的测量数据块，输入kalman滤波-rts平滑模块启动第l 1个周期的状态估计。
[0124]
在实验中，本发明所述的线上贝叶斯压缩成像方法被简称为online-ksbl。为了便于讨论，不失一般性地，将转移矩阵简化为d＝ρi，其中ρ为相关系数，并且ρ＝0.9表示高相关性的情况。图2(a)-图2(c)为各超参数的收敛性曲线。在δ
′
＝3,5，7三种条件下，在测量快拍数足够的情况下，能够保证各参数估计的收敛性。此外，不同δ
′
取值的情况下曲线偏差是非常明显的。δ
′
更大时，各参数将更频繁地更新和跟踪参数的真值，实现加速收敛。
[0125]
对于doa应用，在仿真中考虑角度和源振幅估计精度，使用以下指标进行了衡量：
ⅰ
.支撑域重构率(srr):
[0126][0127]
ⅱ
.相对均方误差(rmse):
[0128][0129]
其中，和x
t
分别表示第t个快拍源向量的估计值真实值。
[0130]
考虑到数值精度和噪声，(20)中的范数通过预定义的阈值η来近似判断，并在模拟中设置η＝10-2
。根据(20)和(21)的定义，srr是源向量支撑域的重构正确率，与角度精度有关。srr＝1表示源向量完美恢复。rmse表示源幅值的估计误差。此外，以程序的运行时间表征算法的复杂性。所有方法的程序重复100次，取平均值作为最终结果。
[0131]
仿真采用了两种初始化方案，即最优初始化和确定数值初始化。最优初始化为msbl采用前δ 1个测量数据估计超参数的初值；确定数值初始化使用如下设置：γ＝1,β＝std(y
t
)，d＝0.5
·
1，其中1表示全1向量，std(
·
)是标准差算子。
[0132]
以msbl作为比较方法，但是需要使用周期性更新处理。具体说，msbl按批量周期性处理，每次处理的数据块为δ
′
。因此，msbl的更新速率可以与online-ksbl同步。图3(a)
‑‑
图3(l)为online-ksbl和msbl在rmse、srr和运行时间三个方面的重构性能曲线。随着阵列阵元数目m的增加，所有方法的性能都得到了提高，包括rmse和srr两种性能。另一方面，online-ksbl的运行时间与m呈正相关。值得注意的是，msbl的运行时间曲线出现了一种违反直觉的情况。这可以解释为:虽然单次迭代的计算量较高，但m越大收敛速度越快。
[0133]
从图3(d)-(f)所示的稀疏性曲线来看，其趋势基本符合预期。online-ksbl和msbl的性能随着稀疏度的增加而恶化。由于其非迭代性质，online-ksbl的运行时间与稀疏度无关。而msbl需要更多的迭代次数才能达到收敛。
[0134]
对于确定数值初始化的online-ksbl，如图3(g)和(h)所示，由于参数逐渐收敛，随着测量值t的增大，估计精度逐渐提高；最优初始化online-ksbl的rmse和srr曲线变化不大，因为在初始化阶段就完成了对超参数的高精度估计。而msbl的重构性能与t无关。由于源向量按时序进行更新，online-ksbl的计算复杂度与t近似成正比。
[0135]
online-ksbl和msbl在低信噪比情况下的性能快速下降，如图3(j)-(l)所示。由于在低信噪比的情况下需要增加迭代次数，msbl的计算成本更高。然而，由于online-ksbl具
有非迭代性质，其运行时间与信噪比无关。
[0136]
总结图3(a)
‑‑
图3(l)的一些共性结论。最优初始化online-ksbl的重构精度优于msbl，而确定数值初始化online-ksbl的计算效率更高。online-ksbl的计算复杂度仅与测量数据t成正比，而msbl的计算量依赖于问题本身和源向量，具有不可预测性且计算量更高。
[0137]
下面通过水下多波束测深仪实验验证算法的成像性能。该实验在湖泊中进行，一个ping的回声数据由430个快拍组成。实验中，将两个直径为20cm的球形聚乙烯目标分别固定在深度约27m和28m处。图4(a)
‑‑
图4(d)为四种对比方法的成像结果，包括本发明提出的online-ksbl，多快拍msbl，动态稀疏编码(dsc)，和时间多快照稀疏贝叶斯学习(tmsbl)。其中，取δ
′
＝3和δ＝150。从msbl的成像结果看，图像在各快拍是分层的，这是由于在重建源向量时以δ
′
为周期独立地估计超参数，结果导致重构结果不一致性而且噪声更大。dsc算法重构的图像呈粒状结构，噪声增大。而tmsbl无法重建点目标，而且丢失目标细节。通过比较，证明了online-ksbl在成像质量方面是最好的。
[0138]
本发明还提出一种电子设备，包括存储器和处理器，所述存储器存储有计算机程序，所述处理器执行所述计算机程序时实现所述线上贝叶斯压缩水下成像方法的步骤。
[0139]
本发明还提出一种计算机可读存储介质，用于存储计算机指令，所述计算机指令被处理器执行时实现所述线上贝叶斯压缩水下成像方法的步骤。
[0140]
以上对本发明所提出的一种线上贝叶斯压缩水下成像方法及系统、设备和介质进行了详细介绍，本文中应用了具体个例对本发明的原理及实施方式进行了阐述，以上实施例的说明只是用于帮助理解本发明的方法及其核心思想；同时，对于本领域的一般技术人员，依据本发明的思想，在具体实施方式及应用范围上均会有改变之处，综上所述，本说明书内容不应理解为对本发明的限制。