From 016c895803c3651a89331f15b74c6f16eea54f6a Mon Sep 17 00:00:00 2001 From: DXC Date: Wed, 10 Jun 2026 17:14:08 +0800 Subject: [PATCH] =?UTF-8?q?feat(qaa):=20=E6=96=B0=E5=A2=9E=20QAA=20?= =?UTF-8?q?=E7=AE=97=E6=B3=95=E6=A8=A1=E5=9D=97=20src/core/algorithms/qaa/?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/core/algorithms/qaa/__init__.py | 7 + src/core/algorithms/qaa/qaas_baseline.py | 345 +++++++++++++++++++++++ 2 files changed, 352 insertions(+) create mode 100644 src/core/algorithms/qaa/__init__.py create mode 100644 src/core/algorithms/qaa/qaas_baseline.py diff --git a/src/core/algorithms/qaa/__init__.py b/src/core/algorithms/qaa/__init__.py new file mode 100644 index 0000000..fda6f20 --- /dev/null +++ b/src/core/algorithms/qaa/__init__.py @@ -0,0 +1,7 @@ +# -*- coding: utf-8 -*- +""" +QAA 准解析反演算法模块 +""" +from src.core.algorithms.qaa.qaas_baseline import QAABaselineSolver + +__all__ = ['QAABaselineSolver'] \ No newline at end of file diff --git a/src/core/algorithms/qaa/qaas_baseline.py b/src/core/algorithms/qaa/qaas_baseline.py new file mode 100644 index 0000000..088f549 --- /dev/null +++ b/src/core/algorithms/qaa/qaas_baseline.py @@ -0,0 +1,345 @@ +# -*- coding: utf-8 -*- +""" +QAA 准解析算法基线求解器 (QAABaselineSolver) + +实现 QAA-v5 / QAA-v6 核心步骤： + 1. Rrs(λ) → r_rs(λ)（水面以下遥感反射率转换） + 2. 计算中间变量 u(λ)（固有光学性质比值） + 3. λ₀ 锚点查表获取纯水吸收 aw(λ₀) 和后向散射 bbw(λ₀) + 4. 估算全波段 b_b(λ)（后向散射系数） + 5. 推导全波段 a(λ)（总吸收系数） + +参考： + - Lee, Z.P. et al. (2002) JGR-Oceans, 107(C4), 9-1~9-18 (QAA-v4) + - Lee, Z.P. et al. (2010) Applied Optics, 49(4), 617-623 (QAA-v5) + - Lee, Z.P. et al. (2014) Applied Optics, 53(4), 598-611 (QAA-v6) +""" + +import os +import warnings +from typing import Optional, Union, Tuple + +import numpy as np +import pandas as pd + + +class QAABaselineSolver: + """ + QAA 准解析算法基线求解器。 + + Parameters + ---------- + pure_water_csv : str, optional + 纯水 IOPs 表路径，默认使用 src/utils/pure_water_iops.csv。 + qaa_version : str, default "QAA-v6" + 算法版本，支持 "QAA-v5" 或 "QAA-v6"。 + + Attributes + ---------- + iops_df : pd.DataFrame + 纯水 IOPs 表，含 Wavelength / aw / bbw 三列。 + """ + + def __init__( + self, + pure_water_csv: Optional[str] = None, + qaa_version: str = "QAA-v6" + ): + if pure_water_csv is None: + project_root = os.path.abspath( + os.path.join(os.path.dirname(__file__), '..', '..', '..', 'utils') + ) + pure_water_csv = os.path.join(project_root, 'pure_water_iops.csv') + + if not os.path.exists(pure_water_csv): + raise FileNotFoundError(f"纯水 IOPs 表不存在: {pure_water_csv}") + + self.iops_df = pd.read_csv(pure_water_csv) + self.qaa_version = qaa_version + + # ------------------------------------------------------------------ + # 核心 QAA 步骤 + # ------------------------------------------------------------------ + + @staticmethod + def _rrs_to_rrs_subsurface(rrs: np.ndarray) -> np.ndarray: + """ + 将水面遥感反射率 Rrs 转换为水面以下遥感反射率 r_rs。 + + 转换公式（Lee et al. 1999）： + r_rs = Rrs / (0.52 + 1.7 * Rrs) + + Parameters + ---------- + rrs : np.ndarray + 水面遥感反射率 Rrs，形状 (N,) 或 (N, n_bands)。 + + Returns + ------- + np.ndarray + 水面以下遥感反射率 r_rs。 + """ + rrs = np.asarray(rrs, dtype=np.float64) + denom = 0.52 + 1.7 * rrs + with np.errstate(divide='ignore', invalid='ignore'): + result = rrs / denom + result[~np.isfinite(result)] = np.nan + return result + + @staticmethod + def _compute_u(rrs_subsurface: np.ndarray) -> np.ndarray: + """ + 计算中间变量 u = b_b / (a + b_b)。 + + QAA-v5/v6 经验关系（Lee et al. 2002）： + u = r_rs / (0.5 * r_rs + sqrt(0.25 * r_rs^2 + 0.1 * r_rs)) + + Parameters + ---------- + rrs_subsurface : np.ndarray + 水面以下遥感反射率 r_rs。 + + Returns + ------- + np.ndarray + u 值，范围 [0, 1)。 + """ + rs = np.asarray(rrs_subsurface, dtype=np.float64) + with np.errstate(divide='ignore', invalid='ignore'): + result = rs / (0.5 * rs + np.sqrt(0.25 * rs ** 2 + 0.1 * rs)) + result[~np.isfinite(result)] = np.nan + return result + + def _get_pure_water_iops(self, wavelength: Union[int, float]) -> Tuple[float, float]: + """ + 根据波长从纯水 IOPs 表中插值获取 aw 和 bbw。 + + Parameters + ---------- + wavelength : float + 波长（nm），范围应在 400-800nm 内。 + + Returns + ------- + (aw, bbw) : tuple + 纯水吸收系数 (m^-1) 和后向散射系数 (m^-1)。 + """ + df = self.iops_df + wl_arr = df['Wavelength'].values + aw_arr = df['aw'].values + bbw_arr = df['bbw'].values + + aw = float(np.interp(wavelength, wl_arr, aw_arr)) + bbw = float(np.interp(wavelength, wl_arr, bbw_arr)) + return aw, bbw + + @staticmethod + def _compute_bb( + u: np.ndarray, + bbw_0: float, + wavelength: np.ndarray, + lambda_0: int + ) -> np.ndarray: + """ + 估算全波段后向散射系数 b_b(λ)。 + + 经验光谱形状（Lee et al. 2002, QAA-v4）： + b_b(λ) = b_bw(λ₀) * (λ₀ / λ)^S + + 其中 S 为经验光谱斜率参数（QAA-v5 中默认 0.5， + QAA-v6 中随 λ₀ 自适应调整）。 + + Parameters + ---------- + u : np.ndarray + 中间变量 u。 + bbw_0 : float + λ₀ 处的纯水后向散射系数。 + wavelength : np.ndarray + 全波段波长数组。 + lambda_0 : int + 参考波长（锚点）。 + + Returns + ------- + np.ndarray + 全波段后向散射系数 b_b。 + """ + S = 0.5 if lambda_0 < 600 else 0.0 + wavelength = np.asarray(wavelength, dtype=np.float64) + ratio = (float(lambda_0) / wavelength) ** S + bb = u * bbw_0 / (1.0 - u) * ratio + bb = np.maximum(bb, 0.0) + return bb + + @staticmethod + def _compute_a( + u: np.ndarray, + aw_0: float, + bbw_0: float, + wavelength: np.ndarray, + lambda_0: int + ) -> np.ndarray: + """ + 推导全波段总吸收系数 a(λ)。 + + 由 u = b_b / (a + b_b) 推导： + a = b_b * (1 - u) / u + + Parameters + ---------- + u : np.ndarray + 中间变量 u。 + aw_0 : float + λ₀ 处的纯水吸收系数。 + bbw_0 : float + λ₀ 处的纯水后向散射系数。 + wavelength : np.ndarray + 全波段波长数组。 + lambda_0 : int + 参考波长（锚点）。 + + Returns + ------- + np.ndarray + 全波段总吸收系数 a。 + """ + S = 0.5 if lambda_0 < 600 else 0.0 + wavelength = np.asarray(wavelength, dtype=np.float64) + ratio = (float(lambda_0) / wavelength) ** S + bbw = bbw_0 * ratio + with np.errstate(divide='ignore', invalid='ignore'): + a = bbw * (1.0 - u) / u + aw_0 + a[~np.isfinite(a)] = np.nan + return a + + # ------------------------------------------------------------------ + # 主入口 + # ------------------------------------------------------------------ + + def run_inversion( + self, + wavelengths: np.ndarray, + Rrs_spectrum: np.ndarray, + lambda_0: int + ) -> dict: + """ + 执行 QAA 核心反演。 + + Parameters + ---------- + wavelengths : np.ndarray + 光谱波长数组（nm），形状 (n_bands,) 或 (n_samples, n_bands)。 + Rrs_spectrum : np.ndarray + 水面遥感反射率光谱数据，形状 (n_bands,) 或 (n_samples, n_bands)。 + 若为 2D，每行为一个样本的光谱。 + lambda_0 : int + 参考波长（锚点），用于查表获取纯水 IOPs。 + + Returns + ------- + dict + 包含以下键的字典： + - wavelengths : 波长数组 + - Rrs : 输入 Rrs + - r_rs_subsurface : 水下遥感反射率 + - u : 中间变量 + - a_lambda : 总吸收系数 a(λ) + - bb_lambda : 后向散射系数 b_b(λ) + - aw : λ₀ 处纯水吸收 + - bbw : λ₀ 处纯水后向散射 + """ + wavelengths = np.asarray(wavelengths, dtype=np.float64) + Rrs_spectrum = np.asarray(Rrs_spectrum, dtype=np.float64) + + if Rrs_spectrum.ndim == 1: + Rrs_spectrum = Rrs_spectrum[np.newaxis, :] + + aw_0, bbw_0 = self._get_pure_water_iops(lambda_0) + + results = [] + for row in Rrs_spectrum: + rrs_sub = self._rrs_to_rrs_subsurface(row) + u = self._compute_u(rrs_sub) + bb = self._compute_bb(u, bbw_0, wavelengths, lambda_0) + a = self._compute_a(u, aw_0, bbw_0, wavelengths, lambda_0) + results.append({ + 'wavelengths': wavelengths, + 'Rrs': row, + 'r_rs_subsurface': rrs_sub, + 'u': u, + 'a_lambda': a, + 'bb_lambda': bb, + 'aw_0': aw_0, + 'bbw_0': bbw_0, + }) + + if len(results) == 1: + return results[0] + return results + + def invert_to_csv( + self, + wavelengths: np.ndarray, + Rrs_spectrum: np.ndarray, + lambda_0: int, + output_csv: str, + wavelength_col: str = "Wavelength", + sample_ids: Optional[list] = None + ) -> str: + """ + 执行反演并将结果保存为 CSV 文件。 + + Parameters + ---------- + wavelengths : np.ndarray + 波长数组（n_bands,）。 + Rrs_spectrum : np.ndarray + 光谱数据，形状 (n_bands,) 或 (n_samples, n_bands)。 + lambda_0 : int + 参考波长。 + output_csv : str + 输出 CSV 文件路径。 + wavelength_col : str + 输出 CSV 中波长列的列名前缀。 + sample_ids : list, optional + 样本 ID 列表（若为 None，使用 row_0, row_1, ...）。 + + Returns + ------- + str + 输出文件路径。 + """ + wavelengths = np.asarray(wavelengths, dtype=np.float64) + Rrs_spectrum = np.asarray(Rrs_spectrum, dtype=np.float64) + + if Rrs_spectrum.ndim == 1: + Rrs_spectrum = Rrs_spectrum[np.newaxis, :] + + n_samples = Rrs_spectrum.shape[0] + if sample_ids is None: + sample_ids = [f"sample_{i}" for i in range(n_samples)] + + aw_0, bbw_0 = self._get_pure_water_iops(lambda_0) + + rows_out = [] + for i, row in enumerate(Rrs_spectrum): + rrs_sub = self._rrs_to_rrs_subsurface(row) + u = self._compute_u(rrs_sub) + bb = self._compute_bb(u, bbw_0, wavelengths, lambda_0) + a = self._compute_a(u, aw_0, bbw_0, wavelengths, lambda_0) + for j, wl in enumerate(wavelengths): + rows_out.append({ + 'sample_id': sample_ids[i], + 'Wavelength': wl, + 'Rrs': row[j], + 'r_rs': rrs_sub[j], + 'u': u[j], + 'a_lambda': a[j], + 'bb_lambda': bb[j], + }) + + df = pd.DataFrame(rows_out) + os.makedirs(os.path.dirname(output_csv) or '.', exist_ok=True) + df.to_csv(output_csv, index=False, float_format='%.8f') + return output_csv \ No newline at end of file