Add validated land clutter models

Author: SpaceEngineerSS

docs/MODEL_FIDELITY.md

@@ -45,7 +45,8 @@ runtime model:
 - PRF, pulse width, receiver noise bandwidth, noise figure, system temperature and losses;
 - probability of false alarm and number of integrated pulses;
 - atmospheric temperature, pressure and water-vapour density;
-- clutter enablement, terrain class, sea state and rain rate;
+- clutter enablement, terrain class, land model and its calibrated surface parameters, sea
+  state and rain rate;
 - target kinematics, mean RCS, Swerling case and onboard jammer power.
 
 Values omitted by a scenario use documented defaults. A saved run must record both supplied and
@@ -83,15 +84,22 @@ validation, input pedigree, uncertainty characterization, and result robustness.
 - V. Gregers-Hansen and R. Mital, "An Improved Empirical Model for Radar Sea Clutter
   Reflectivity," IEEE Transactions on Aerospace and Electronic Systems, 2012,
   DOI: 10.1109/TAES.2012.6324732.
+- Y. Oh, K. Sarabandi, and F. T. Ulaby, "An Empirical Model and an Inversion Technique for
+  Radar Scattering from Bare Soil Surfaces," IEEE Transactions on Geoscience and Remote
+  Sensing, 1992, DOI: 10.1109/36.134086.
+- J. B. Billingsley, *Low-Angle Radar Land Clutter: Measurements and Empirical Models*,
+  William Andrew Publishing, 2002.
 
 ## Known baseline limitations
 
 - Gaseous attenuation uses the P.676-13 Annex 1 line-by-line equations for a homogeneous
   terrestrial path. Slant paths through height-dependent atmospheric profiles are not yet
   represented.
-- Land-clutter terrain-class coefficients remain an empirical screening model and require
-  site-specific measurement or calibrated soil/vegetation inputs for quantitative use. Sea clutter
-  uses the NRL model within its stated frequency, grazing-angle, and sea-state domain.
+- The constant-gamma land model is a screening model. Terrain-class gamma values are nominal
+  priors; quantitative work requires measured or site-calibrated gamma. Oh-1992 is available for
+  bare soil only within its measured 10-70 degree incidence, L/C/X-band, and 0.1 <= ks <= 6
+  domain. It must not be substituted for a low-grazing Billingsley model. Sea clutter uses the NRL
+  model within its stated frequency, grazing-angle, and sea-state domain.
 - Scenario-level ECM does not yet model receiver saturation or waveform-dependent processing.
 - Pulse-Doppler synthesis starts from a range-compressed representation.
 - SAR/ISAR processing is not yet quantitatively validated.

docs/physics.md

@@ -230,15 +230,39 @@ $$ V_{blind} = n \cdot \frac{\lambda \cdot PRF}{2} $$
 
 ## Environmental Clutter
 
-RadarSim models clutter as a statistical process that degrades the effective SNR: $SNR_{eff} = SNR - CNR$.
+Surface clutter is formed from normalized backscatter and the common range-antenna resolution
+cell:
 
-### Clutter Models
+$$\sigma_c = \sigma^0 A_c$$
 
-| Type | Distribution | Typical $\sigma^0$ |
-|------|--------------|-------------------|
-| **Ground** | Weibull | -15 to -25 dB |
-| **Sea** | GIT Model (Douglas Sea State) | -20 to -45 dB |
-| **Rain** | Marshall-Palmer | Variable ($Z=200R^{1.6}$) |
+Thermal noise and clutter are combined as powers rather than subtracting CNR in decibels:
+
+$$\frac{1}{\mathrm{SINR}} = \frac{1}{\mathrm{SNR}} + \frac{\sigma_c}{\sigma_t}$$
+
+The result record exposes the selected model, $\sigma^0$, cell area, clutter RCS, and resulting
+SINR loss.
+
+### Land clutter
+
+The default constant-gamma screening model is
+
+$$\sigma^0 = \gamma\sin\psi$$
+
+where $\psi$ is grazing angle. Terrain-name gamma values are nominal priors and do not encode
+site moisture, season, cultivation, frequency, or spatial resolution. A measured `gamma_db`
+should be supplied for quantitative studies.
+
+The optional Oh-Sarabandi-Ulaby 1992 model computes HH, VV, and HV normalized backscatter from
+complex relative permittivity and RMS surface height. Its enforced domain is bare soil, L/C/X
+band, 10-70 degree incidence, and $0.1 \le ks \le 6$. It is not used below that incidence-domain
+boundary and is not presented as a low-angle land-clutter model.
+
+### Sea and rain clutter
+
+Sea reflectivity uses the NRL five-parameter model for HH/VV polarization, 0.5-35 GHz,
+0.1-60 degree grazing, Douglas sea states 0-6. Rain volume reflectivity uses the
+Marshall-Palmer $Z=200R^{1.6}$ relationship. These amplitude models are separate from
+ITU-R P.838 path attenuation.
 
 **Implementation:** `src/physics/clutter.py`
 
