阵列声波测井仪接收的波形受井孔周围岩石性质的影响, 其慢度和频散特征等均与地层性质密切相关^[1]。准确、有效地提取模式波中的频散信息, 能够更好地获取与地层性质相关的信息。然而, 实测波形数据受噪声和模式波相互干扰等影响并不严格满足数学模型假设, 相比于参数化频散分析算法对理论模型误差的敏感性, 非参数频散分析算法无需假设信号模型且具有更好的鲁棒性, 更适合处理复杂的实测波形数据^[2]。

1898年, SCHUSTER^[3]利用傅里叶变换算法进行谱分析时提出了周期图的概念, 但该算法存在主瓣能量泄露和方差较大的问题。NOLTE等^[4]基于傅里叶变换理论提出了处理阵列谱相干函数的加权频谱相干法(WSS), 采用多点加权策略提高了算法精度及抗噪性, 但该方法不适用于处理单一频率下多个模式波共存的波形数据。前人对周期图进行了大量研究, 其中包括带通滤波器组的处理算法, 如最小方差无失真波束形成算法(Capon), 通过设计数据相关的滤波器进行谱估计^[5]。LI等^[6]提出了一种自适应有限脉冲滤波器算法, 也称为振幅相位估计法(APES), 最初用于目标范围的特征估计与合成孔径雷达成像。与最小方差无失真波束形成算法相比, 振幅相位估计法算法对幅度估计更加准确。曹正良等^[7]利用模拟和实测阵列声波测井数据比较了普罗尼(Prony)方法、同态处理方法和加权频谱相干法3种频散分析方法。王瑞甲等^[8]将矩阵束算法与加权频谱相干法相结合, 给出了一种多道波形信号的频散分析方案, 压制了加权频谱相干法处理中出现的周期性伪解。刘航等^[9-10]通过纵波的频散分析提出纵波频散谱处理方法, 并形成基于频散谱的碳酸盐岩储层评价方法。

阵列声波测井波场分离技术用于分离全波列中的不同模式波^[11]。宫昊等^[12]用慢度-时间相干法与线性预测滤波法的组合从全波列中提取较弱的反射波, 有效地解决了反射声波测井中成像结果的多解性和可靠性差的问题。实测数据处理中, 由于模式波相互干涉及噪声影响等, 在慢度-频谱图中存在模式波慢度-频率相关曲线不连续或无法有效分辨等问题。因此, 有必要将时、频域波场分离与频散提取方法相结合以实现波形数据的频散分析。首先, 采用长短时方差比值法^[13]准确提取初至, 然后, 将其作为时域滤波窗的起止时刻, 并实施时域波场分离。

虽然加权频谱相干法、振幅相位估计法和最小方差无失真波束形成法等算法的性能在以往文献^{[7, 14-15]}中已有阐述, 但几种算法的应用效果以及性能的综合对比尚不全面, 实测数据的频散信息提取效果也不理想。本文利用模拟数据与实测数据详细比较了傅里叶变换法(FTM)、加权频谱相干算法、振幅相位估计算法和最小方差无失真波束形成算法及其前后向算法等共6种频散信息提取算法的性能及应用效果, 并认识到将基于长短时方差比值的时、频域滤波的波场分离算法与非参数频散估计算法相结合, 可以显著改善频散信息的提取效果。

1 频散信息提取方法的对比分析

考虑由N个等间距d接收器组成的声波阵列, 它接收的波形信号是平面波与噪声线性叠加后的信号。若以第一个阵列声波接收器为参考, 则频率域波形信号的解析解可以表示为^[16]:

$y_n(\omega)=\sum\limits_{p=1}^P \alpha_p(\omega) \exp \left[-\mathrm{j} \omega s_p(n-1) d\right]+e_n(\omega)$

(1)

式中: n=1, 2, …, N表示接收器的编号; α_p和s_p分别为第p个模式波的复幅度及慢度; ω=2πf是圆频率, f为频率; e_n(ω)表示接收信号中的噪声以及第n个接收器中模式波之间的相互干涉。

1.1 频散提取算法的基本关系

