极化SAR特征提取：从相干矩阵到H/α工程实践

简介：面向遥感与SAR图像处理方向的学习者，这份代码包聚焦全极化数据的特征提取，围绕H、A、α三分量分解展开，可用于地物分类、地表参数估计、变化检测与灾害监测等场景的入门练习和算法复现，适合具备一定C语言基础、希望理解极化分解实现细节的读者。压缩包共17个文件，约29KB，以C源文件和头文件为主体，配合txt说明与工程配置文件，构成一套可编译的T3矩阵H/A/α分解程序骨架，目录简洁，便于按模块查阅极化数据处理流程。已有1658人浏览学习。读者可从中了解Cloude-Pottier分解的基本环节，理解H、A、α分量分别对应的相干性、散射不对称性与极化方向含义，并获得矩阵运算、配置读取等辅助代码及说明文档中的使用提示，为后续特征提取实验、课程设计和二次开发提供可直接参考的代码基础。

1. 极化SAR特征提取：多两个极化通道，信息量真的翻倍了吗

全极化 SAR 相比单极化只多了一个 HV 通道，但真正被打开的并不是通道数量，而是 HH 与 VV 之间的相对相位。可这个相位恰恰是整条链路里最容易被破坏的量。极化SAR特征提取之所以难，难在它不像光学纹理那样“看到什么算什么”：成像几何、极化定标、斑点噪声、多视窗口，每一层都会先改变相位，再让你去提特征。见过不少把特征堆到几十维、分类精度反而不如三维的案例，问题往往不在分类器，而在信息源头已经被污染，后面再精细的调参也只是安慰剂。把这件事的底层逻辑理清，比多背几个分解算法有用得多。

2. 散射矩阵到相干矩阵：特征的信息源头在哪里

极化 SAR 的原始观测是一个 2×2 的复散射矩阵，最终用于提特征的却几乎都是 3×3 的 Hermitian 矩阵。这个转换过程决定了后面所有特征的上限。如果转换方式错了，后面算熵、算 α 角、做分解，得到的数值看起来仍然很“正常”，但它们描述的已经不是地物，而是处理流程本身引入的偏差。

2.1 单视相干矩阵为什么是“废”的：秩 1 与多视平均

在互易假设下（单站 SAR 里 Shv ≈ Svh），一个像素的散射信息可以用三维 Pauli 向量表示：

k = 1/√2 · [Shh+Svv, Shh−Svv, 2Shv]ᵀ

对这个向量做外积 k·kᴴ，得到的就是单视相干矩阵 T3。它是一个 3×3 的 Hermitian 矩阵，有 9 个实自由度，但因为是单个向量的外积，秩恒为 1，9 个实参数里只有 5 个是独立的（6 个实分量减去一个全局相位）。这意味着：单视 T3 的两个特征值必然为 0，特征值谱退化成一根独苗，基于特征值分布的熵 H 永远等于 0，各向异性 A 永远是 0/0 的不定式。

常见的实现思路是：先做空间多视平均，用邻域内 L 个独立样本估计二阶统计量。平均之后 T3 才可能满秩，特征值谱才有形状可谈。代价同样明显：多视窗口越大，等效视数 L 越高，斑点噪声压得越狠，但空间分辨率丢得越多，边缘和小目标被窗口“化掉”。多视这一步本质上是在用空间分辨率换取统计稳定性；它要解决的是“要不要用空间换统计”这一权衡，并非要不要加一个滤波步骤的问题。

这里有个容易忽略的细节：多视窗口必须落在同一类地物内部才成立。跨越道路、水陆边界的窗口，会把两种散射机制平均成一个物理上不存在的中间态，这类像素在 H/α 平面上常常落在两个聚类之间，成为分类器最不稳定的样本。

2.2 Pauli 基与 Lexicographic 基：换坐标系不换信息

同一份散射信息，有几套常用的向量化方式。Lexicographic 基直接取 [Shh, √2·Shv, Svv]ᵀ，Pauli 基取 [Shh+Svv, Shh−Svv, 2Shv]ᵀ/√2。两套基之间相差一个固定的 3×3 酉矩阵，属于相似变换，因此特征值谱完全一致。

之所以工程上偏爱 Pauli 基，是因为它的三个分量各自对应一种典型散射机制：Shh+Svv 对应奇次散射（地表、水面），Shh−Svv 对应偶次散射（建筑二面角），Shv 对应体散射（植被冠层）。这样算出来的相干矩阵，对角线元素天然就是三种机制的功率，解释成本低。

