From 7d8d848e58bf03aeceb5c3e32bd99c87d477c3d1 Mon Sep 17 00:00:00 2001 From: duxin Date: Wed, 29 Jul 2026 10:02:31 +0800 Subject: [PATCH] =?UTF-8?q?feat:=20=E6=96=B0=E5=A2=9E=20MNFTransformer=20+?= =?UTF-8?q?=20PhysicalFeatureExtractor=20+=20DualStream=5FMNF?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit spectral_Preprocessing.py: - MNFTransformer: 纯 numpy/scipy 实现的 MNF 降维 fit: 噪声协方差→白化→PCA 旋转，导出 self.mean_/self.W_mnf_ transform: (X - mean) @ W_mnf_[:, :n_components] - PhysicalFeatureExtractor: 从原始波段提取 NDVI/NDWI/蓝藻比值 支持 DataFrame(列名波长) 和 ndarray(wavelengths参数) 输入 - get_preprocessing_transformer: DualStream_MNF 返回 FeatureUnion 并行连接 PhysicalFeatureExtractor(3维) + MNFTransformer(10维) --- src/preprocessing/spectral_Preprocessing.py | 148 +++++++++++++++++++- 1 file changed, 141 insertions(+), 7 deletions(-) diff --git a/src/preprocessing/spectral_Preprocessing.py b/src/preprocessing/spectral_Preprocessing.py index b7e839f..e7f1381 100644 --- a/src/preprocessing/spectral_Preprocessing.py +++ b/src/preprocessing/spectral_Preprocessing.py @@ -355,13 +355,140 @@ class WVAETransformer(TransformerMixin, BaseEstimator, _ArrayAsFloat64): # ============================================================================ -# 工厂函数：根据方法名返回对应的 sklearn 兼容 Transformer（None 表示无预处理） +# 双流特征提取：MNF 降维（纯 numpy/scipy 实现，可导出矩阵用于 C++/Rust 部署） +# ============================================================================ + +class MNFTransformer(TransformerMixin, BaseEstimator, _ArrayAsFloat64): + """Minimum Noise Fraction (MNF) 变换 — 纯 numpy/scipy 实现。 + + 不依赖 spectral 库，fit 阶段计算白化矩阵 + 旋转矩阵， + transform 阶段进行线性变换降维。 + + 为 C++/Rust 部署设计：训练完成后可从实例提取 self.mean_ 和 + self.W_mnf_，在脱离 Python/ sklearn 的环境完成预测前置处理。 + + Parameters + ---------- + n_components : int, default=10 + 保留的 MNF 成分数量（降维后的特征维度） + """ + + def __init__(self, n_components: int = 10): + self.n_components = n_components + + def fit(self, X, y=None): + X = self._to_ndarray(X) + n_samples, n_features = X.shape + + # 1) 噪声矩阵：相邻行差分（假设输入按空间/时间连续排列） + noise = X[1:] - X[:-1] # (n-1, p) + Cn = np.cov(noise, rowvar=False) # (p, p) + Cn += np.eye(n_features) * 1e-9 # 防止奇异 + + # 2) 噪声白化：Cn = U diag(S) Vt — 对称矩阵用 eigh + eigvals_n, eigvecs_n = np.linalg.eigh(Cn) + # eigh 返回升序排列；白化矩阵 = V @ diag(1/sqrt(s)) @ V.T + inv_sqrt = 1.0 / np.sqrt(np.maximum(eigvals_n, 1e-12)) + Wn = eigvecs_n @ np.diag(inv_sqrt) @ eigvecs_n.T # (p, p) + + # 3) 均值 + self.mean_ = X.mean(axis=0, keepdims=True) # (1, p) + + # 4) 白化后的数据协方差 → PCA + X_w = (X - self.mean_) @ Wn # (n, p) + Cw = np.cov(X_w, rowvar=False) # (p, p) + eigvals_w, eigvecs_w = np.linalg.eigh(Cw) + # 降序排列（eigh 升序，翻转） + Wp = eigvecs_w[:, ::-1] # (p, p) + + # 5) 最终变换矩阵 + self.W_mnf_ = Wn @ Wp # (p, p) + + return self + + def transform(self, X): + X = self._to_ndarray(X) + return (X - self.mean_) @ self.W_mnf_[:, :self.n_components] + + +# ============================================================================ +# 物理特征提取器：从原始波段计算水质/植被指数 +# ============================================================================ + +class PhysicalFeatureExtractor(TransformerMixin, BaseEstimator, _ArrayAsFloat64): + """从原始高光谱波段中提取物理意义明确的比值/指数特征。 + + 输入 X 为 pandas DataFrame 时，列名应为波长字符串（如 "800.000000"）； + 输入为 numpy ndarray 时，需在 fit 中传入 wavelengths 参数。 + + 计算的特征包括： + - NDVI: (R800 - R670) / (R800 + R670) + - NDWI: (R550 - R800) / (R550 + R800) + - 蓝藻比值: R705 / R670 + + Parameters + ---------- + wavelengths : list of float or None + ndarray 输入时需显式传入波长列表；DataFrame 输入时自动从列名解析。 + """ + + _TARGET_NAMES = ('NDVI', 'NDWI', 'AlgaeRatio') + + # 默认目标波长 (nm) — 通过 fit 时扫描实际波长列自动修正为最接近值 + _TARGET_WL = { + 'NDVI': (800.0, 670.0), + 'NDWI': (550.0, 800.0), + 'AlgaeRatio': (705.0, 670.0), + } + + def __init__(self, wavelengths=None): + self.wavelengths = wavelengths + + def _find_nearest_wl(self, wl_array, target): + """返回最接近 target 的波长值""" + idx = np.argmin(np.abs(np.asarray(wl_array) - target)) + return float(wl_array[idx]), int(idx) + + def fit(self, X, y=None): + if isinstance(X, pd.DataFrame): + self.wavelengths = [float(str(c)) for c in X.columns] + elif self.wavelengths is None: + raise ValueError( + "ndarray 输入时必须在 __init__ 中提供 wavelengths 参数" + ) + # 预计算每个特征对应的列索引 + self._feat_cols_ = [] + wl_arr = self.wavelengths + for name, (wl_a, wl_b) in self._TARGET_WL.items(): + _, ia = self._find_nearest_wl(wl_arr, wl_a) + _, ib = self._find_nearest_wl(wl_arr, wl_b) + self._feat_cols_.append((name, ia, ib)) + return self + + def transform(self, X): + X = self._to_ndarray(X) + feats = [] + for name, ia, ib in self._feat_cols_: + a = X[:, ia] + b = X[:, ib] + denom = a + b + denom = np.where(np.abs(denom) < 1e-12, np.sign(denom) * 1e-12, denom) + if name == 'AlgaeRatio': + denom2 = np.where(np.abs(b) < 1e-12, np.sign(b) * 1e-12, b) + feats.append((a / denom2).reshape(-1, 1)) + else: + feats.append(((a - b) / denom).reshape(-1, 1)) + return np.hstack(feats) + + +# ============================================================================ +# 工厂函数：根据方法名返回对应的 sklearn 兼容 Transformer # ============================================================================ _PREPROCESSING_TRANSFORMERS = { 'None': IdentityTransformer, - 'MMS': MinMaxScaler, # sklearn 自带 - 'SS': StandardScaler, # sklearn 自带 + 'MMS': MinMaxScaler, + 'SS': StandardScaler, 'CT': CTTransformer, 'SNV': SNVTransformer, 'MA': MATransformer, @@ -377,19 +504,26 @@ _PREPROCESSING_TRANSFORMERS = { def get_preprocessing_transformer(method: str): """根据预处理方法名返回 sklearn 兼容的 Transformer 实例。 - - method 为 "None" 或 None：返回 IdentityTransformer（等价于无处理） + - method 为 "None" 或 None：返回 IdentityTransformer - method 为 "MMS"/"SS"：直接返回 sklearn 自带 MinMaxScaler/StandardScaler - - method 为 "CT"/"SNV"/"MA"/"SG"/"MSC"/"D1"/"D2"/"DT"/"WVAE"：返回对应包装类 - - method 不识别：返回 IdentityTransformer + 打印警告（与原 Preprocessing 行为一致） + - method 为 "DualStream_MNF"：返回 FeatureUnion，并行连接 + PhysicalFeatureExtractor(3 维物理指数) + MNFTransformer(10 维降维特征) + - method 不识别：返回 IdentityTransformer + 打印警告 Args: - method: 预处理方法名（大小写敏感，与 Preprocessing() 一致） + method: 预处理方法名 Returns: sklearn 兼容的 Transformer 实例（可直接放入 Pipeline） """ if method is None: return IdentityTransformer() + if method == "DualStream_MNF": + from sklearn.pipeline import FeatureUnion + return FeatureUnion([ + ('physical', PhysicalFeatureExtractor()), + ('mnf', MNFTransformer(n_components=10)), + ]) if method not in _PREPROCESSING_TRANSFORMERS: print(f"未知预处理方法 '{method}'，回退为 IdentityTransformer") return IdentityTransformer()