2. 海杂波幅度分布拟合与参数估计

为什么海杂波不能只用瑞利分布

海杂波幅度分布随海况、擦地角、极化方式和雷达分辨率变化显著。低海况、低分辨率、大擦地角时，回波包络接近瑞利分布；当海尖峰（sea spike）出现时，尾部变厚，瑞利分布会严重低估大幅度回波概率，导致按瑞利假设设计的门限在实测数据上虚警率飙升。对数正态、韦布尔和K分布分别用不同的尾部建模能力来刻画这种非均匀性，其中K分布由散斑（快速起伏的复高斯分量）与调制分量（慢变的Gamma分布）复合而成，最贴合高分辨率海杂波物理机理。

四类分布的幅度PDF与参数含义

设回波幅度为 x ≥ 0，各分布的概率密度函数如下（均为幅度域形式）：

• 瑞利：p(x) = (x/σ²)·exp(−x²/(2σ²))，单参数 σ，σ² 为同相/正交分量方差，无形状参数。
• 对数正态：p(x) = 1/(x·σ·√(2π))·exp(−(ln x − μ)²/(2σ²))，μ 为 ln x 的均值（尺度），σ 为 ln x 的标准差（形状）。
• 韦布尔：p(x) = (k/λ)·(x/λ)^(k−1)·exp(−(x/λ)^k)，k 为形状参数（k=2 退化为瑞利），λ 为尺度参数。
• K分布：p(x) = (2/(a·Γ(ν+1)))·(x/(2a))^(ν+1)·K_ν(x/a)，ν 为形状参数（越小尾部越厚），a 为尺度参数，K_ν 为第二类修正贝塞尔函数。

工程上更常用的是强度域（功率）数据，此时瑞利对应指数分布，K分布的散斑分量对应Gamma分布。拟合前必须明确数据是幅度还是功率，否则参数会差一个平方关系。

参数估计：矩估计与最大似然

瑞利与对数正态有闭式矩估计；韦布尔可用对数矩法快速给出初值；K分布参数估计较难，工程上常用一阶/二阶矩加对数矩的混合矩法，或对数值似然做数值优化。下面给出可直接运行的Python实现，采用矩估计初值 + 最大似然精化的策略，并用K-S统计量与AIC做拟合优度比较。

import numpy as np
from scipy import stats, optimize, special

def fit_rayleigh(x):
    sigma = np.sqrt(np.mean(x**2) / 2.0)
    return {"sigma": sigma}

def fit_lognormal(x):
    lx = np.log(x)
    return {"mu": lx.mean(), "sigma": lx.std(ddof=0)}

def fit_weibull(x):
    # 对数矩法初值，再用MLE精化
    lx = np.log(x)
    m1, m2 = lx.mean(), lx.var()
    k0 = np.pi / np.sqrt(6*m2)          # 形状初值
    lam0 = np.exp(m1 + 0.5772/k0)       # 尺度初值
    def nll(p):
        k, lam = p
        if k <= 0 or lam <= 0: return 1e12
        return -np.sum(np.log(k/lam) + (k-1)*np.log(x/lam) - (x/lam)**k)
    res = optimize.minimize(nll, [k0, lam0], method="Nelder-Mead")
    return {"k": res.x[0], "lambda": res.x[1]}

def fit_k(x):
    # 混合矩法初值：二阶矩/四阶矩 + 对数矩
    m2 = np.mean(x**2); m4 = np.mean(x**4)
    ratio = m4 / m2**2
    nu0 = max(0.1, 2.0/(ratio - 2.0)) if ratio > 2 else 0.5
    a0 = np.sqrt(m2 / (4*nu0))
    def nll(p):
        nu, a = p
        if nu <= 0.05 or a <= 0: return 1e12
        logp = (np.log(2) - np.log(a) - special.gammaln(nu+1)
                + (nu+1)*np.log(x/(2*a))
                + np.log(special.kv(nu, x/a) + 1e-300))
        return -np.sum(logp)
    res = optimize.minimize(nll, [nu0, a0], method="Nelder-Mead",
                            options={"maxiter": 2000})
    return {"nu": res.x[0], "a": res.x[1]}

def ks_and_aic(x, cdf, n_params):
    xs = np.sort(x)
    ecdf = np.arange(1, len(xs)+1) / len(xs)
    D = np.max(np.abs(ecdf - cdf(xs)))
    ll = np.sum(np.log(np.clip(cdf(xs+1e-9)-cdf(xs-1e-9), 1e-300, None)))
    aic = 2*n_params - 2*ll
    return D, aic

拟合完成后，用K-S统计量D衡量经验分布与理论分布的最大偏差，用AIC在拟合优度与参数个数之间折中。判读原则：D越小越好；AIC越小越好；若K分布的D显著低于其余三者且AIC最小，则说明海尖峰尾部被有效刻画，应选K分布作为CFAR设计的杂波模型。

参数估计结果判读与CFAR衔接

以某X波段雷达、HH极化、3级海况、擦地角约5°的一段实测幅度数据为例（示例数据，非实测）：瑞利拟合得 σ≈0.42，K-S D≈0.118；对数正态 μ≈−0.95、σ≈0.78，D≈0.061；韦布尔 k≈1.35、λ≈0.92，D≈0.047；K分布 ν≈1.8、a≈0.55，D≈0.029，AIC最小。判读结论：尾部由海尖峰主导，K分布最优，ν≈1.8 表明调制分量起伏较强，后续CFAR应按K分布杂波背景设计门限，而不能沿用瑞利假设的门限系数。

参数估计的稳定性同样重要：把数据分成若干独立段重复拟合，观察 ν 与 a 的段间离散度。若 ν 的段间标准差超过其均值的20%，说明海况在扫描周期内发生变化，应缩短参数更新周期或改用自适应参考窗，而不是用单一参数覆盖整个扫描。