Train Neural Network for Single Station

Introduction

This project trains a neural network to predict the Maximum Usable Frequency (MUF) 30 minutes ahead for a single ionosonde station. The MUF is the highest radio frequency that can be reflected back to Earth by the ionosphere at a given time and path; accurate short-term forecasts are operationally valuable for HF communication planning.

The model does not predict the absolute MUF directly. Instead, it learns to correct a persistence baseline — the linear extrapolation MUF + 30·slope, where slope is a least-squares fit over the four most recent MUF samples (0, −20, −40, −60 min). The network therefore predicts the residual

Δresidual_30min = ΔMUF_30min − 30 · slope

and the final forecast is reconstructed as MUF_pred = MUF + 30·slope + Δresidual_30min. Framing the problem this way centers the target distribution tightly around zero, so the network only needs to learn the correction on top of what linear inertia already provides — storm onsets, sporadic-E events, and sunrise transitions are the dominant non-zero signal it must capture.

Training data and inference inputs are read from a SQLite database that exposes two tables: muf(timestamp, MUF, ...) and space_indices(timestamp, Dst, Kp, F10obs, F10avg). Pre-trained models are serialized to nn-<station>.ser (e.g. nn-bvj03.ser, nn-caj2m.ser) and can be used directly with predict_muf_30min without retraining.

---

Design

Feature Engineering

Raw observations are joined on timestamp and enriched with seven engineered features before being fed to the network. The full input vector has 14 features divided into two groups:

Cyclic features — encoded as (sin, cos) pairs so that day 365 and day 1, or 23:00 and 00:00, are treated as adjacent rather than maximally distant:

Feature	Description
doy_sin, doy_cos	Day-of-year cyclic encoding, period = 365 days
hour_sin, hour_cos	Hour-of-day cyclic encoding (decimal hours), period = 24 h
cos_sza	Cosine of the solar zenith angle at the station (low-order approximation)

Scalar features — z-scored at training time using per-feature mean and standard deviation computed on the training split; the cyclic rows are left on their natural [−1, 1] scale:

Feature	Description
Dst	Dst geomagnetic index at timestamp (nT)
Dst_lag_60min	Dst nearest-sample lookup 60 min in the past (nT)
Kp	Kp geomagnetic index at timestamp
F10obs	Observed F10.7 solar flux (sfu)
F10avg	81-day average F10.7 solar flux (sfu)
MUF	Current MUF at timestamp (MHz)
MUF_trend_20min	MUF − MUF_lag_20min (MHz)
MUF_trend_40min	MUF − MUF_lag_40min (MHz)
MUF_trend_60min	MUF − MUF_lag_60min (MHz)

MUF lag values are obtained by linear interpolation over the observed MUF series (gap tolerance: 30 min), so the feature is still defined when an ionosonde measurement falls on a non-standard timestamp. The Dst_lag_60min feature uses a nearest-sample lookup (gap tolerance: 1 h) instead of interpolation, matching the hourly update cadence of the Dst index.

The data pipeline is illustrated below:

@startuml
skinparam backgroundColor #FAFAFA
skinparam ArrowColor #555555
skinparam ActivityBorderColor #555555

|Database|
start
:Read **muf** table\n(timestamp, MUF);
:Read **space_indices** table\n(timestamp, Dst, Kp, F10obs, F10avg);
:Read full **Dst** series\n(for 60 min lag lookup);

|Feature Engineering|
:Inner-join muf ∩ space_indices\non timestamp;
:Build MUF linear interpolator;
:Compute cyclic encodings\ndoy_sin, doy_cos, hour_sin, hour_cos;
:Compute cos(SZA)\nfrom lat, lon, timestamp;
:Interpolate MUF lags\n(−20, −40, −60 min);
:Nearest-sample Dst lag\n(−60 min, max gap = 1 h);
:Compute MUF trend features\n(MUF − MUF_lag_*min);
:Compute least-squares MUF slope\nover 4 evenly-spaced samples;
:Compute target\nΔresidual_30min = ΔMUF_30min − 30·slope;
:Drop rows with any missing feature;