@@ -260,6 +284,13 @@ RadarSim models clutter as a statistical process that degrades the effective SNR
 4. **Blake, L.V.** (1980). *Radar Range-Performance Analysis*. Artech House.
    - 4/3 Earth model origins
 
+5. **Oh, Y., Sarabandi, K., & Ulaby, F.T.** (1992). "An Empirical Model and an
+   Inversion Technique for Radar Scattering from Bare Soil Surfaces." IEEE TGRS.
+   DOI: 10.1109/36.134086.
+
+6. **Gregers-Hansen, V., & Mital, R.** (2012). "An Improved Empirical Model for Radar
+   Sea Clutter Reflectivity." IEEE TAES. DOI: 10.1109/TAES.2012.6324732.
+
 ---
 
 *Document generated for RadarSim v1.0.0*

src/io/exporter.py

@@ -150,6 +150,17 @@ def _extract_environment(engine) -> Dict[str, Any]:
             "terrain_type": str(getattr(engine, "terrain_type", "rural")),
             "sea_state": int(getattr(engine, "sea_state", 0)),
         },
+        "ground_surface": {
+            "model": str(getattr(engine, "ground_model", "gamma")),
+            "gamma_db": getattr(engine, "land_gamma_db", None),
+            "relative_permittivity_real": float(
+                getattr(engine, "ground_relative_permittivity", complex(8.0, -0.8)).real
+            ),
+            "relative_permittivity_loss": float(
+                -getattr(engine, "ground_relative_permittivity", complex(8.0, -0.8)).imag
+            ),
+            "rms_height_m": float(getattr(engine, "ground_rms_height_m", 0.01)),
+        },
     }
 
     # ECM settings

src/io/scenario_loader.py

@@ -91,6 +91,11 @@ class EnvironmentConfig:
     terrain_type: str = "rural"
     sea_state: int = 0
     rain_rate_mm_hr: float = 0.0
+    ground_model: str = "gamma"
+    land_gamma_db: Optional[float] = None
+    ground_relative_permittivity_real: float = 8.0
+    ground_relative_permittivity_loss: float = 0.8
+    ground_rms_height_m: float = 0.01
 
     def __post_init__(self) -> None:
         if self.pressure_hpa <= 0.0:
@@ -99,6 +104,14 @@ def __post_init__(self) -> None:
             raise ValueError("Water-vapor density and rain rate cannot be negative")
         if not 0 <= self.sea_state <= 6:
             raise ValueError("sea_state must be between 0 and 6")
+        if self.ground_model not in {"gamma", "oh1992"}:
+            raise ValueError("ground_model must be 'gamma' or 'oh1992'")
+        if self.ground_relative_permittivity_real <= 1.0:
+            raise ValueError("ground_relative_permittivity_real must be greater than one")
+        if self.ground_relative_permittivity_loss < 0.0:
+            raise ValueError("ground_relative_permittivity_loss cannot be negative")
+        if self.ground_rms_height_m <= 0.0:
+            raise ValueError("ground_rms_height_m must be greater than zero")
 
 
 @dataclass
@@ -287,14 +300,26 @@ def _parse_targets(self) -> List[TargetConfig]:
     def _parse_environment(self) -> EnvironmentConfig:
         """Parse environment configuration."""
         env = self.data.get("environment", {})
