简介：面向遥感与雷达信号处理学习者的极化SAR特征提取工程源码，聚焦全极化数据H/A/α分解的实现与复现，适合希望从代码层面理解极化分解流程、补齐SAR特征提取实践环节的读者。压缩包共17个文件、约29KB，以c与h源文件为主体，包含分解主程序、矩阵运算、ENVI配置读取与图像输出等模块，配合txt说明和VC6工程及调试配置文件，可直接打开编译并对照理解算法结构。程序基于T3矩阵计算H、A、α分量，并将结果整理为可用于后续分析的影像数据，可服务于地物分类、地表参数估计、变化检测及灾害监测等场景的前期实验。已有1658人学习，适合作为极化SAR特征提取的入门参考与二次开发基础。

1. 一景四极化数据摆上桌：极化SAR特征提取真正难的是算之前

拿到一景四极化单视复数数据，很多人的第一反应是找公式、凑特征，把 H/A/α、三分量分解、各种极化比值一口气堆到四五十维，再丢进分类器。跑完发现精度还不如只用三个强度通道——极化SAR特征提取里最常见的一类失败就是这个。问题几乎不出在特征公式上，而出在它上游：定标系数有没有套、C3 和 T3 用哪个基、多视够不够、滤波之后协方差矩阵还半不半正定。下面按一条能跑的链路拆开：从散射矩阵构建 C3/T3，到特征值分解和三分量分解，再到滤波、归一化、特征筛选，最后用空间分块验证这套特征到底值不值。适合已经能读写复数 SAR 数据、准备自己做特征工程的人，也适合带课题、需要一个能交付流程的从业者。

2. 数据入口决定上限：散射矩阵到 C3/T3 的构建与定标

特征提取的所有下游问题，几乎都能追溯到矩阵构建这一步。同一景数据，C3 用词典基、T3 用 Pauli 基，两者差一个酉变换；酉变换不改变特征值，所以 H/A/α 是一样的，但矩阵元素本身作为特征喂给分类器时，含义完全不同。把 C3 的九个实数当成特征、又把 T3 的九个实数拼进去，等于把同一组信息换了个坐标重复注入两遍，还会让后面基于相关性的特征筛选彻底失效。所以中间章从最不起眼的地方开始：先把矩阵建对。

2.1 先确认三件事：定标、基、元素排列

在写任何一行特征代码之前，值得花十分钟翻元数据。常见做法是核对下面四项，缺一项后面就会以“精度莫名其妙掉”的形式还回来。

第三条最容易被忽略。词典基下的向量是 $[S_{HH}, \sqrt{2}S_{HV}, S_{VV}]^T$，Pauli 基下是 $\frac{1}{\sqrt{2}}[S_{HH}+S_{VV}, S_{HH}-S_{VV}, 2S_{HV}]^T$。两者的第三个分量一个是交叉极化本身，一个是交叉极化的两倍，量纲差一个常数，伸进分类器后权重会被这个常数带偏。

第四条则是纯工程坑。有些导出把复数当成两个波段存成 (实部, 虚部) 交错，有些则按“6 个上三角元素的实部、6 个虚部”平铺成 12 个波段。后者读进来如果直接reshape，矩阵会静默错位——不报错，但相位特征全乱。

提示：先取一小块（比如 200×200 像元）做全流程验证，确认矩阵 Hermitian、特征值非负、σ0 量级合理，再跑全图。整景跑一遍再回头查错，成本高得多。

2.2 复散射矩阵到 C3/T3：一段能跑通的构建代码

假设已经拿到定标后的四个复通道，且满足互易条件（$S_{HV} = S_{VH}$）。下面这段把逐像元的散射矩阵组织成 C3，再变换到 T3。

import numpy as np def scattering_to_c3(hh, hv, vh, vv, reciprocal=True): """hh/hv/vh/vv: 复数二维数组，已定标并开方，单位是复幅度""" if reciprocal: # 互易介质下两个交叉极化理论上相等，取平均可以压掉通道间相位偏差 hv = 0.5 * (hv + vh) # 词典基：sqrt(2) 是为了让总功率等于四个通道功率之和 s = np.stack([hh, np.sqrt(2.0) * hv, vv], axis=0) # (3, H, W) c3 = s[:, None, :, :] * np.conj(s[None, :, :, :]) # 外积 -> (3, 3, H, W) return c3 # T3 = U3 @ C3 @ U3^H，U3 为 Pauli 变换矩阵 U3 = np.array([[1, 0, 1], [1, 0, -1], [0, np.sqrt(2), 0]], dtype=np.complex128) / np.sqrt(2) def c3_to_t3(c3): t3 = np.einsum('ij,jkab,lk->ilab', U3, c3, np.conj(U3)) # 数值上强制 Hermitian，避免后续 eigh 因微小非对称报错 return 0.5 * (t3 + np.conj(np.transpose(t3, (1, 0, 2, 3)))) def boxcar_multilook(mat, looks_r, looks_a): """mat: (3, 3, H, W) 复数；末尾不足一块的像元直接裁掉""" _, _, h, w = mat.shape ha, wa = h // looks_a, w // looks_r m = mat[:, :, :ha * looks_a, :wa * looks_r] m = m.reshape(3, 3, ha, looks_a, wa, looks_r).mean(axis=(3, 5)) return m