但“特征值不变”这件事经常被误读成“特征不变”。相似变换保住的是谱，不保住特征向量的物理含义。如果拿 Pauli 基下算出的特征向量去做极化旋转不变性分析，再和 Lexicographic 基下的结果混着比较，就会得到一堆互相矛盾的数字。

2.3 特征值分解：3 个实数把 9 个自由度压成可解释量

相干矩阵的特征值分解 T3 = Σ λᵢ·uᵢuᵢᴴ，把 9 个实自由度整理成“3 个强度 + 3 个方向”的组合。λ₁ ≥ λ₂ ≥ λ₃ ≥ 0，特征值的大小顺序本身就代表了主导散射机制和次要机制的强弱关系；特征向量的第一个分量 |uᵢ₁| 则编码了该机制偏离奇次散射的程度。

import numpy as np def build_t3(shh, shv, svv, looks=1): """由单视 Pauli 向量构造 T3，并按 looks 做不重叠多视平均""" # Pauli 基：k = 1/sqrt(2) * [Shh+Svv, Shh-Svv, 2*Shv] k = np.stack([shh + svv, shh - svv, 2.0 * shv], axis=0) / np.sqrt(2.0) # 单视外积：秩 1，此时特征值只有一个是非零的 t_single = k @ k.conj().T if looks <= 1: return t_single, 1 # 多视：沿方位/距离方向做块平均，等效视数 = looks # 实际工程里通常用滑窗，这里用不重叠块保证 L 可控 h, w = shh.shape t_multi = np.zeros((3, 3, h // looks, w // looks), dtype=np.complex128) for i in range(looks): for j in range(looks): t_multi += t_single[:, :, i::looks, j::looks][:, :, :h // looks, :w // looks] return t_multi / (looks * looks), looks * looks def eigen_decompose(t3): """特征值降序排列，返回 (lambda, U)，U 的每一列是对应特征向量""" lam, u = np.linalg.eigh(t3) # eigh 返回升序，且保证是 Hermitian idx = np.argsort(lam)[..., ::-1] # 反转为降序 lam = np.take_along_axis(lam, idx, axis=-1) u = np.take_along_axis(u, idx[..., None, :], axis=-1) return np.clip(lam, 0.0, None), u

build_t3里的looks参数就是等效视数 L 的平方根，传入 3 表示 3×3 的块平均，得到 L=9。np.linalg.eigh比eig更合适，因为它强制按 Hermitian 处理，数值上不会跑出复数特征值。eigen_decompose里做降序重排是必须的：H/α 的定义依赖 λ₁ 最大、λ₃ 最小，顺序错了熵值不受影响，但各向异性 A 会变成负数。

调试时有个快速自检：对单视数据做特征值分解，如果非零特征值超过一个，说明输入数据不是单视的，或者某个环节已经混入了平均；如果多视后三个特征值全是零，那基本可以确定数据里存在全零块或已损坏的像素。

3. 极化特征的三层谱系：功率、相位与统计分布

极化特征不是一堆并列的公式，它有三层清晰的递进关系：第一层是功率和相位这类直接可测的量，第二层是从特征值谱里导出的分布形状量，第三层是带物理假设的模型分解量。层次越往下，物理可解释性越强，但对假设的依赖也越重。搞不清这三层的边界，就会把模型假设的产物当成观测事实来用。

3.1 功率与相位：Span、HH-VV 相位差与复相关系数

最基础的一组特征不需要任何分解。Span（总散射功率）等于相干矩阵的迹，也等于 |Shh|² + 2|Shv|² + |Svv|²，它对地物整体的后向散射强度敏感，是所有特征里最稳定的一个，也是最容易被其他特征“淹没”的一个。

比功率更有价值的是相位类特征。HH 与 VV 的相位差 φ = arg(⟨Shh·Svv*⟩) 直接区分奇次散射与偶次散射：水面、裸地通常在 0° 附近，建筑二面角在 ±180° 附近，植被则因为体散射的随机性而分布弥散。同极化复相关系数 ρ = ⟨Shh·Svv*⟩ / √(⟨|Shh|²⟩⟨|Svv|²⟩) 的模值则刻画了这两种散射的相干程度，植被区通常显著低于建筑区。

这组特征有个天然弱点：对目标的取向角极为敏感。绕视线方向旋转一个角度，功率和相关系数模值不变，但相位差会整体平移。同一栋楼，成像几何差一点，相位差就可能从 +170° 变成 −170°，分类器看到的是两个完全不同的数。常见的做法是先做取向角补偿（deorientation），把相位旋转对齐，再提取相位特征。