振幅相位估计算法(APES)是根据噪声协方差矩阵的特征值分布与大小自适应地设计相应有限脉冲滤波器, 使得慢度为s的信号无失真地通过, 并对其它信号及噪声进行压制, 振幅相位估计算法及其前后向(FB)算法估计所得复振幅分别为^[6]:

$\alpha_{\mathrm{APES}}(s)=\frac{\boldsymbol{a}^H(s) \overline{\boldsymbol{Q}}^{-1}(s) \overline{\boldsymbol{g}}(s)}{\boldsymbol{a}^H(s) \overline{\boldsymbol{Q}}^{-1}(s) \boldsymbol{a}(s)}$

(2)

$\alpha_{\mathrm{FBAPES}}(s)=\frac{\boldsymbol{a}^H(s) \boldsymbol{Q}^{-1}(s) \overline{\boldsymbol{g}}(s)}{\boldsymbol{a}^H(s) \boldsymbol{Q}^{-1}(s) \boldsymbol{a}(s)}$

(3)

式中: Q, Q分别为仅前向、前后向数据噪声协方差矩阵; a(s)为导向向量; g(s)为数据向量的归一化傅里叶变换, 算法细节见附录A及附录B。

最小方差无失真波束形成算法与振幅相位估计算法均为匹配滤波器频谱估计算法, 振幅相位估计法算法是在保证对期望信号响应一定(W^Ha(s)=1)的前提下, 最小化输出总功率, 其中, W为滤波器向量。令$\widetilde{\boldsymbol{R}}$表示数据向量协方差矩阵, 则滤波器输出功率为$W^H \widetilde{\boldsymbol{R}} W$, 具体求解过程可参考附录B。利用拉格朗日乘子法求得最小方差无失真波束形成算法及其前后向算法(FBCapon)估计的复幅度分别为^[17]:

$\alpha_{\text {Capon }}(s)=\frac{\boldsymbol{a}^H(s) \overline{\boldsymbol{R}}^{-1}(s) \overline{\boldsymbol{g}}(s)}{\boldsymbol{a}^H(s) \overline{\boldsymbol{R}}^{-1}(s) \boldsymbol{a}(s)}$

(4)

$\alpha_{\text {FBCapon }}(s)=\frac{\boldsymbol{a}^H(s) \boldsymbol{R}^{-1}(s) \overline{\boldsymbol{g}}(s)}{\boldsymbol{a}^H(s) \boldsymbol{R}^{-1}(s) \boldsymbol{a}(s)}$

(5)

式中: R, R分别为仅利用前向或前后向数据协方差矩阵。

加权频谱相干法包括频率域相关及加权两个步骤^[4]。将全波列波形数据d_n(t)变换到频率域{|D_n(ω)|}(n代表接收器编号), 第n个接收器与第n+1个接收器之间的关系为:

$D_{n+1}(\omega)=D_n(\omega) \exp \left(-a^0 d\right) \exp (-\mathrm{i} \omega s d)$

(6)

式中: a⁰表示频率ω下模式波的衰减系数; s表示其模式波慢度; d为接收器间距。(6)式中, 第一个指数部分代表衰减, 第二个指数部分代表波的传播。若定义向量D(ω)与X(ω, s)为:

$\boldsymbol{D}(\omega)=\left[\begin{array}{llll}D_1(\omega) & D_2(\omega) & \cdots & D_n(\omega)\end{array}\right]^{\mathrm{T}}$

(7)

$\boldsymbol{X}(\omega, s)=\left[\begin{array}{llll}1 & \mathrm{e}^{-\mathrm{i} \omega s d} & \cdots & \mathrm{e}^{-\mathrm{i} \omega s(n-1) d}\end{array}\right]^{\mathrm{T}}$

(8)

式中: X(ω, s)中的元素代表各接收器相对于第一个接收器的相移。则频率慢度相关归一化函数为:

$f(\omega, s)=\frac{\left|\boldsymbol{D}^*(\omega) \boldsymbol{X}(\omega, s)\right|}{\sqrt{\boldsymbol{D}^*(\omega) \boldsymbol{D}(\omega)} \sqrt{\boldsymbol{X}^*(\omega, s) \boldsymbol{X}(\omega, s)}}$

(9)

式中: 符号*表示共轭转置。对(9)式进行加权得到加权相干系数, 即