scattering_to_c3里的sqrt(2)不是装饰。词典基向量的模平方决定了 C3 对角线是不是总功率：$C_{11}=|S_{HH}|^2$，$C_{22}=2|S_{HV}|^2$，$C_{33}=|S_{VV}|^2$，三者相加正好是 SPAN。少了这个系数，后面用 C3 元素算出来的交叉极化比例会整体偏高或偏低，是那种“看着像对的、但阈值永远调不准”的错误。

c3_to_t3用einsum一次性完成逐像元的相似变换，写法上把通道维放在最后两个轴，是为了避免把 (H, W) 展平后再还原的样板代码。这种写法的代价是内存：双精度复数 3×3×H×W，一景 5000×5000 的图大约占 3.6 GB。工程上一般用complex64存中间结果，在做特征值分解前再转complex128。

boxcar_multilook的looks_r/looks_a是距离向和方位向的视数。多视不是可选项：单视复数数据的相干斑服从 Rayleigh 分布，直接把单视协方差矩阵做特征值分解，H 会系统性地偏低，α 会向 45° 附近塌缩。用视数 L 的块平均后，等效视数 ENL ≈ L，斑点标准差下降约 $1/\sqrt{L}$。

2.3 多视数与像元间距：别让特征图和标签错位

多视比例怎么定，取决于原始像元间距，而不是凭经验拍一个 5×5。机载或星载单视复数据，方位向像元间距常常远小于距离向（典型是 1.5 m 对 5 m），所以多视矩阵一般不平权，常见配比是距离向多看、方位向少看，比如 1×3、2×5。若强行做 5×5，结果是距离向被过度平滑、方位向分辨率白白丢掉。

这里有个血泪经验：多视之后像元尺寸变成间距 × 视数，如果标签是从另一份底图矢量化来的，配准必须走地理坐标，不能按行列号直接叠加。行列号对齐在原始分辨率下勉强能看，多视之后偏移会被放大几倍，最后表现为“分类结果整体往一个方向飘了一条带”，而精度表上只掉两三个点，很难察觉。

另一个值得记住的边界是 ENL 与特征稳定性的关系。ENL 低于 4 时，H 的直方图会明显向低端偏移，A 在植被区几乎全是噪声。此时要么继续多视，要么改用对相干斑不敏感的特征（比值类、强度类）打底，把 H/A/α 当辅助，别让它主导分类。

注意：多视会把空间分辨率永久降下来，且不可逆。如果后续还想做细小地物（道路、单栋建筑）的分割，宁可保留一份单视数据，多视结果只用来提统计特征，两者通过坐标索引关联。

3. 极化特征的三条主线：H/A/α、三分量分解与比值特征

矩阵建好之后，特征大致分三档：特征值分解类、模型分解类、比值与相关系数类。三档的稳定性、物理含义、计算代价完全不同，混在一起堆维度是最常见的做法，也是最容易吃亏的做法。下面逐档说清楚它们各自擅长什么、什么时候会失效。

3.1 H/A/α：三个特征值撑起的旋转不变特征

对 T3 做特征值分解，得到 $\lambda_1 \ge \lambda_2 \ge \lambda_3 \ge 0$ 和对应特征向量。归一化后得到伪概率 $P_i = \lambda_i / \sum \lambda$，于是：

$$H = -\sum_{i=1}^{3} P_i \log_3 P_i,\quad A = \frac{\lambda_2 - \lambda_3}{\lambda_2 + \lambda_3},\quad \bar{\alpha} = \sum_i P_i \alpha_i$$

其中 $\alpha_i = \arccos(|v_{1i}|)$，$v_{1i}$ 是第 i 个特征向量在 Pauli 基第一个分量（$S_{HH}+S_{VV}$）上的投影。这一组特征的物理意义很干净：H 描述散射随机性，A 描述两个较弱分量的相对强弱，α 描述主导散射机制（0° 附近是表面散射，45° 附近是体散射，90° 附近是二次散射）。