3.2 特征值谱特征：熵 H、平均散射角 α、各向异性 A

把特征值归一化成概率 pᵢ = λᵢ / (λ₁+λ₂+λ₃)，就得到了三个描述“散射随机程度”的量：

α 角之所以稳健，是因为它只依赖特征向量的方向，不依赖特征值的绝对大小，乘性噪声和整体增益误差基本不影响它。A 则相反，分母是两个小量之差，在低熵区域（λ₂、λ₃ 都接近 0）几乎是纯噪声，工程上通常要加一个与 λ₁ 成比例的门限，或者干脆在 H < 0.5 的区域不采用 A。

熵 H 的“系统性偏高”值得多说一句。它是多视不足时的典型表现：样本数不够，协方差矩阵的估计误差会让谱“变平”，本来应该集中的能量被摊到三个特征值上，算出来的 H 就偏高。所以 H 的绝对值在不同数据源之间几乎不可比，只有同一景、同一多视配置下的相对大小才有意义。

3.3 模型分解：Freeman-Durden 三成分与它必须满足的假设

模型分解走的是另一条路：先假设观测到的协方差矩阵是若干种理想散射机制的加权和，再反解权重。三成分分解把 ⟨C₃⟩ 写成表面散射、二面角散射、体散射三部分之和，表面散射用 Bragg 模型，二面角用理想导体二面角，体散射用一组随机取向的偶极子。

这条路能给出物理意义极强的分量功率，但它成立的前提是反射对称性假设：⟨Shh·Shv*⟩ ≈ 0 且 ⟨Shv·Svv*⟩ ≈ 0。城市区域里，建筑朝向杂乱，交叉极化与同极化之间的相关性并不为零，这个假设被破坏，反解出来的负功率就是最直接的证据。四成分分解加入了螺旋体散射来吸收一部分不满足对称性的能量，缓解但不根治。

import numpy as np def h_alpha_a(t3, eps=1e-10): """由 T3 计算熵 H、平均散射角 alpha、各向异性 A""" lam, u = np.linalg.eigh(t3) lam = np.clip(lam[::-1], 0.0, None) # 降序 u = u[:, ::-1] total = lam.sum() + eps p = np.clip(lam / total, eps, 1.0) # 避免 log(0) h = -np.sum(p * np.log(p)) / np.log(3.0) # 用 log3 归一化到 [0,1] # alpha_i = arccos(|u_i1|)，u_i1 是第 i 个特征向量的第一个分量 alpha = np.sum(p * np.degrees(np.arccos(np.abs(u[0, :])))) a = (lam[1] - lam[2]) / (lam[1] + lam[2] + eps) return h, alpha, a def freeman_durden(c3, eps=1e-10): """三成分分解，输入 3x3 协方差矩阵（Lexicographic 基）""" c11, c22, c33 = c3[0, 0].real, c3[1, 1].real, c3[2, 2].real c12, c13, c23 = c3[0, 1], c3[0, 2], c3[1, 2] # 体散射分量：由交叉极化功率反推，系数 4/3 来自随机偶极子模型 fv = 4.0 * c22 # 扣除体散射后剩余的表面/二面角功率 s_rem = c11 - fv / 4.0 d_rem = c33 - fv / 4.0 # 二面角与表面散射的相位关系：实部为正取表面，为负取二面角 if (c13.real > 0): fs = s_rem fd = 0.0 else: fd = d_rem fs = 0.0 return np.clip(np.array([fs, fd, fv]), 0.0, None) + eps

h_alpha_a里用np.log(3.0)而不是自然对数，是为了把 H 严格归一化到 0~1，方便和文献里的 H/α 平面分区阈值（0.5 和 0.9）直接对照。如果用了 ln，同一个像素算出来的 H 会系统性偏大，落进高熵区，分区就全错了。

freeman_durden只是一个够用的简化实现，c13.real的符号判断是区分表面散射与二面角散射的关键。真实数据里，这个符号在城区边界处会来回跳变，导致 fs 和 fd 在相邻像素间剧烈翻转，也就是常说的“负功率”和“分解不稳定”。工程上通常加一步功率非负约束，或者干脆改用基于特征值的分解来绕开这个问题。

4. 手写特征提取流水线：从数据读入到特征堆叠