|Output|
:Feature DataFrame\n(14 inputs + target);
stop
@enduml

Neural Network Architecture

The model is a two-branch residual network implemented with Flux.jl:

• Deep path: three GELU hidden layers (128 → 128 → 64 units) with LayerNorm after each hidden activation and Dropout(0.15) after each norm. The final layer is Dense(64, 1) with bias. LayerNorm is required to keep activations on scale at this depth; without it, three GELU layers tend to drift and require a much smaller learning rate to stabilize.
• Linear skip: a single Dense(N, 1) without bias. The bias on the deep path's output layer already absorbs the mean offset of the residual distribution, so the skip carries only the linear feature response without a redundant second bias term.

Both paths receive the same 14-dimensional standardized input vector; their scalar outputs are summed to produce Δresidual_30min.

@startuml
skinparam backgroundColor #FAFAFA
skinparam rectangle {
    BackgroundColor #EEF4FB
    BorderColor #4A90D9
}
skinparam ArrowColor #555555

rectangle "Input\n14 features" as IN

rectangle "Dense(14→128, GELU)\nLayerNorm(128)\nDropout(0.15)\nDense(128→128, GELU)\nLayerNorm(128)\nDropout(0.15)\nDense(128→64, GELU)\nLayerNorm(64)\nDropout(0.15)\nDense(64→1)" as DEEP #E8F5E9

rectangle "Dense(14→1)\n(no bias)" as SKIP #FFF8E1

rectangle "+" as ADD #F3E5F5
rectangle "Δresidual_30min" as OUT

IN --> DEEP
IN --> SKIP
DEEP --> ADD
SKIP --> ADD
ADD --> OUT
@enduml

The final MUF forecast is assembled outside the network:

MUF_pred(t + 30 min) = MUF(t) + 30·slope(t) + model(x_norm)

where slope(t) is the closed-form least-squares slope over [MUF_lag_60, MUF_lag_40, MUF_lag_20, MUF].

Training Strategy

The network is trained with AdamW (β₁ = 0.9, β₂ = 0.999) and MSE loss on the Δresidual_30min target. MSE is preferred over Huber because large residuals — storm-time drops, sporadic-E enhancements, sunrise jumps — are rare but operationally important; Huber would cap their gradient contribution and cause the network to under-learn those regimes.

Key regularization and optimization choices:

Mechanism	Default	Purpose
Input noise	σ = 0.05 on z-scored scalars	Data augmentation; cyclic features are left unperturbed
EMA weights	decay = 0.999	Validation is evaluated on the EMA shadow, not noisy live weights
LR-on-plateau	factor 0.5, patience 25 epochs	Automatic step-size reduction when smoothed val-loss stalls
Smoothed val-loss	trailing mean over 5 epochs	Prevents stochastic dips from creating phantom "best" epochs
AdamW decay	1×10⁻⁴	Decoupled weight decay
Dropout	0.15 (training only)	Layer-wise regularization within the deep path

Train / Validation Split

Data is partitioned into alternating chronological blocks (3 weeks training, 1 week validation) rather than a single holdout at the end. This ensures the validation set samples the same solar-cycle phase and seasonal patterns as the training set, giving a more honest estimate of generalization. A 2-hour gap is dropped from both ends of every block so that no training row's 30-minute target falls inside a validation block, and no validation row's 60-minute input window reaches into a training block.

Warm Restart Escalation

When the smoothed validation loss stalls, the training loop applies a two-tier restart strategy. A hall of fame keeps the top-5 models seen during training, sorted by smoothed validation loss, to seed restarts from a known-good region of weight space.

@startuml
skinparam backgroundColor #FAFAFA
skinparam ActivityBorderColor #555555
skinparam ArrowColor #555555

start

:Train epoch;
:Update EMA shadow;
:Evaluate smoothed val-loss;