def h_alpha_anisotropy(t3): """t3: (3, 3, H, W) 复数 Hermitian 矩阵，建议先多视""" t3 = 0.5 * (t3 + np.conj(np.transpose(t3, (1, 0, 2, 3)))) lam, vec = np.linalg.eigh(t3.transpose(2, 3, 0, 1)) # 升序，形状 (H, W, 3) lam = np.clip(lam, 0.0, None) lam = lam[..., ::-1] # 降序 l1 >= l2 >= l3 vec = vec[..., ::-1] total = lam.sum(axis=-1, keepdims=True) + 1e-12 p = lam / total # p 可能为 0，取对数前先夹住，否则整幅图出 NaN h = -(p * np.log(np.clip(p, 1e-12, None))).sum(-1) / np.log(3.0) with np.errstate(invalid='ignore', divide='ignore'): a = (lam[..., 1] - lam[..., 2]) / (lam[..., 1] + lam[..., 2] + 1e-12) a = np.where(lam[..., 1] + lam[..., 2] > 1e-10, a, 0.0) # alpha_i 由特征向量首分量的模决定，取绝对值是因为特征向量符号可翻转 alpha_i = np.degrees(np.arccos(np.clip(np.abs(vec[..., 0, :]), 0.0, 1.0))) alpha = (p * alpha_i).sum(-1) return h, a, alpha, lam

np.linalg.eigh只对 Hermitian 矩阵保证实数特征值，所以第一行强制对称化不是保险，是必需。滤波、多视、重采样之后矩阵会有 1e-7 量级的非对称，eigh可能直接抛错，也可能悄悄用下三角算出一个不完整的结果。

np.clip(lam, 0, None)处理的是另一种常见状况：协方差矩阵在理论上半正定，但经过逐通道滤波、整型量化、辐射定标误差之后，最小特征值可能变成 -1e-5 这样的小负数。夹到 0 之后 H 仍然有意义，否则 $\log$ 直接崩。

α 的计算里取np.abs是必要的。特征向量差一个相位因子不改变物理结果，但首分量可能整体反号，不取模会让 α 在 180° − α 之间跳变，图上表现为大片黑白翻转的斑块。

这组特征的短板在低 ENL 区很突出。ENL 小于 4 时，λ2 和 λ3 已经和数值噪声同量级，A 基本是随机数；此时 A 应该在特征集里被降权甚至剔除，而不是当成一个和 H 平权的特征。判断方法很土但有效：把 A 单独画出来看，如果在同质农田里它也没有空间结构，说明它在这个数据上就是噪声。

3.2 Freeman-Durden 三分量：残差功率才是诊断量

三分量模型把总散射拆成表面散射、二次散射和体散射，前提是反射对称（$\langle S_{HH}S_{HV}^\rangle = \langle S_{VV}S_{HV}^\rangle = 0$）。在反射对称下，交叉极化只来自体散射，于是可以直接反解：

def freeman_durden(c3): """c3: (3, 3, H, W) 复数；返回三个分量功率与模型残差""" c11 = np.real(c3[0, 0]) c22 = np.real(c3[1, 1]) c33 = np.real(c3[2, 2]) c13 = c3[0, 2] span = c11 + c22 + c33 + 1e-12 hv2 = c22 / 2.0 # C3 用 sqrt(2) 加权，这里还原 <|S_HV|^2> fv = 3.0 * hv2 # 体散射系数 fv = np.clip(fv, 0.0, np.minimum(c11, c33)) # 物理约束：体积分量不能超过总功率 half = 0.5 * (c11 + c33) - fv # fs + fd diff = np.real(c13) - fv / 3.0 # fs - fd（取 alpha=-1, beta=1 的经典假设） fs = np.clip(0.5 * (half + diff), 0.0, None) fd = np.clip(0.5 * (half - diff), 0.0, None) ps = 2.0 * fs # 1 + |beta|^2 pd = 2.0 * fd # 1 + |alpha|^2 pv = 8.0 / 3.0 * fv # 体散射在球面上积分后的总功率 resid = np.abs(c11 - c33) / span # 模型不自洽程度的代理量 return ps, pd, pv, resid

hv2 = c22 / 2.0这一步是把自己坑过好几次的地方。C3 的第二对角元是 $2\langle|S_{HV}|^2\rangle$，而模型里 $\langle|S_{HV}|^2\rangle = f_v/3$。少了这个除以 2，体散射分量会整体偏大一倍，植被区 Ps/Pd/Pv 三者比例完全变样，但图上看起来还挺“像那么回事”。

