简介面向遥感与雷达信号处理学习者的极化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 先确认三件事定标、基、元素排列在写任何一行特征代码之前值得花十分钟翻元数据。常见做法是核对下面四项缺一项后面就会以“精度莫名其妙掉”的形式还回来。检查项常见表现处理方式幅度量纲数据是 16 位无符号整型直接当线性功率用σ0 差几十 dB乘定标系数开方后才是复幅度定标类型σ0、β0、γ0 混用地形起伏区差异被当成特征明确用的是哪一种全流程统一矩阵基给的是四个极化通道、C3、还是 T3用 U3 变换统一到 T3 再提特征元素排列6 元素上三角实虚交错、还是 9 元素全存按元数据顺序重建别按经验猜第三条最容易被忽略。词典基下的向量是 $[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, reciprocalTrue): hh/hv/vh/vv: 复数二维数组已定标并开方单位是复幅度 if reciprocal: # 互易介质下两个交叉极化理论上相等取平均可以压掉通道间相位偏差 hv 0.5 * (hv vh) # 词典基sqrt(2) 是为了让总功率等于四个通道功率之和 s np.stack([hh, np.sqrt(2.0) * hv, vv], axis0) # (3, H, W) c3 s[:, None, :, :] * np.conj(s[None, :, :, :]) # 外积 - (3, 3, H, W) return c3 # T3 U3 C3 U3^HU3 为 Pauli 变换矩阵 U3 np.array([[1, 0, 1], [1, 0, -1], [0, np.sqrt(2), 0]], dtypenp.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 mscattering_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_{i1}^{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, keepdimsTrue) 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(invalidignore, divideignore): 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, lamnp.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, beta1 的经典假设 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, residhv2 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 比值与指数特征便宜且抗噪的那一档第三档特征是比值和相关系数。它们不需要特征值分解不需要迭代求解对相干斑的敏感度也低得多在特征集里应该占据基础位置。特征表达式抗噪性主要用途后向散射 σ0$10\log_{10}\langleS_{pq}^2\rangle$HH/VV 比$10\log_{10}(\sigma_{HH}/\sigma_{VV})$高表面粗糙度、二面角主导区交叉极化比$10\log_{10}(\sigma_{HV}/\sigma_{VV})$高体散射强度植被生物量代理去极化度$1 -\rho_{HHVV}$RVI$8\sigma_{HV}/(\sigma_{HH}\sigma_{VV}2\sigma_{HV})$高植被区分对定标误差不敏感相关系数相位$\arg\langle S_{HH}S_{VV}^*\rangle$低二面角/表面散射判别需高 ENLdef ratio_features(c3, floor1e-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, win7): 退化版矩阵域均值滤波 半正定投影够用作对比基线 k np.ones((1, 1, win, win), dtypenp.float32) / (win * win) sm np.stack([np.stack([ np.convolve(c3[i, j].ravel(), k.ravel(), modesame).reshape(c3.shape[2:]) for j in range(3)], axis0) for i in range(3)], axis0) 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, floor1e-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_estimators300, min_samples_leaf2, class_weightbalanced_subsample, random_state0, n_jobs-1)), ]) pipe.fit(X_train, y_train)floor1e-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, thr0.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_featuresFalse, random_state0) order np.argsort(mi)[::-1] keep, dropped prune_by_corr(X_train, feat_names, thr0.95) X_sel X_train[:, keep]thr0.95是经验值。0.9 以下会开始误删——H 和交叉极化比在植被区的相关系数能到 0.85但它们提供的信息并不完全重叠。0.98 以上则几乎不删东西。0.95 附近是常见取舍点。mutual_info_classif的discrete_featuresFalse必须显式写。默认值在不同版本间变过如果被当成离散特征处理连续的 dB 值会被分箱互信息估计会严重失真。剪完之后通常还剩十几到二十维。这时候再回头看消融实验就会发现某些看着很“高级”的特征其实贡献为零。哪些能删、怎么验证放到最后一章用数据说。5. 五个常见失败场景从负特征值到精度虚高下面这几条都是能复现的、有明确现象和确定处理方式的失败场景。按出现频率从高到低排。5.1 特征图出现 NaN 或整片黑现象H 或 α 的图上出现成片的零值或 NaN位置集中在山体阴影、水体、以及图像边缘。原因主要有两个。一是总功率为 0 或接近 0 的像元无效值填充、定标本文还有配套的精品资源点击获取 
长河君 · 2026/10/27 15:07:53 · 阅读 1.2k 极化SAR特征提取全链路：从C3/T3构建到特征筛选与验证
简介面向遥感与雷达信号处理学习者的极化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 先确认三件事定标、基、元素排列在写任何一行特征代码之前值得花十分钟翻元数据。常见做法是核对下面四项缺一项后面就会以“精度莫名其妙掉”的形式还回来。检查项常见表现处理方式幅度量纲数据是 16 位无符号整型直接当线性功率用σ0 差几十 dB乘定标系数开方后才是复幅度定标类型σ0、β0、γ0 混用地形起伏区差异被当成特征明确用的是哪一种全流程统一矩阵基给的是四个极化通道、C3、还是 T3用 U3 变换统一到 T3 再提特征元素排列6 元素上三角实虚交错、还是 9 元素全存按元数据顺序重建别按经验猜第三条最容易被忽略。词典基下的向量是 $[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, reciprocalTrue): hh/hv/vh/vv: 复数二维数组已定标并开方单位是复幅度 if reciprocal: # 互易介质下两个交叉极化理论上相等取平均可以压掉通道间相位偏差 hv 0.5 * (hv vh) # 词典基sqrt(2) 是为了让总功率等于四个通道功率之和 s np.stack([hh, np.sqrt(2.0) * hv, vv], axis0) # (3, H, W) c3 s[:, None, :, :] * np.conj(s[None, :, :, :]) # 外积 - (3, 3, H, W) return c3 # T3 U3 C3 U3^HU3 为 Pauli 变换矩阵 U3 np.array([[1, 0, 1], [1, 0, -1], [0, np.sqrt(2), 0]], dtypenp.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 mscattering_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_{i1}^{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, keepdimsTrue) 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(invalidignore, divideignore): 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, lamnp.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, beta1 的经典假设 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, residhv2 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 比值与指数特征便宜且抗噪的那一档第三档特征是比值和相关系数。它们不需要特征值分解不需要迭代求解对相干斑的敏感度也低得多在特征集里应该占据基础位置。特征表达式抗噪性主要用途后向散射 σ0$10\log_{10}\langleS_{pq}^2\rangle$HH/VV 比$10\log_{10}(\sigma_{HH}/\sigma_{VV})$高表面粗糙度、二面角主导区交叉极化比$10\log_{10}(\sigma_{HV}/\sigma_{VV})$高体散射强度植被生物量代理去极化度$1 -\rho_{HHVV}$RVI$8\sigma_{HV}/(\sigma_{HH}\sigma_{VV}2\sigma_{HV})$高植被区分对定标误差不敏感相关系数相位$\arg\langle S_{HH}S_{VV}^*\rangle$低二面角/表面散射判别需高 ENLdef ratio_features(c3, floor1e-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, win7): 退化版矩阵域均值滤波 半正定投影够用作对比基线 k np.ones((1, 1, win, win), dtypenp.float32) / (win * win) sm np.stack([np.stack([ np.convolve(c3[i, j].ravel(), k.ravel(), modesame).reshape(c3.shape[2:]) for j in range(3)], axis0) for i in range(3)], axis0) 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, floor1e-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_estimators300, min_samples_leaf2, class_weightbalanced_subsample, random_state0, n_jobs-1)), ]) pipe.fit(X_train, y_train)floor1e-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, thr0.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_featuresFalse, random_state0) order np.argsort(mi)[::-1] keep, dropped prune_by_corr(X_train, feat_names, thr0.95) X_sel X_train[:, keep]thr0.95是经验值。0.9 以下会开始误删——H 和交叉极化比在植被区的相关系数能到 0.85但它们提供的信息并不完全重叠。0.98 以上则几乎不删东西。0.95 附近是常见取舍点。mutual_info_classif的discrete_featuresFalse必须显式写。默认值在不同版本间变过如果被当成离散特征处理连续的 dB 值会被分箱互信息估计会严重失真。剪完之后通常还剩十几到二十维。这时候再回头看消融实验就会发现某些看着很“高级”的特征其实贡献为零。哪些能删、怎么验证放到最后一章用数据说。5. 五个常见失败场景从负特征值到精度虚高下面这几条都是能复现的、有明确现象和确定处理方式的失败场景。按出现频率从高到低排。5.1 特征图出现 NaN 或整片黑现象H 或 α 的图上出现成片的零值或 NaN位置集中在山体阴影、水体、以及图像边缘。原因主要有两个。一是总功率为 0 或接近 0 的像元无效值填充、定标本文还有配套的精品资源点击获取