+        clutter = env.get("clutter", {})
+        ground = env.get("ground_surface", {})
+        gamma_value = ground.get("gamma_db")
 
         return EnvironmentConfig(
             temperature_c=float(env.get("temperature_c", 15.0)),
             pressure_hpa=float(env.get("pressure_hpa", 1013.25)),
             water_vapor_gpm3=float(env.get("water_vapor_gpm3", 7.5)),
-            terrain_type=str(env.get("terrain_type", "rural")),
-            sea_state=int(env.get("sea_state", 0)),
+            terrain_type=str(env.get("terrain_type", clutter.get("terrain_type", "rural"))),
+            sea_state=int(env.get("sea_state", clutter.get("sea_state", 0))),
             rain_rate_mm_hr=float(env.get("rain_rate_mm_hr", 0.0)),
+            ground_model=str(ground.get("model", "gamma")).lower(),
+            land_gamma_db=float(gamma_value) if gamma_value is not None else None,
+            ground_relative_permittivity_real=float(
+                ground.get("relative_permittivity_real", 8.0)
+            ),
+            ground_relative_permittivity_loss=float(
+                ground.get("relative_permittivity_loss", 0.8)
+            ),
+            ground_rms_height_m=float(ground.get("rms_height_m", 0.01)),
         )
 
     def get_config(self) -> Optional[SimulationConfig]:
@@ -411,6 +436,13 @@ def create_simulation_engine(self):
             terrain_type=self._config.environment.terrain_type,
             sea_state=self._config.environment.sea_state,
             rain_rate_mm_hr=self._config.environment.rain_rate_mm_hr,
+            ground_model=self._config.environment.ground_model,
+            land_gamma_db=self._config.environment.land_gamma_db,
+            ground_relative_permittivity=complex(
+                self._config.environment.ground_relative_permittivity_real,
+                -self._config.environment.ground_relative_permittivity_loss,
+            ),
+            ground_rms_height_m=self._config.environment.ground_rms_height_m,
         )
 
         return engine

src/physics/clutter.py

@@ -15,7 +15,7 @@
 """
 
 from enum import Enum
-from typing import Dict, Tuple
+from typing import Dict, Optional
 
 import numba
 import numpy as np
@@ -46,17 +46,13 @@ class SeaState(Enum):
     HIGH = 6  # Very large waves
 
 
-# Terrain parameters: (A, B) for σ0 = A + B*sin(ψ) [dB]
-# Reference: Nathanson, "Radar Design Principles", Table 7.1
-TERRAIN_PARAMETERS: Dict[str, Tuple[float, float]] = {
-    "urban": (-15, 15),
-    "suburban": (-20, 12),
-    "rural": (-25, 10),
-    "forest": (-20, 12),
-    "desert": (-30, 8),
-    "mountains": (-18, 14),
-    "sea_calm": (-40, 5),
-    "sea_rough": (-25, 12),
+LAND_GAMMA_PRIORS_DB: Dict[str, float] = {
+    "urban": -5.0,
+    "suburban": -12.0,
+    "rural": -20.0,
+    "forest": -15.0,
+    "desert": -30.0,
+    "mountains": -8.0,
 }
 
 
@@ -130,11 +126,12 @@ def ground_clutter_sigma0(
         terrain_type: str = "rural",
         frequency_ghz: float = 10.0,
         polarization: str = "HH",
+        gamma_db: Optional[float] = None,
     ) -> float:
         """
         Ground clutter backscatter coefficient (σ0).
 
-        Uses empirical model: σ0 = A + B*sin(ψ) [dB]
+        Uses the constant-gamma engineering model σ0 = γ sin(ψ).
 
         Args:
             grazing_angle_rad: Grazing angle [rad]
