feat(qaa): 新增 QAA 算法模块 src/core/algorithms/qaa/

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']

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