if (improved?) then (yes)
    :Update best model;
    :Reset patience counters;
    :Update Hall of Fame\n(if min-gap elapsed);
else (no)
    :Increment since_improved;
endif

if (since_improved ≥ patience\nAND hall of fame non-empty\nAND not yet mutated) then (yes)
    :Tier 1 — clone random\nHoF member + Gaussian\nmutation (σ = 0.1·std(w));
    :Reset LR, EMA, val history;
    note right: since_improved keeps climbing
else if (since_improved ≥ 2·patience\nOR (≥ patience AND HoF empty)) then (yes)
    :Tier 2 — full weight\nre-initialization;
    :Reset all counters;
else (no)
    :Continue;
endif

:Next epoch;

stop
@enduml

Tier 1 (mutation) fires once at patience epochs without improvement and clones a randomly chosen hall-of-fame member, perturbing its weights with Gaussian noise scaled by each layer's own weight standard deviation. Tier 2 (full reset) fires if the mutated network still has not recovered after another full patience window, or immediately if the hall of fame is empty. The best model across all restart cycles is always retained.

---

Results for Boa Vista (BVJ03)

The neural network for Boa Vista is trained on the bvj03.db dataset, which contains approximately 12 years of ionosonde observations from the station located at latitude 2.87°N, longitude −60.71°W. The metrics for the model are shown as follows. All errors are in MHz and defined as predicted − observed.

Metric	Model	Baseline
Bias (mean error)	0.0148	0.1010
MAE	1.8464	2.9536
RMSE	2.9604	4.8680
Std of error	2.9604	4.8670
Median error	−0.0573	0.0016
|err| 90th percentile	4.2437	6.8618
|err| 95th percentile	6.2392	10.0604
|err| 99th percentile	11.6146	19.7868
Max |err|	48.6860	79.9616
R²	0.8989	0.7268
Pearson correlation	0.9641	0.8866
Skill score (1 − MSE_model / MSE_baseline)	+0.6302	—

The skill score of +0.63 means the network eliminates 63% of the residual variance left by the linear persistence baseline. The near-zero bias (0.015 MHz) confirms the residual framing works as intended — the model adds a well-calibrated correction on top of the extrapolation rather than a biased shift. The roughly 2× reduction in MAE (1.85 vs 2.95 MHz) and RMSE (2.96 vs 4.87 MHz) holds all the way into the tail percentiles, suggesting the network has learned to anticipate the large-excursion events (storms, sporadic-E, sunrise transitions) that dominate the 95th–99th percentile range.

---

Results for Cachoeira Paulista (CAJ2M)

The neural network for Cachoeira Paulista is trained on the caj2m.db dataset. CAJ2M is a mid-latitude station in southeastern Brazil (approximately 22.7°S, 45.0°W), in contrast to the equatorial BVJ03, so the ionosphere here is dominated by different dynamics.

Metric	Model	Baseline
Bias (mean error)	0.0462	−0.2371
MAE	1.6158	2.7934
RMSE	2.4660	4.6410
Std of error	2.4656	4.6349
Median error	0.0094	−0.1006
|err| 90th percentile	3.4996	5.9250
|err| 95th percentile	4.7175	8.3035
|err| 99th percentile	8.9350	19.9379
Max |err|	38.0016	61.4681
R²	0.9298	0.7512
Pearson correlation	0.9723	0.8995
Skill score (1 − MSE_model / MSE_baseline)	+0.7177	—

The skill score of +0.72 is notably higher than Boa Vista's +0.63, indicating the model captures an even larger fraction of the variance left by the persistence baseline at this mid-latitude site. The bias remains near zero (0.046 MHz), and the baseline's systematic negative bias (−0.24 MHz) is essentially eliminated. The ~2.9× reduction in max absolute error (38.0 vs 61.5 MHz) and the consistently tighter tail percentiles suggest the network generalizes well across the more varied ionospheric regimes present at Cachoeira Paulista.