$F\left(\omega_k, s\right)=\sum\limits_{j=k-l}^{k+l} W\left(\omega_j\right) f\left(\omega_j, s\right)$

(10)

式中: W(ω_j)为权函数; l为加权的波谱点半径。可使用高斯函数作权函数, 在被加权的频率点处具有最大的权重, 即

$W\left(\omega_j\right)=\exp \left[-\frac{\left(\omega_k-\omega_j\right)^2}{2 \sigma^2}\right]$

(11)

式中: σ=N_wΔω, N_w表示参与加权的数据点数, Δω代表频率间隔。仅单点加权(N_ω=1)时, 加权频谱相干算法将退化为(9)式。

傅里叶变换法(FTM)利用傅里叶变换算法得到频率域波形信号, 消除由于传播带来的相位差异即可得到近似的复幅度估计。在实际数据处理中, (9)式分子部分即为FTM。

1.2 提取方法性能对比分析 1.2.1 保幅性及分辨率

根据阵列声波波形信号的解析解(1)式, 设计慢度分别为50 μs/ft, 60 μs/ft, 80 μs/ft 3个模式波的模型(1 ft≈30.48 cm)。模拟数据处理结果如图 1a所示, 黑色圆点为模式波实际慢度及振幅, 星标(曲线极值点)为算法估计的慢度及幅度; 振幅相位估计算法、前后向振幅相位估计算法以及前后向最小方差无失真波束形成算法估计的80 μs/ft模式波幅度相近; 加权频谱相干算法及傅里叶变换算法不能有效识别邻近较弱模式波, 弱模式主瓣被强模式旁瓣所淹没, 两者均是基于傅里叶变换理论来估计信号所包含的慢度成分, 它是周期函数截断的结果, 造成截断边界处存在能量泄漏, 导致估计结果失真、分辨率下降。振幅相位估计类算法及最小方差无失真波束形成类算法以目标慢度信号保真为前提设计的滤波器, 能准确地估计模式波的慢度。振幅相位估计算法及前后向振幅相位估计算法均可准确估计信号幅度和相位。最小方差无失真波束形成算法估计的信号幅度普遍偏低, 前后向最小方差无失真波束形成算法比仅利用前向数据的最小方差无失真波束形成算法幅度估计精度更高。

对辽河某地实测数据进行频散分析, 6种算法的实际数据处理(2.41 kHz处)的最大分辨率对比结果如图 1b所示。此处仅进行分辨率对比, 所以将所有扫描慢度点的幅度值均除以幅度最大值进行归一化。最小方差无失真波束形成类频散分析算法无旁瓣、能量集中、分辨率最佳, 振幅相位估计类算法分辨率次之, 傅里叶变换法及加权频谱相干算法分辨率最低。

1.2.2 抗噪能力

声波测井实测数据不可避免地含有噪声。在单频(8 kHz)单模式波(80 μs/ft)波形解析解(1)式信号中加入随机高斯白噪声, 以慢度估计的相对误差衡量算法抗噪性。模拟数据接收器总数设为13个, 最小方差无失真波束形成算法、振幅相位估计算法及其前后向算法的滤波器权系数个数均为6, 曲线每点均为1 000次重复实验的相对误差平均结果。图 2显示了非参数算法的抗噪能力对比结果。由图 2可见, 最小方差无失真波束形成算法和振幅相位估计算法对应的前后向算法抗噪性能得到提升; 傅里叶变换算法与单点加权的加权频谱相干算法具有极强的抗噪性, 两者慢度估计相对误差随信噪比变化曲线几乎重合, 抗噪性能相当; 多频点数据加权的频谱相干算法增大了数据量, 有效降低了噪声的影响。

1.3 提取方法应用效果

利用辽河某地实测单极子声波波形数据对比6种频散提取方法的应用效果。每个波列共有672个时间采样点, 采样间隔为16 μs, 8个接收器, 振幅相位估计算法与最小方差无失真波束形成算法以及它们的前后向算法滤波器权系数个数均为4。处理结果如图 3所示。由图 3可见, 傅里叶变换算法(图 3c)与加权频谱相干算法(图 3f)的处理结果中能量较强的旁瓣会对频散信息准确提取造成干扰, 但其模式波的慢度-频率曲线连续性好, 能量分布较为均衡。振幅相位估计算法的处理结果(图 3a)受噪声等影响较大, 存在如白圈中心处的能量局部异常点。前后向算法使得能量更加聚集, 多模式分析方法中的前后向最小方差无失真波束形成算法(图 3e)与前后向振幅相位估计算法(图 3d)处理效果相近, 单模式分析方法中的加权频谱相干算法处理效果更好(图 3f)。6种频散提取方法处理结果中模式波及伪解特征相似, 但均未有效识别出纵波以及伪瑞利波。