np.clip(fv, 0, min(c11, c33))是第二个约束。原始模型反解出的 fv 在强点目标或定标异常处可能超过总功率，导致 fs、fd 变负。强行截断不是数学上最优的做法，但能保证输出的三个分量都是非负功率，可以直接当特征用。

真正值得当特征用的，反而是那个resid。反射对称假设在朝向建筑、陡坡、城市街区这些地方必然失效，此时 C11 和 C33 差异显著，残差抬升。把残差当成一个“模型可信度”通道喂进分类器，比单纯给三个分量功率更有信息量——它某种程度上编码了地物的朝向结构。

提示：如果你发现某类地物的 Pv 占比高得离谱而实际是裸地，先查残差图。残差高的区域，三分量的物理解释本身就不可靠，这时候该做的是朝向角补偿或换四分量模型，不是继续调分类器参数。

3.3 比值与指数特征：便宜且抗噪的那一档

第三档特征是比值和相关系数。它们不需要特征值分解，不需要迭代求解，对相干斑的敏感度也低得多，在特征集里应该占据基础位置。

def ratio_features(c3, floor=1e-4): c11, c22, c33 = np.real(c3[0, 0]), np.real(c3[1, 1]), np.real(c3[2, 2]) c13 = c3[0, 2] hh_vv = 10 * np.log10(np.clip(c11 / (c33 + floor), floor, None)) hv_vv = 10 * np.log10(np.clip(c22 / (c33 + floor), floor, None)) rvi = 8 * c22 / (c11 + c33 + 2 * c22 + floor) rho = np.abs(c13) / np.sqrt((c11 + floor) * (c33 + floor)) rho = np.clip(rho, 0.0, 1.0) phase = np.angle(c13) # 弧度，-pi 到 pi return hh_vv, hv_vv, rvi, 1.0 - rho, phase

比值类特征最大的好处是系统增益自动抵消。辐射定标如果整体差一个常数因子，比值完全不受影响；这也是为什么在定标信息存疑的数据上，先看比值特征往往比先看绝对 σ0 更可靠。

但phase是个例外。相位对相干斑极其敏感，ENL 低于 8 时基本是均匀分布在 $[-\pi, \pi]$ 的噪声，喂进分类器只会稀释有效维度。要用相位，先保证多视足够，或者用多视后的复相干矩阵做极化相位优化，而不是直接取逐像元相位。

到这一步，特征池大概长这样：σ0 三个、比值三个、RVI 一个、相关系数模一个、H/A/α 三个、三分量三个加残差一个、T3 的九个实元素。总共二十多维，已经足够。后面两章的任务是把这二十多维收拾干净，并确认哪些是真的有用。

4. 斑点与维度：把特征矩阵收拾干净

特征算出来只是半成品。斑点噪声会让同一地物的特征值分布严重重叠，量纲差异会让基于距离的分类器直接退化，冗余维度则会让小样本训练出的模型方差爆炸。这三件事各自有成熟做法，但都有细节。

4.1 极化滤波不能逐通道做：半正定投影是必要的补救措施

最省事的滤波当然是对每个通道单独做 Lee 或中值，但这样做会破坏矩阵的代数结构。逐通道滤波后，C12 和 C21 不再互为共轭，对角线也不再满足功率约束，特征值分解出来的 λ 可能为负，H 直接溢出。常见的“滤波后特征图出现雪花点”就是这么来的。

正确的做法只有两个方向：要么在矩阵域做多视或极化 Lee，滤波对象始终是整个 3×3 矩阵；要么滤波后强制投影回半正定锥。

def psd_project(t3): """把数值上可能非半正定的矩阵投影回半正定锥""" t3 = 0.5 * (t3 + np.conj(np.transpose(t3, (1, 0, 2, 3)))) lam, vec = np.linalg.eigh(t3.transpose(2, 3, 0, 1)) lam = np.clip(lam, 0.0, None) # 负特征值直接砍掉 out = np.einsum('...ij,...j,...kj->...ik', vec, lam, np.conj(vec)) return out.transpose(2, 3, 0, 1) def refined_lee_c3(c3, win=7): """退化版：矩阵域均值滤波 + 半正定投影，够用作对比基线""" k = np.ones((1, 1, win, win), dtype=np.float32) / (win * win) sm = np.stack([np.stack([ np.convolve(c3[i, j].ravel(), k.ravel(), mode='same').reshape(c3.shape[2:]) for j in range(3)], axis=0) for i in range(3)], axis=0) return psd_project(sm)

np.einsum('...ij,...j,...kj->...ik', vec, lam, np.conj(vec))实现的是 $V \Lambda V^H$。写成einsum而不是matmul，是为了让中间维（也就是 3）不参与广播，避免把 (H, W, 3) 当成批量维度后再转置回来时出错。

