diff --git a/src/core/modeling/modeling_batch.py b/src/core/modeling/modeling_batch.py index 83f4eda..c61737d 100644 --- a/src/core/modeling/modeling_batch.py +++ b/src/core/modeling/modeling_batch.py @@ -774,7 +774,7 @@ class WaterQualityModelingBatch: _deploy = { 'mnf_mean': np.asarray(_mnf.mean_).ravel(), 'mnf_W': np.asarray(_mnf.W_mnf_), - 'mnf_n_components': int(_mnf.n_components), + 'mnf_n_components': int(getattr(_mnf, 'n_components_', _mnf.n_components)), } # SVR 部署矩阵（rbf kernel 需要 support_vectors_） _deploy['svr_dual_coef'] = np.asarray(_svr.dual_coef_) diff --git a/src/preprocessing/spectral_Preprocessing.py b/src/preprocessing/spectral_Preprocessing.py index 8637b4e..e0624c7 100644 --- a/src/preprocessing/spectral_Preprocessing.py +++ b/src/preprocessing/spectral_Preprocessing.py @@ -369,11 +369,11 @@ class MNFTransformer(TransformerMixin, BaseEstimator, _ArrayAsFloat64): Parameters ---------- - n_components : int, default=10 - 保留的 MNF 成分数量（降维后的特征维度） + n_components : int or float, default=0.95 + int 时保留固定数量的成分；float(0~1) 时根据累积方差比自动确定 """ - def __init__(self, n_components: int = 10): + def __init__(self, n_components=0.95): self.n_components = n_components def fit(self, X, y=None): @@ -382,18 +382,16 @@ class MNFTransformer(TransformerMixin, BaseEstimator, _ArrayAsFloat64): n_samples, n_features = X.shape # 1) 噪声矩阵：相邻行差分 - noise = X[1:] - X[:-1] # (n-1, p) - Cn = np.cov(noise, rowvar=False) # (p, p) - # 正则化：按总方差的 1% 加对角项，强于固定 1e-9 + noise = X[1:] - X[:-1] + Cn = np.cov(noise, rowvar=False) _reg = np.trace(Cn) / n_features * 0.01 Cn += np.eye(n_features) * max(_reg, 1e-6) - # 2) 噪声白化：仅对显著特征值做白化，忽略噪声主导的成分 + # 2) 噪声白化 eigvals_n, eigvecs_n = np.linalg.eigh(Cn) - # 截断：特征值 < max * 1e-4 的视为噪声，不参与白化（置单位矩阵） _keep_mask = eigvals_n > eigvals_n.max() * 1e-4 if _keep_mask.sum() < 2: - _keep_mask[:2] = True # 至少保留 2 个成分 + _keep_mask[:2] = True _kept_vals = np.where(_keep_mask, eigvals_n, 1.0) inv_sqrt = np.where(_keep_mask, 1.0 / np.sqrt(np.maximum(_kept_vals, 1e-12)), 1.0) Wn = eigvecs_n @ np.diag(inv_sqrt) @ eigvecs_n.T @@ -405,9 +403,21 @@ class MNFTransformer(TransformerMixin, BaseEstimator, _ArrayAsFloat64): X_w = (X - self.mean_) @ Wn Cw = np.cov(X_w, rowvar=False) eigvals_w, eigvecs_w = np.linalg.eigh(Cw) - Wp = eigvecs_w[:, ::-1] + # 降序排列 + _order = np.argsort(eigvals_w)[::-1] + eigvals_w = eigvals_w[_order] + eigvecs_w = eigvecs_w[:, _order] + Wp = eigvecs_w - # 5) 最终变换矩阵 + # 5) 自动确定成分数 + if isinstance(self.n_components, float) and 0 < self.n_components <= 1: + _cumvar = np.cumsum(np.maximum(eigvals_w, 0)) / np.sum(np.maximum(eigvals_w, 0)) + self.n_components_ = int(np.searchsorted(_cumvar, self.n_components)) + 1 + self.n_components_ = max(2, min(self.n_components_, n_features)) + else: + self.n_components_ = int(self.n_components) + + # 6) 最终变换矩阵 self.W_mnf_ = Wn @ Wp return self @@ -415,8 +425,7 @@ class MNFTransformer(TransformerMixin, BaseEstimator, _ArrayAsFloat64): def transform(self, X): X = self._to_ndarray(X) X = np.nan_to_num(X, nan=0.0, posinf=0.0, neginf=0.0) - result = (X - self.mean_) @ self.W_mnf_[:, :self.n_components] - # 防 SVR 拒绝 inf：clip 到 float32 安全范围 + result = (X - self.mean_) @ self.W_mnf_[:, :self.n_components_] result = np.nan_to_num(result, nan=0.0, posinf=0.0, neginf=0.0) result = np.clip(result, -1e15, 1e15) return result @@ -427,39 +436,17 @@ class MNFTransformer(TransformerMixin, BaseEstimator, _ArrayAsFloat64): # ============================================================================ class PhysicalFeatureExtractor(TransformerMixin, BaseEstimator, _ArrayAsFloat64): - """从原始高光谱波段中提取物理意义明确的比值/指数特征。 + """从原始高光谱波段计算全部水质/植被指数（63 个公式）。 - 输入 X 为 pandas DataFrame 时，列名应为波长字符串（如 "800.000000"）； - 输入为 numpy ndarray 时，需在 fit 中传入 wavelengths 参数。 + 利用项目已有的 WaterQualityIndexCalculator，自动根据 wavelength 列名 + 匹配各公式所需的最近波段，计算涵盖 Chl/BGA/Turb/TSM 的完整指数集。 - 计算的特征包括： - - NDVI: (R800 - R670) / (R800 + R670) - - NDWI: (R550 - R800) / (R550 + R800) - - 蓝藻比值: R705 / R670 - - Parameters - ---------- - wavelengths : list of float or None - ndarray 输入时需显式传入波长列表；DataFrame 输入时自动从列名解析。 + 输入 X 为 DataFrame 时自动解析列名；ndarray 时需在 __init__ 传入 wavelengths。 """ - _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] @@ -467,29 +454,23 @@ class PhysicalFeatureExtractor(TransformerMixin, BaseEstimator, _ArrayAsFloat64) 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)) + # 预加载公式列表（fit 时只做一次） + from src.utils.water_index import WaterQualityIndexCalculator + self._calc = WaterQualityIndexCalculator() + self._formulas = self._calc.list_available() + self._n_features_out_ = len(self._formulas) 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) + # 重建 DataFrame（列名 = 波长），供 WaterQualityIndexCalculator 使用 + col_names = [f'{wl:.6f}' for wl in self.wavelengths] + df = pd.DataFrame(X, columns=col_names) + result_df = self._calc.calculate_many(self._formulas, df, fast=True) + out = np.asarray(result_df, dtype=np.float64) + out = np.nan_to_num(out, nan=0.0, posinf=0.0, neginf=0.0) + out = np.clip(out, -1e15, 1e15) + return out # ============================================================================ @@ -535,7 +516,7 @@ def get_preprocessing_transformer(method: str, wavelengths=None): from sklearn.pipeline import FeatureUnion return FeatureUnion([ ('physical', PhysicalFeatureExtractor(wavelengths=wavelengths)), - ('mnf', MNFTransformer(n_components=10)), + ('mnf', MNFTransformer()), # n_components=0.95 自动确定 ]) if method not in _PREPROCESSING_TRANSFORMERS: print(f"未知预处理方法 '{method}'，回退为 IdentityTransformer")