表 1为数值模拟计算采用的井孔流体和地层参数, 所用声源函数选择主频为12 kHz的Ricker子波, 最小源距为3.66 m, 13个换能接收器等间隔(0.152 4 m)排列。

衰减模型为^[18]:

$v(\omega)=v\left(\omega_{\text {ref }}\right)\left(1+\frac{1}{\pi Q} \log \frac{\omega}{\omega_{\text {ref }}}-\frac{\mathrm{i}}{2 Q}\right)$

(12)

式中: Q为60, ω_ref为12 kHz, 频散提取方法模拟数据应用效果对比如图 4所示。结果中虚线分别为井内流体慢度(203.2 μs/ft)以及地层纵、横波理论频散曲线。

由图 4可见, 傅里叶变换算法(图 4c)及加权频谱相干算法(图 4f)由于强信号能量的泄露会严重干扰较弱模式的识别, 所以更适用于单一模式波识别。即使在不含噪声的情况下, 由于算法分辨率及不同模式波之间的相互干涉等影响, 6种非参数频散信息提取方法也都不能同时且完整地提取所有模式波。

针对实际应用中频散信息提取效果变差的问题, 我们认为, 首先, 应进行基于长短时方差比值的时域滤波与频域滤波相结合实现波场分离。然后, 再进行频散提取处理。波场分离将多种模式波叠加的信号分解成单一模式, 降低模式波间的相互干涉及噪声等干扰对处理效果的影响; 同时, 也能使单模式频散分析方法较为准确地提取纵波、横波及伪瑞利波的频散信息。

2 模式波分离方法

不同模式波出现时会表现出明显的能量变化。定义反映瞬态能量变化和整体能量变化的“短时方差”及“长时方差”, 两者的比值可以突出初至到达时的能量变化, 从而检测纵、横波波至^[19](图 5)。设信号f(t)在第n点处的短时窗能量与长时窗能量的比值为^[13]:

$R_{n+1}=\frac{n \sum\limits_{i=0}^{l-1} f_{n+i}^2}{\sum\limits_{m=1}^n \sum\limits_{i=0}^{l-1} f_{m+i}^2}$

(13)

式中: l为短时窗的长度。

声波全波列中, 各模式波传播速度不同, 同时仪器较大源距使得时域波形分离较好, 利用长短时方差比值法提取初至从而完成时域波场分离, 再结合频域分离即可较好地分离目标模式波。波场分离使得单模式频散分析方法能够识别出被较强模式旁瓣淹没的较弱模式波, 图 6a为图 4对应的地层中单极子激发所得完整波形数据纵波分离后的频散分析结果, 即图 4f纵波分离后的结果, 白色虚线为纵波理论频散曲线。图 6b为与图 3相同的实测数据体的处理结果, 即图 3f纵波分离后的结果, 白色虚线为该深度的声波时差参考值(54.8 μs/ft)。处理结果对比表明, 基于长短时方差比值提取模式波初至与时域滤波结合, 再进行频域滤波的波场分离不改变目标模式波慢度及频谱分布特征。

3 频散信息提取

基于频散分析方法的性能及应用效果的综合对比可知, 最小方差无失真波束形成类算法抗噪能力强、模式波慢度估计精度及分辨率高, 受旁瓣影响极小, 且能够有效识别同一频率下多个强弱不同的模式波。前后向最小方差无失真波束形成算法相对仅利用前向最小方差无失真波束形成算法综合性能更好, 所以选择前后向最小方差无失真波束形成算法提取频散信息较为理想。