refined_lee_c3只是基线，真正的极化 Lee 需要按局部统计量自适应加权，但它已经能说明问题：只要滤波在矩阵域做，输出就还是一个合法的协方差矩阵，后面的特征值分解不会出异常值。我一般会先用这个基线跑一遍，确认特征正常，再考虑换更精细的滤波器。

注意：滤波窗口不是越大越好。7×7 在 5 m 分辨率下相当于 35 m，已经开始跨越田块边界；10 m 分辨率下 7×7 还算克制。窗口超过地物尺度，边缘特征就会被抹平，分类图上田埂、道路这类线状目标直接消失。

4.2 归一化：dB 域、标准化与 H/A/α 的例外

σ0 在典型场景里能跨四到五个数量级，水体和建筑区的线性功率差着十万倍。直接标准化（减均值除标准差）会被少数极值点主导，均值和方差都不可信。常规处理是先转 dB，再标准化。

import numpy as np from sklearn.pipeline import Pipeline from sklearn.preprocessing import StandardScaler from sklearn.ensemble import RandomForestClassifier def to_db(x, floor=1e-6): return 10.0 * np.log10(np.clip(x, floor, None)) # X: (n_samples, n_features)，前 8 列是线性功率域的量 X = X.copy() X[:, 0:8] = to_db(X[:, 0:8]) # H/A/α、相关系数模、残差本身已经在有界区间，保持不变 # 只在训练集上 fit，避免测试集信息泄漏 pipe = Pipeline([ ('scaler', StandardScaler()), ('clf', RandomForestClassifier(n_estimators=300, min_samples_leaf=2, class_weight='balanced_subsample', random_state=0, n_jobs=-1)), ]) pipe.fit(X_train, y_train)

floor=1e-6是为了对付零值像元（无效值填充、水体极端平滑区）。不加这个下限，log10(0)会产生-inf，标准化后整列变成NaN，而很多分类器对NaN不会报错，只会静默给出一个很差的模型。

StandardScaler一定要在管道里、只用训练集 fit。手动对整幅图做(X - mean) / std是最常见的精度虚高来源，尤其在样本量小的时候能虚出十几个点。

H/A/α 这一组不需要额外处理。H 和 A 天然在 [0, 1]，α 在 [0, 90]，标准化只会引入无谓的数值抖动。如果你的特征集把有界特征和无界特征混在一起标准化，问题不大，但解释特征重要性时要注意尺度已经被改过了。

4.3 特征筛选：先砍相关性，再看互信息

二十多维持平，但真正独立的信息维度远没有这么多。T3 的九个元素和 H/A/α、三分量之间高度共线，硬塞进去只会让树模型在相近特征间随机选一个，重要性排序变得不稳定，换一批样本结论就变。

筛选分两步：先按相关系数做硬剪枝，再用互信息排序。

import numpy as np from sklearn.feature_selection import mutual_info_classif def prune_by_corr(X, names, thr=0.95): """贪心保留：按与标签的互信息降序，剔除高度相关的后来者""" keep, dropped = [], [] for idx in order: # order 为互信息降序的列索引 if all(abs(np.corrcoef(X[:, idx], X[:, k])[0, 1]) < thr for k in keep): keep.append(idx) else: dropped.append(idx) return keep, dropped mi = mutual_info_classif(X_train, y_train, discrete_features=False, random_state=0) order = np.argsort(mi)[::-1] keep, dropped = prune_by_corr(X_train, feat_names, thr=0.95) X_sel = X_train[:, keep]

thr=0.95是经验值。0.9 以下会开始误删——H 和交叉极化比在植被区的相关系数能到 0.85，但它们提供的信息并不完全重叠。0.98 以上则几乎不删东西。0.95 附近是常见取舍点。

mutual_info_classif的discrete_features=False必须显式写。默认值在不同版本间变过，如果被当成离散特征处理，连续的 dB 值会被分箱，互信息估计会严重失真。

剪完之后通常还剩十几到二十维。这时候再回头看消融实验，就会发现某些看着很“高级”的特征其实贡献为零。哪些能删、怎么验证，放到最后一章用数据说。

5. 五个常见失败场景：从负特征值到精度虚高

下面这几条都是能复现的、有明确现象和确定处理方式的失败场景。按出现频率从高到低排。

5.1 特征图出现 NaN 或整片黑

现象：H 或 α 的图上出现成片的零值或 NaN，位置集中在山体阴影、水体、以及图像边缘。

原因：主要有两个。一是总功率为 0 或接近 0 的像元（无效值填充、定标

本文还有配套的精品资源，点击获取