feat: 新增 MNFTransformer + PhysicalFeatureExtractor + DualStream_MNF

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维)
This commit is contained in:
duxin
2026-07-29 10:02:31 +08:00
parent c4c70b2ff6
commit 7d8d848e58

View File

@ -355,13 +355,140 @@ class WVAETransformer(TransformerMixin, BaseEstimator, _ArrayAsFloat64):
# ============================================================================
# 工厂函数:根据方法名返回对应的 sklearn 兼容 TransformerNone 表示无预处理
# 双流特征提取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()