实测数据处理时, 若斯通利波能量轴连续、峰值强, 可直接提取; 否则, 利用长短时方差比值法提取纵、横以及斯通利波初至分离纵横波, 然后分别提取频散信息; 最后提取斯通利波及横波慢度范围内能量较强的多阶伪瑞利波(p-R), 得到所有模式波频散信息如图 7a所示。从图 7a中看到, 提取的频散曲线与细线的理论频散曲线吻合较好。图 7b给出了图 3中使用的辽河某地实测数据得到的频散曲线, 图中看到纵波、横波和斯通利波的频散曲线均连续性好, 图 3中不能提取的纵波模式也得到了较好的提取结果。

4 结论

基于非参数估计算法能够较好地处理模式波数量未知且慢度频散规律复杂的实测数据, 利用模拟数据及实测数据综合对比了傅里叶变换算法、加权频谱相干算法、振幅相位估计算法、前后向振幅相位估计算法、最小方差无失真波束形成算法和前后向最小方差无失真波束形成算法共6种非参数频散信息提取算法的保幅性、分辨率、抗噪性以及实际应用效果, 并与波场分离相结合提取了不同模式波的频散信息, 得到以下认识:

1) 加权频谱相干算法和傅里叶变换算法具有良好的抗噪能力, 但旁瓣大、保幅能力差且分辨率低; 振幅相位估计类算法具有较高分辨率, 振幅、相位估计优势明显, 但抗噪能力弱; 最小方差无失真波束形成算法振幅估计值偏低, 但其模式波慢度估计精度高, 同时具有良好的分辨率及抗噪能力。与仅使用前向数据的算法相比, 前后向算法的抗噪能力有所提高; 前后向最小方差无失真波束形成算法的幅值估计精度比最小方差无失真波束形成算法的估计精度高近一倍。因此, 建议在实测数据的速度频散信息分析中选用前后向最小方差无失真波束形成算法。

2) 在处理实测数据时, 振幅相位估计算法和最小方差无失真波束形成算法及其前后向4种多模式分析方法处理结果的慢度-频谱曲线连续性差, 峰值能量变化大; 单模式分析方法仅能分辨能量较强的模式波, 难以分辨慢度间隔较小的相邻模式且旁瓣能量强; 由于噪声、模式波的相互干涉及算法分辨率等影响, 6种频散信息提取方法均难以同时有效且完整地分辨出所有模式波。

3) 基于长短时方差比值提取初至并结合频域滤波的波场分离方法不改变模式波慢度及频散特征, 与频散分析方法组合使用可以降低噪声及模式波相互干涉的影响, 提升频散分析效果。而且, 波场分离后用单模式频散分析方法也能较好提取纵波和斯通利波等模式波的频散信息。

由于伪瑞利波与斯通利波和横波在时频域均有混叠, 基于长短时方差比值的时域与频域相结合的波场分离方法仍然不能有效分离伪瑞利波, 且单模式分析方法仅能识别同一频率下较强的信号, 所以加权频谱相干算法无法从模拟数据中提取得到伪瑞利波, 而前后向最小方差无失真波束形成算法由于其较高的分辨率可以提取伪瑞利波频散信息。

附录A 振幅相位估计法 A1 仅前向振幅相位估计算法(APES)

构建信号向量Y_n(1≤n≤N－M+1), 即

$\boldsymbol{Y}_n=\left[y_n, y_{n+1}, \cdots, y_{n+M-1}\right]^{\mathrm{T}}$

(A1)

将频域波形信号解析解(1)式代入上式, 有:

$\boldsymbol{Y}_n=\boldsymbol{A} s(n)+\boldsymbol{e}(n)$

(A2)

式中: A=[a(s₁)  a(s₂)  …  a(s_p)], 向量a(s), s(n)及e(n)分别定义为:

$\boldsymbol{a}(\boldsymbol{s})=\left[\begin{array}{c}1 \\ \mathrm{e}^{-\mathrm{j} \omega s d} \\ \vdots \\ \mathrm{e}^{-\mathrm{j} \omega s(M-1) d}\end{array}\right] \quad \boldsymbol{s}(\boldsymbol{n})=\left[\begin{array}{c}\alpha_1 \mathrm{e}^{-\mathrm{j} \omega s 1(n-1) d} \\ \alpha_2 \mathrm{e}^{-\mathrm{j} \omega s 2(n-1) d} \\ \vdots \\ \alpha_p \mathrm{e}^{-\mathrm{j} \omega s p(n-1) d}\end{array}\right]\quad \boldsymbol{e}(\boldsymbol{n})=\left[\begin{array}{c}e(n) \\ e(n+1) \\ \vdots \\ e(n+M-1)\end{array}\right]$