把前面这些特征串成一条可复现的流水线，中间有几处容易出问题的地方，而且出问题时程序不会报错，只会安静地产出一张看起来很正常的错误结果。这类“静默失败”比崩溃难查得多。

4.1 定标与数据格式：为什么必须先做辐射定标

原始散射矩阵的四个通道各有独立的增益和相位偏置，通道间还存在交叉极化泄漏（cross-talk）。如果不做极化定标就直接算 HH-VV 相位差，得到的相位里混着系统相位偏差，水面可能算出 30° 而不是 0°。

常见的处理顺序是：先做辐射定标把 DN 值转成后向散射系数，再做极化定标（用已知的角反射器或者分布目标法）校正通道不平衡和交叉极化。这两步顺序不能颠倒，因为交叉极化泄漏的量级和绝对功率相关，先归一化再校正会更稳定。

还有一个细节是数据的方向。不同产品对复数数据的存放顺序不一致（有的按 [real, imag] 交错，有的分成两个实数组），读错顺序不会报错，但会让所有相位特征变成噪声。验证方法很土但有效：找一块平静水面，看 HH-VV 相位差的均值是否接近 0；如果结果在 0~360° 上均匀分布，那就先回去检查数据读取，别急着调算法。

4.2 滤波与等效视数：斑点噪声如何污染 H/α

斑点噪声不是加性噪声，它是相干成像的乘性干涉结果，单视数据的强度服从指数分布。对复数通道直接做均值滤波，会把相位平均掉，等于人为制造低相干，ρ 会被系统性地压低。

正确的做法是在协方差矩阵或相干矩阵域滤波：先构造 T3，再对 T3 的每个元素做空间平均。精极化 Lee 滤波在此基础上加了边缘对齐的方向窗和同质区选择，能在压噪的同时保住点目标和边缘。工程上如果只需要快速验证，boxcar 平均在 T3 域上做也够用，但窗口别开太大。

等效视数 L 会同时影响多个特征，而且方向一致：L 越小，H 越高、A 越抖、ρ 越低。这意味着用一景数据调出来的 H/α 分区阈值，换一景不同多视配置的数据就不能直接复用。比较稳妥的做法是记录每景数据的 L 值，作为特征的一部分或者至少作为元数据保存下来。

4.3 特征归一化与堆叠：维度上去了，判别力未必上去

把 Span、H、α、A、三个分解分量、相位差、相关系数全部堆起来，很容易到十几二十维。这些特征的量纲和分布差异极大：Span 可以跨好几个数量级，相位差是环形量，α 是 0~90 度的角度，熵是 0~1。直接丢进 SVM 或随机森林，Span 会独占分裂增益，其余特征几乎不起作用。

处理方式因模型而异。树模型对单调变换不敏感，真正需要处理的是 Span 的长尾，通常取对数；SVM 和神经网络则需要标准化，但标准化的统计量必须在训练集上算，再应用到验证集，否则会引入信息泄漏。相位差这类环形量不能简单标准化，要先转成 sin/cos 两个分量，或者先做取向角补偿把它拉回 0 附近。

特征筛选上，一个实用经验是：先把 20 维特征做主成分分析，看前三个主成分解释了多少方差。如果超过 90%，说明特征之间存在大量冗余，此时降维后的判别力往往不降反升；如果前三个主成分只能解释 50%，那说明特征之间信息互补，但也要警惕其中混入了噪声主导的分量。

5. 常见问题排查：H/α 分类图为什么一片糊

5.1 现象：整景熵值偏高，几乎所有像素都落在高熵区

现象很直观：算出来的 H 中位数在 0.8 以上，H/α 平面分区图几乎全部塌到高熵那三格里，本该是水面的区域也显示为高熵。

原因是多视不足加上协方差估计误差。样本数 L 太小时，T3 的谱被噪声抹平，λ₂、λ₃ 被抬高，归一化后的 p 分布趋于均匀，H 自然趋近 1。另一个常见诱因是滤波做在了强度域而不是 T3 域，相位被破坏，同样表现为谱变平。

解决方法是先确认数据的实际等效视数，把 L 提到 9 以上再重算；如果数据本身是多视产品，就不要重复多视；滤波统一切到 T3 域，或者在 T3 域上做精极化 Lee。改完之后先用已知地物（水面、裸地）验证 H 是否落回 0.2 以下，再去看整景。

5.2 现象：训练集精度 95%，换一景数据分类结果全崩

现象是典型的过拟合外观：模型在训练景上表现完美，换一景后精度掉到随机猜测附近。