@@ -145,24 +142,71 @@ def ground_clutter_sigma0(
         Returns:
             σ0 in dB (dB relative to 1 m²/m²)
 
-        Reference: Nathanson, "Radar Design Principles", Table 7.1
+        Terrain-name values are nominal priors, not site calibration. Pass gamma_db from
+        measured clutter whenever quantitative accuracy is required.
         """
-        sin_psi = np.sin(grazing_angle_rad)
-
-        # Get terrain parameters
-        A, B = TERRAIN_PARAMETERS.get(terrain_type, (-25, 10))
-
-        sigma0_db = A + B * sin_psi
-
-        # Frequency adjustment (σ0 increases ~3 dB per octave above X-band)
-        if frequency_ghz > 10:
-            sigma0_db += 3 * np.log2(frequency_ghz / 10)
-
-        # VV polarization typically 2-4 dB higher than HH
-        if polarization == "VV":
-            sigma0_db += 2.5
+        if not 0.0 < grazing_angle_rad <= np.pi / 2.0:
+            raise ValueError("grazing_angle_rad must be between 0 and pi/2")
+        if frequency_ghz <= 0.0:
+            raise ValueError("frequency_ghz must be positive")
+        if polarization.upper() not in {"HH", "VV"}:
+            raise ValueError("polarization must be 'HH' or 'VV'")
+        if gamma_db is None:
+            try:
+                gamma_db = LAND_GAMMA_PRIORS_DB[terrain_type.lower()]
+            except KeyError as error:
+                raise ValueError(f"unknown terrain type: {terrain_type}") from error
+        return float(gamma_db + 10.0 * np.log10(np.sin(grazing_angle_rad)))
 
-        return sigma0_db
+    @staticmethod
+    def bare_soil_oh1992_sigma0(
+        grazing_angle_rad: float,
+        frequency_ghz: float,
+        relative_permittivity: complex,
+        rms_height_m: float,
+        polarization: str = "HH",
+    ) -> float:
+        """Oh-Sarabandi-Ulaby (1992) bare-soil normalized backscatter."""
+        incidence = np.pi / 2.0 - grazing_angle_rad
+        incidence_deg = float(np.degrees(incidence))
+        if not 10.0 <= incidence_deg <= 70.0:
+            raise ValueError("Oh-1992 model requires 10-70 degree incidence")
+        if not 1.0 <= frequency_ghz <= 10.0:
+            raise ValueError("Oh-1992 measurement domain is L-, C-, and X-band")
+        if relative_permittivity.real <= 1.0 or relative_permittivity.imag > 0.0:
+            raise ValueError("permittivity must use the passive convention eps'-j eps''")
+        if rms_height_m <= 0.0:
+            raise ValueError("rms_height_m must be positive")
+
+        wavelength = SPEED_OF_LIGHT / (frequency_ghz * 1e9)
+        ks = 2.0 * np.pi / wavelength * rms_height_m
+        if not 0.1 <= ks <= 6.0:
+            raise ValueError("Oh-1992 model requires 0.1 <= k*s <= 6")
+
+        root_eps = np.sqrt(relative_permittivity)
+        gamma_0 = abs((1.0 - root_eps) / (1.0 + root_eps)) ** 2
+        root_term = np.sqrt(relative_permittivity - np.sin(incidence) ** 2)
+        r_h = (np.cos(incidence) - root_term) / (np.cos(incidence) + root_term)
+        r_v = (
+            relative_permittivity * np.cos(incidence) - root_term
+        ) / (relative_permittivity * np.cos(incidence) + root_term)
+        gamma_h = abs(r_h) ** 2
+        gamma_v = abs(r_v) ** 2
+
+        sqrt_p = 1.0 - (2.0 * incidence / np.pi) ** (1.0 / (3.0 * gamma_0)) * np.exp(
+            -ks
+        )
+        if sqrt_p <= 0.0:
+            raise ValueError("Oh-1992 co-polarization ratio is outside its physical domain")
+        p = sqrt_p**2
+        q = 0.23 * np.sqrt(gamma_0) * (1.0 - np.exp(-ks))
+        g = 0.7 * (1.0 - np.exp(-0.65 * ks**1.8))
+        sigma_vv = g * np.cos(incidence) ** 3 * (gamma_v + gamma_h) / sqrt_p
+        sigma = {"VV": sigma_vv, "HH": p * sigma_vv, "HV": q * sigma_vv}
+        try:
+            return float(10.0 * np.log10(sigma[polarization.upper()]))
+        except KeyError as error:
+            raise ValueError("polarization must be HH, VV, or HV") from error
 
     @staticmethod
     def sea_clutter_sigma0(
@@ -409,9 +453,8 @@ def generate_clutter_map(
         Returns:
             2D array of clutter power [linear]
         """
-        # Create range and azimuth arrays
+        # Create range array
         ranges = np.linspace(100, max_range_m, range_bins)
-        azimuths = np.linspace(0, 2 * np.pi, azimuth_bins)
 
         # Calculate grazing angles
         clutter_map = np.zeros((range_bins, azimuth_bins))