feat: PhysicalFeature 扩展到全部63个水色指数 + MNF成分数自动确定
PhysicalFeatureExtractor: - 废弃硬编码3指数(NDVI/NDWI/AlgaeRatio)，改为调用 WaterQualityIndexCalculator.calculate_many() 计算全部63个公式 - 涵盖 BGA(25)/Chl(16)/Turb(10)/TSM(3) 等完整指数集 MNFTransformer: - n_components 支持 float(0~1): 根据累积方差比自动确定成分数 n_components=0.95 → 保留95%方差的成分 - fit 后保存 n_components_ (实际成分数) 供部署导出使用 - Wp 显式按特征值降序排列
This commit is contained in:
duxin
@ -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_)
|
||||
|
||||
@ -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")
|
||||
|
||||
Reference in New Issue
Block a user