原因通常不是模型复杂度，而是特征分布漂移。不同景的入射角、多视配置、定标精度都不同，而 H、A、ρ 这些特征对这三者都很敏感。训练集里模型学到的可能是“这一景的熵值整体偏高”，而不是“高熵对应植被”。

解决方法是把与成像几何强相关的特征做归一化或去除，优先保留 α 这种对增益不敏感的量；训练集和验证集必须来自不同景、不同时相，不能在同一景上随机划分；如果条件允许，在特征里显式加入入射角作为协变量，让模型自己学补偿。

5.3 现象：相位差图上出现规则条纹，沿方位向周期性抖动

现象是 HH-VV 相位差图里出现等间距的条纹，条纹方向与航向垂直，周期和脉冲重复频率相关。

原因多半出在数据读取的块边界上，或者是相位解缠的参考点选在了变化剧烈的区域。还有一种情况是多视窗口的步长与图像的行列宽不整除，导致最后一块和前面块的处理方式不一致，边界处出现跳变。

解决方法是先检查数组切片的边界处理，把不整除的余数单独处理或者直接裁掉；相位差不要做解缠，它是环形量，直接取复数的辐角即可；如果是块处理的产物，把窗口改成滑动窗口并保证边缘 padding 的一致性。

5.4 现象：特征从 8 维堆到 40 维，验证精度反而下降

现象是加特征时训练精度一路涨，验证精度在中途达到峰值后开始下滑，继续加特征甚至不如只用 Span + H + α 三维。

原因是特征之间的强共线性把距离度量搞坏了，也有一部分是噪声特征稀释了有效信息。极化特征里重复度极高：Span 和三个分解分量之和高度相关，H 和 A 都来自同一组特征值，相位差和相关系数实部也部分重叠。

解决方法是先做相关性筛查，把相关系数超过 0.95 的特征对删掉一个；再用包裹式方法（例如递归特征消除）在验证集上评估子集，而不是凭直觉堆。实践里 8~12 维往往就是拐点，超过这个数量带来的边际收益，基本被噪声抵消了。

6. 设计取舍：为什么不存在“通用最优特征集”

极化特征的选择本质上是在三种东西之间做权衡：物理可解释性、对成像条件的鲁棒性、以及对具体任务的判别力。这三者很少能同时最优，而且权衡点随任务变化。做地物分类，α 角加熵的两维组合常常已经够用，因为它们刻画的是散射机制的类型，与任务目标直接对齐；做建筑提取，HH-VV 相位差和二次散射分量的贡献远大于熵，因为判别依据是几何结构而不是随机程度；做农作物分类，时序上的相位和相干性变化比任何单景特征都更有信息量。

一个绕不开的现实是：基于模型分解的特征物理意义强，但对反射对称性假设的依赖也强；基于特征值谱的特征假设少、鲁棒性好，但可解释性弱一些，熵和 α 都是统计量而非物理量；纯数据驱动的深度特征判别力可能最强，但需要大量标注，而且换一景就重新训练的成本很高。选哪条路，取决于你有多少标注、数据的一致性有多好、以及结果是否需要向非技术人员解释。

验证方法上，我倾向于先做一个最小基线：只用 Span、H、α 三维跑一遍分类，把这个结果作为下限。任何新增特征都必须在这个基线上有可复现的提升，而且提升要在不同景、不同时相上都成立，不能只在单景的验证集上好看。如果一组特征在跨景验证时提升消失，那它捕捉到的就是这一景的成像条件，不是地物本身。

还有个容易被忽略的取舍是多视窗口的大小。窗口开大，H 和 A 的估计更稳，但小目标会被抹掉；窗口开小，边缘保住了，但低熵区的 A 基本是噪声。我的做法是按任务定窗口：目标检测用小窗口加单特征，地物分类用大窗口加多维特征，同一份数据出两套特征，分别服务两个任务，而不是强行用一套特征打通所有场景。

这些经验基本都是踩出来的。早些年我把一景数据上调好的 20 维特征直接搬到另一景，看到精度从 0.9 掉到 0.4，还以为是代码有 bug，查了两天才发现是两景的等效视数差了四倍，H 的分布整体偏移。从那以后我每次提特征之前，都强制先跑一遍单景自检：确认定标状态、确认 L 值、确认水面相位差和熵值落在合理区间，三样不对就不往下走。这套资料里把矩阵构造、特征值分解和分解算法的实现细节都摊开写了，对着代码走一遍，比只看公式更容易发现自己在哪一步理解偏了。希望帮到你。

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