(A3)

将信号向量通过权系数个数为M的滤波器W=[w₁  w₂  …  w_M]^T, 则输出为:

$X_n(\omega)=\boldsymbol{W}^H \boldsymbol{Y}_n=\boldsymbol{W}^H \boldsymbol{a}\left(s_1\right) \alpha_1 \mathrm{e}^{-\mathrm{j} \omega s 1(n-1) d}+\boldsymbol{W}^H \cdot z_n$

(A4)

式中: z_n为其它慢度分量和噪声的叠加。

设计滤波器系数使慢度为s₁的信号无失真地通过, 同时尽可能压制Y_n中其它慢度分量和噪声, 使W^Hz_n平均功率最小, 即

$\min\limits_{\boldsymbol{W}} E\left\{\left|\boldsymbol{W}^H \cdot z_n\right|^2\right\} \quad \text { st. } \quad \boldsymbol{W}^H \boldsymbol{a}\left(s_1\right)=1$

(A5)

求解上述约束最小化问题(见附录B)可以得到滤波器矢量及信号幅度的估计值, 满足下式:

$W(s)=\frac{\bar{Q}^{-1}(s) \boldsymbol{a}(s)}{\boldsymbol{a}^H(s) \bar{Q}^{-1}(s) \boldsymbol{a}(s)}$

(A6)

$\alpha(s)=\frac{\boldsymbol{a}^H(s) \bar{Q}^{-1}(s) \overline{\boldsymbol{g}}(s)}{\boldsymbol{a}^H(s) \bar{Q}^{-1}(s) \boldsymbol{a}(s)}$

(A7)

$\bar{\boldsymbol{g}}(s)=\frac{1}{L} \sum\limits_{n=1}^L \boldsymbol{Y}_n \mathrm{e}^{\mathrm{j} \omega s(n-1) d}$

(A8)

式中: L=N－M+1, g(s)是数据向量的归一化傅里叶变换。也可以看作信号向量在基向量a(s)上的展开。令R为信号协方差矩阵:

$\overline{\boldsymbol{R}}=\frac{1}{L} \sum\limits_{n=1}^{N-M+1} \boldsymbol{Y}_n \boldsymbol{Y}_n^H$

(A9)

令Q(ω)为噪声协方差矩阵:

$\overline{\boldsymbol{Q}}(s)=\overline{\boldsymbol{R}}-\overline{\boldsymbol{g}}(s) \overline{\boldsymbol{g}}^H(s)$

(A10)

A2 前后向振幅相位估计算法(FBAPES)

对于后向振幅相位估计滤波器(BAPES), 只需要将公式(A1)信号向量Y_n改变为:

$\widetilde{\boldsymbol{Y}}_n=\left[y_{N-n+1}^*, y_{N-n}^*, \cdots, y_{N-n-M+2}^*\right]^{\mathrm{T}}$

(A11)

式中: [y]^*表示y的复共轭, [y]^T表示y的转置, 1≤n≤N－M+1。则后向信号协方差矩阵为:

$\tilde{\boldsymbol{R}}=\frac{1}{L} \sum\limits_{n=1}^{N-M+1} \widetilde{\boldsymbol{Y}}_n \tilde{\boldsymbol{Y}}_n^H$

(A12)

前后向信号协方差矩阵R为:

$\boldsymbol{R}=\frac{1}{2}(\overline{\boldsymbol{R}}+\tilde{\boldsymbol{R}})$

(A13)

定义后向数据向量的归一化傅里叶变换g(s)为:

$\tilde{g}(s)=\frac{1}{L} \sum\limits_{n=1}^{N-M+1} \widetilde{\boldsymbol{Y}}_n \mathrm{e}^{\mathrm{j} \omega s(n-1) d}$

(A14)

前后向噪声协方差矩阵与前向噪声协方差矩阵相似, 由(A10)式我们可以获得前后向噪声协方差矩阵Q:

$\boldsymbol{Q}(s)=\boldsymbol{R}-\boldsymbol{G}(s) \boldsymbol{G}^H(s)$

(A15)

其中,

$\boldsymbol{G}(s)=\frac{1}{\sqrt{2}}[\overline{\boldsymbol{g}}(s) \quad \tilde{\boldsymbol{g}}(s)]$

(A16)

附录B 约束优化问题的求解方法

用样本序列的时间平均代替(A5)式中的数学期望并将问题一般化, 可以得到:

$\min\limits_{\boldsymbol{W} \cdot a}\left\{H=\frac{1}{L} \sum\limits_{n=1}^L\left|\boldsymbol{W}^H \boldsymbol{Y}_n-\alpha \mathrm{e}^{-\text{j} \omega s(n-1) d}\right|^2\right\} \quad \text { st. } \quad \boldsymbol{W}^H \boldsymbol{a}(s)=1$

(B1)

令L=N－M+1, 将上式展开得到:

$H=\frac{1}{L} \sum\limits_{n=1}^L\left[\boldsymbol{W}^H \boldsymbol{Y}_n \boldsymbol{Y}_n^H \boldsymbol{W}-\boldsymbol{W}^H \boldsymbol{Y}_n \alpha^* \mathrm{e}^{\text{j} \omega s(n-1) d}-\alpha \mathrm{e}^{-\text{j} \omega s(n-1) d} \boldsymbol{Y}_n^H \boldsymbol{W}+|\alpha|^2\right]$

(B2)

利用等式$\boldsymbol{W}^H \overline{\boldsymbol{g}}=\overline{\boldsymbol{g}}^{\mathrm{T}} \boldsymbol{W}^*$可以将上式改为:

$H=\left|\alpha-\boldsymbol{W}^H \overline{\boldsymbol{g}}(s)\right|^2+\boldsymbol{W}^H\left[\overline{\boldsymbol{R}}-\overline{\boldsymbol{g}}(s) \overline{\boldsymbol{g}}^H(s) \boldsymbol{W}\right]$

(B3)

为使得H最小, 则在约束条件下W^Ha(s)=1有:

$\alpha=\boldsymbol{W}^H \overline{\boldsymbol{g}}(s)$

(B4)

以及$\boldsymbol{W}^H \overline{\boldsymbol{Q}} \boldsymbol{W}=\min$, 其中Q=R-g(s)g^H(s)可以看作噪声协方差矩阵。则(A5)式中最小化问题可以简化为条件极值问题, 如下所示:

$\min\limits_{\boldsymbol{W}} \boldsymbol{W}^{\mathrm{H}} \overline{\boldsymbol{Q}} \boldsymbol{W} \quad \text { st. } \quad \boldsymbol{W}^{\mathrm{H}} \boldsymbol{a}(s)=1$

(B5)

利用拉格朗日乘子法, 构造代价函数:

$J(\boldsymbol{W}, \lambda)=\boldsymbol{W}^H \overline{\boldsymbol{Q}} \boldsymbol{W}+\lambda\left(1-\boldsymbol{W}^{\mathrm{H}} \boldsymbol{a}\right)$

(B6)

式中: λ称为拉格朗日乘子。求代价函数梯度并令$\nabla J(\boldsymbol{W}, \lambda)=0$可得:

$\nabla J(\boldsymbol{W}, \lambda)=2 \overline{\boldsymbol{Q}} \boldsymbol{W}-2 \lambda \boldsymbol{a}=0$

(B7)

假设噪声和干扰的协方差矩阵非奇异, 则有

$\boldsymbol{W}=\lambda \overline{\boldsymbol{Q}}^{-1} \boldsymbol{a}$

(B8)

将上式代入约束条件W^Ha=1, Q为对称矩阵, 则可求得:

$\lambda=\frac{1}{\boldsymbol{a}^{\mathrm{H}} \bar{\boldsymbol{Q}}^{-1} \boldsymbol{a}}$

(B9)

则滤波器权系数有:

$W=\frac{\bar{\boldsymbol{Q}}^{-1} \boldsymbol{a}}{\boldsymbol{a}^H \bar{\boldsymbol{Q}}^{-1} \boldsymbol{a}}$

(B10)