<figure><img src="image_1.png" /></figure>

HAL
open science

# Ionospheric tomography by Over-The-Horizon radar

Corinna Roy

► To cite this version:

Corinna Roy. Ionospheric tomography by Over-The-Horizon radar. Sciences of the Universe [physics]. Institut de Physique du Globe de Paris, 2014. English. NNT: tel-04021027

HAL Id: tel-04021027

https://hal.science/tel-04021027

Submitted on 9 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>
On the front page: Electron density perturbation over Europe obtained by in-version of real data measured by the OTH radar Nostradamus.
<figure><img src="image_4.png" /></figure>
# Contents

<table><tr><td>1</td><td>Résumé étendue en français</td><td>1</td></tr><tr><td>1.1</td><td>Introduction</td><td>1</td></tr><tr><td>1.2</td><td>Ionosphere</td><td>2</td></tr><tr><td>1.2.1</td><td>Les caractéristiques de l'ionosphère</td><td>2</td></tr><tr><td>1.3</td><td>Les modèles ionosphériques</td><td>4</td></tr><tr><td>1.4</td><td>L'ionosphère et le sondage ionosphérique</td><td>4</td></tr><tr><td>1.4.1</td><td>L'inversion des ionogrammes de rétrodiffusion</td><td>5</td></tr><tr><td>1.5</td><td>La tomographie de l'ionosphère</td><td>6</td></tr><tr><td>1.5.1</td><td>La tomographie de l'ionosphère par GNSS</td><td>6</td></tr><tr><td>1.6</td><td>La tomographie de l'ionosphère par radar transhorizon</td><td>8</td></tr><tr><td>1.6.1</td><td>Le tracé de rayons</td><td>8</td></tr><tr><td>1.6.2</td><td>La théorie de la tomographie de l'ionosphère</td><td>9</td></tr><tr><td>1.6.3</td><td>L'inversion par les moindres carrés</td><td>10</td></tr><tr><td>1.6.4</td><td>Les méthodes "Velocity" et "Velocity&amp;Ray"</td><td>10</td></tr><tr><td>1.7</td><td>Résultats obtenus</td><td>11</td></tr><tr><td>1.7.1</td><td>Tests de résolutions</td><td>12</td></tr><tr><td>1.7.2</td><td>L'approche itérative</td><td>14</td></tr><tr><td>1.7.3</td><td>L'inversion de données réelles</td><td>15</td></tr><tr><td>1.8</td><td>Conclusion et Perspectives</td><td>17</td></tr></table>

<table><tr><td>2</td><td>Earth's ionosphere and ionospheric sounding</td><td>19</td></tr><tr><td>2.1</td><td>Introduction</td><td>19</td></tr><tr><td>2.2</td><td>Discovery of the ionosphere: A brief history</td><td>20</td></tr><tr><td>2.3</td><td>Characteristics of the ionosphere</td><td>21</td></tr><tr><td>2.4</td><td>Ionospheric models</td><td>27</td></tr><tr><td>2.5</td><td>Propagation of electromagnetic waves in the ionosphere</td><td>32</td></tr><tr><td>2.5.1</td><td>From waves to rays: Mathematical approximation</td><td>33</td></tr><tr><td>2.6</td><td>Ionospheric sounding</td><td>35</td></tr><tr><td>2.6.1</td><td>Ionosondes</td><td>35</td></tr><tr><td>2.6.2</td><td>Inversion of backscatter ionograms</td><td>38</td></tr></table>
Contents

<table>
  <tbody>
    <tr>
      <td>6.3.1</td>
      <td>Effect of the feedback on the iterative v&r;-method</td>
      <td>112</td>
    </tr>
    <tr>
      <td>6.3.2</td>
      <td>Effect of the regularization on the iterative v&r;-method</td>
      <td>115</td>
    </tr>
    <tr>
      <td>6.3.3</td>
      <td>Iterative Damped least squares inversion with constraints</td>
      <td>119</td>
    </tr>
    <tr>
      <td>6.4</td>
      <td>Iterating the v-method</td>
      <td>122</td>
    </tr>
    <tr>
      <td>6.4.1</td>
      <td>Effect of the feedback on iterative v-method</td>
      <td>122</td>
    </tr>
    <tr>
      <td>6.4.2</td>
      <td>Effect of regularization on the iterative v-method</td>
      <td>125</td>
    </tr>
    <tr>
      <td>6.5</td>
      <td>Summary & Conclusions</td>
      <td>128</td>
    </tr>
    <tr>
      <td><b>7</b></td>
      <td><b>Real data inversion</b></td>
      <td><b>131</b></td>
    </tr>
    <tr>
      <td>7.1</td>
      <td>Data Preprocessing</td>
      <td>131</td>
    </tr>
    <tr>
      <td>7.2</td>
      <td>Inversion results</td>
      <td>138</td>
    </tr>
    <tr>
      <td>7.2.1</td>
      <td>Spring data set</td>
      <td>138</td>
    </tr>
    <tr>
      <td>7.2.2</td>
      <td>Winter data set</td>
      <td>142</td>
    </tr>
    <tr>
      <td>7.3</td>
      <td>Perspectives</td>
      <td>145</td>
    </tr>
    <tr>
      <td>7.3.1</td>
      <td>Comparison with ionosonde data</td>
      <td>145</td>
    </tr>
    <tr>
      <td>7.3.2</td>
      <td>Combination with GPS</td>
      <td>149</td>
    </tr>
    <tr>
      <td>7.4</td>
      <td>Conclusions</td>
      <td>151</td>
    </tr>
    <tr>
      <td><b>8</b></td>
      <td><b>Conclusions and Perspectives</b></td>
      <td><b>155</b></td>
    </tr>
    <tr>
      <td>8.1</td>
      <td>Conclusion</td>
      <td>155</td>
    </tr>
    <tr>
      <td>8.2</td>
      <td>Perspectives</td>
      <td>158</td>
    </tr>
    <tr>
      <td colspan="2"><b>Acknowledgments</b></td>
      <td><b>162</b></td>
    </tr>
    <tr>
      <td colspan="2"><b>List of Figures</b></td>
      <td><b>166</b></td>
    </tr>
    <tr>
      <td colspan="2"><b>Appendices</b></td>
      <td><b>169</b></td>
    </tr>
    <tr>
      <td><b>A</b></td>
      <td><b>Article: New constraints on the 3D shear wave velocity structure of the upper mantle underneath Southern Scandinavia revealed from non-linear tomography</b></td>
      <td><b>171</b></td>
    </tr>
    <tr>
      <td><b>B</b></td>
      <td><b>Article: Complex deep seismic anisotropy below the Scandinavian Mountains</b></td>
      <td><b>191</b></td>
    </tr>
    <tr>
      <td><b>C</b></td>
      <td><b>Article: Effect of ray and speed perturbations on Ionospheric Tomography by Over-the-horizon radar: A new method</b></td>
      <td><b>217</b></td>
    </tr>
    <tr>
      <td colspan="2"><b>Bibliographie</b></td>
      <td><b>237</b></td>
    </tr>
    <tr>
      <td colspan="2"><b>Acronyms</b></td>
      <td><b>247</b></td>
    </tr>
  </tbody>
</table>
# Chapter 1

# Résumé étendue en français

## 1.1 Introduction

Bien que l'ionosphère représente moins de 1% de la masse de l'atmosphère au-dessus de 100 km, son importance résulte de son influence sur les ondes électromagnétiques (EM) à haute fréquence (HF). Tous les signaux des satellites qui traversent l'ionosphère (télécommunications, altimétrie, radar, GNSS, etc.) sont affectés par sa présence. Par conséquent, la connaissance exacte de la structure de l'ionosphère reste un sujet important en Science et pour les applications militaires. La connaissance fine de l'ionosphère joue un rôle clef pour comprendre les processus tels que, par exemple, l'évolution spatiale et temporelle du plasma et pour explorer la propagation d'ondes comme les vagues planétaires, les marées et les perturbations ionosphériques itinérantes (TID)s, et leur influence dans les mesures des satellites scientifiques, des satellites de télécommunication, ou des satellites de service.

Dans le domaine de la défense où les cibles détectées par le radar doivent être localisées précisément, l'effet de l'ionosphère sur les mesures doit être connu. Pour les applications de grande précision (du centimètre au millimètre) et les applications en temps réel, la correction de ces erreurs est un défi important. Réduire les effets de l'ionosphère permettra d'améliorer les précisions pour différentes applications: navigation par satellites, orbitographie, altimétrie satellitaire, télécommunications...

Néanmoins, les effets de l'ionosphère, s'ils ne sont pas corrigés, peuvent réduire la précision de positionnement, d'estimation de vitesse, et diverses informations estimées au cours du temps, qui sont calculées à partir des signaux reçus. Les erreurs de positionnement peuvent varier de quelques mètres à plusieurs dizaines de mètres.

La précision dans le positionnement par Global Positioning System (GPS) et dans la navigation est profondément affectée par l'ionosphère, car les signaux émis par les satellites sont modifiés sur leur chemin dans la densité de son plasma. Le retard causé
par l'ionosphère est devenu la plus grande source d'erreur dans le positionnement par GPS et dans la navigation.

## 1.2 Ionosphere

L'ionosphère a été découverte par Appleton en 1924: un signal continu émis avec un émetteur de la BBC (British Broadcasting Corporation) à Bournemouth près d'Oxford a été reçu à Cambridge. En mesurant la différence en temps entre l'onde directe, se propageant au sol et l'onde réfléchie dans l'atmosphère, il a pu estimer la hauteur de la couche réfléchissante à 100 km. Appleton l'appelle couche E, et pour cette découverte gagne le prix Nobel de physique en 1947. Dans les années suivantes, avec le développement des radars, l'atmosphère a été explorée de plus en plus. Cette exploration a mené à la découverte d'autres couches et finalement à la définition de l'ionosphère. Le nom ionosphère a été proposé la première fois par le physicien écossais Robert Watson-Watt.

### 1.2.1 Les caractéristiques de l'ionosphère

L'ionosphère est définie par la partie supérieure de l'atmosphère entre 60 et 1000 km. Elle est constituée des molécules neutres, d'ions et d'électrons libres formant le plasma ionosphérique. Les molécules principales existantes sont l'azote (N2), l'oxygène (O2), et l'argon, qui forment 99% de l'air. Les ions et électrons sont générés en continu soit par photo-ionisation, soit par l'impact d'atomes et de molécules neutres, ou indirectement par des réactions chimiques. Cet effet s'équilibre avec la recombinaison des ions avec des électrons, menant à la production des photons ou la décomposition d'une molécule dans des atomes excités.

Le niveau d'ionisation varie avec l'altitude, en fonction du type de rayonnement, des différents types de recombinaison des ions et des différents procédés de transport. Aux altitudes basses, la densité de l'air est haute, ainsi les effets de collision et de recombinaison prédominent, et par conséquent, l'ionisation est faible. (Figure 1.1). Avec l'attitude croissante, la densité de l'air diminue, le collision sont plus rares, mais d'un autre côté, le rayonnent solaire plus fort. Dans cet zone une compétition entre ionisation et recombinaison existe. Encore plus haute, l'ionisation par le rayonnement solaire et les rayons X prédominent, et ionisation persiste après le coucher du Soleil à cause de faible taux de recombinaison.

L'ionosphère est décrite en général par plusieurs couches. On distingue plusieurs maximums d'ionisation que l'on désigne par des couches différentes appelées régions D, E et F.
<figure><img src="image_5.png" /><figcaption>Figure 1.1: Densité des ions dans atmosphère et les zones prédominent de l’ionisation et de recombinaison.</figcaption></figure>

## La région D

La région D est située approximativement entre 60 et 90 km d'altitude au-dessus de la surface de la Terre. La couche est faiblement ionisée et disparaît en quelques minutes, lorsque le rayonnement solaire diminue et la recombinaison avec des ions anéantit tous les électrons. A cause de sa faible densité électronique, les ondes électromagnétiques de haute fréquence (3-30 MHz) ne sont pas réfléchies dans la région D, mais subissent une perte d'énergie aboutissant à leur atténuation.

## Les régions E et E sporadique

La région E se situe au-dessus de la région D (90 km and 130 km), et son ionisation est générée par les rayonnements ultraviolets et X sur les atomes d'oxygène. Le taux de recombinaison est plus lent que dans la région D, à cause de la densité d'air réduite (Figure 1.1). La région E diminue pendant la nuit, quand le Soleil disparaît. Dans certains cas, en particulier pendant les mois d'été, des nuages fins avec une ionisation intense peuvent se former sporadiquement, on les nomme E sporadiques. Ils peuvent durer quelques minutes ou plusieurs heures et réfléchissent des ondes avec des fréquences allant jusqu'à 50 MHz ou plus.
## La région F

La couche F (entre 150 km and 800 km) est divisée de jour en deux sous-couches, F1 et F2 qui se rejoignent après le coucher du Soleil. Elle contient la densité la plus élevée en électrons libres. En comparaison des autres couches ionosphériques, la hauteur et la densité électronique sont très variables pour cette couche, en raison des variations diurnes, saisonnières et solaires du rayonnement du soleil. Concernant la propagation des ondes HF, elle joue un rôle primordial, parce qu'elle permet les communications radio à l'échelle mondiale.

## 1.3 Les modèles ionosphériques

Pour décrire l'état de l'ionosphère, des modèles ionosphériques ont été développés qui fournissent des moyennes mensuelles de la densité électronique, ainsi que la température des ions et la composition pour un lieu donné, en fonction de l'heure, la date, l'activité solaire et l'altitude.

Actuellement, environ 170 modèles ionosphériques sont en cours d'utilisation qui diffèrent par leur degré de complexité, le temps de calcul et leur but principal. En général, ils peuvent être classés en modèles empiriques, numériques, paramétrables, et par leur capacité à prendre en compte des données externes.

Les modèles ionosphériques sont constamment améliorés en incluant en outre de nouvelles données ou de meilleures formulations analytiques pour décrire les effets des tempêtes solaires et la variabilité au jour le jour. Les modèles les plus connus sont les modèles empiriques NeQuick et IRI. La structure de l'ionosphère et les modèles sont développés davantage dans le Chapitre 2.

## 1.4 L'ionosphère et le sondage ionosphérique

Lors de la traversée de l'ionosphère, les ondes électromagnétiques transmises interagissent avec les atomes, molécules et les particules chargées qui composent l'ionosphère. Ces interactions modifient les caractéristiques des signaux: retard temporel, déphasage, déviation angulaire, décalage fréquentiel, atténuation. En HF, les ondes électromagnétiques sont réfractées dans l'ionosphère en raison de la présence d'électrons libres et les signaux transmis peuvent être renvoyés vers le sol. Bien que la réfraction des signaux, le temps de retard et le déphasage introduit par l'ionosphère soient des effets indésirables dans la plupart des cas, ces effets peuvent être utilisés pour obtenir des informations sur la densité d'électrons dans l'ionosphère.

Parmi les principaux instruments utilisés pour étudier l'ionosphère et pour obtenir
des informations sur la densité électronique on y trouve les sondeurs verticaux et les radars transhorizon. Tous les deux utilisent les propriétés de réfraction de l'ionosphère sur les ondes électromagnétiques à haute fréquence. Les sondeurs verticaux mesurent un ionogramme vertical, c'est à dire la hauteur virtuelle de réflexion d'un signal en fonction de la fréquence émise. Les ionogrammes verticaux peuvent être inversés pour obtenir un profil vertical de densité d'électrons au-dessus du sondeur.

Les radars trans-horizon (OTH) sondent l'ionosphère suivant une propagation oblique et fournissent des ionogrammes de rétrodiffusion (Chapitre 3). Le signal émis, après déviation dans l'ionosphère, peut atteindre le sol au-delà de l'horizon radioélectrique, généralement à plusieurs milliers de kilomètres de l'émetteur. Le signal est rétrodiffusé et reçu par le récepteur (localisé près de l'émetteur) après propagation par le même trajet (dans le cas le plus courant). Aujourd'hui, il existe peu de radars transhorizon dans le monde : Nostradamus en France ; Jindalee en Australie ; ROTHR, ROTHR-TX, et ROTHR-PR aux Etats-Unis ; un radar en Chine dont les informations sont très réservées, et le réseau SuperDARN localisé près des 2 pôles).

## 1.4.1 L'inversion des ionogrammes de rétrodiffusion

Les études antérieures sur l'analyse de la densité électronique dans l'ionosphère par radar OTH sont toutes basées sur l'inversion des ionogrammes de rétrodiffusion qui représente l'amplitude du signal rétrodiffusé en fonction de la distance et de l'angle d'élévation ou de la fréquence. Pour estimer les trois paramètres ionosphériques principaux de chaque couche (la fréquence critique, la hauteur du maximum d'ionisation, et la demi-épaisseur de chaque couche), deux approches sont possibles (Chapitre 2):

* soit la modélisation analytique de la courbe principale de rétrodiffusion par une couche ionosphérique quasi-parabolique,

* soit l'utilisation d'un tracé de rayons pour simuler numériquement cette courbe.

Les deux approches ne fournissent que des profils verticaux de la densité d'électron, valable à une position latitude / longitude (au-dessus du sondeur vertical ou au point milieu du chemin). On suppose généralement qu'ils sont valables dans un rayon autour de la position. Les sondeurs verticaux sont souvent absents des régions d'intérêt ou alors leur répartition est trop lacunaire pour obtenir une bonne représentation de l'ionosphère sur une large zone.
## 1.5 La tomographie de l'ionosphère

Les méthodes décrites ci-dessus, basées sur les sondeurs verticaux ou les radars trans-horizons, permettent d'obtenir un profil vertical de la densité électronique dans l'ionosphère pour une position locale sur Terre. Le profil est supposé être valable également pour d'autres positions autour du point de mesure.

Il est souvent nécessaire de connaître l'état *a priori* de l'ionosphère non seulement localement, mais aussi sur une échelle mondiale ou régionale. La distribution du plasma peut être calculée par imagerie ionosphérique. Pour y parvenir, des méthodes de tomographie ionosphérique ont été développées.

La tomographie est une technique qui permet de reconstruire le volume d'un objet (d'un corps humain dans le cas de l'imagerie médicale ou d'une structure géologique en géophysique) à partir d'une série de mesures effectuées à l'extérieur de l'objet. Il s'agit d'une technique non invasive qui permet la visualisation des structures internes de l'objet sans modification de l'objet (Chapitre 4).

Ces mesures permettent une reconstruction de certaines propriétés de l'objet en fonction du type d'informations fournies par les capteurs (la pression acoustique, l'atténuation d'un faisceau lumineux, la variation de la vitesse ou de la polarisation des ondes sismiques)

### 1.5.1 La tomographie de l'ionosphère par GNSS

GNSS (Global Navigation Satellite System) est un terme général pour décrire des constellations de satellites de géopositionnement qui sont opérationnels et ont une couverture mondiale. En Avril 2013, le Système de positionnement global américain (GPS) et le système russe GLONASS sont les seuls GNSS opérationnels. Ils permettent de déterminer une position (longitude, latitude, altitude) avec une grande précision à l'aide des signaux radio transmis par des satellites et reçus par des stations au sol.

Les satellites transmettent des signaux bi-fréquences qui se propagent différemment dans l'ionosphère car l'indice de réfraction est à la fois dépendant de la densité du plasma et de la fréquence du signal propageant. La mesure de déphasage entre les deux signaux permet de calculer l'intégral de la densité sur le chemin suivi par les ondes entre le satellite et le récepteur. On s'appelle cet intégral le TEC (Contenu Electronique Total). Le TEC est généralement utilisé pour un représentation 2D de l'ionosphère (Figure 1.2) en supposant que toute la densité électronique soit concentrée au maximum de l'ionisation (hypothèse de couche fine).

La mesure du TEC, conjointement avec certaines informations *a priori*, permet en outre d'estimer la densité électronique locale de l'ionosphère en utilisant la théorie
<figure><img src="image_6.png" /><figcaption>Figure 1.2: Carte global du TEC pour le jour 112 en 2012 (22 avril) à 10 TU, générée sur la base de 200 données GPS/GLONASS. Source: University of Bern</figcaption></figure>

des problèmes inverses (Figure 1.3).

Les premières études de faisabilité de la tomographie ionosphérique par GPS ont été faites par Austen et al. (1988). Grâce à la quantité croissante de données disponibles avec l'installation de réseaux mondiaux et régionaux de récepteurs GPS au sol, cette méthode a subi un développement continu et s'est fortement améliorée au cours des dernières années. Aujourd'hui, c'est la méthode la plus utilisée pour étudier les variations spatiales et temporelles de l'ionosphère en 3D avec affichage de l'évolution dans le temps (3D + 1).

Bien que le GPS soit un outil puissant pour étudier l'ionosphère, il a aussi ses limites. En dehors des limitations spatiales et temporelles associées aux conditions expérimentales, les ensembles de données sont souvent incomplets en raison de l'intervalle d'échantillonnage fini dans le temps et l'angle de vue limité de chaque récepteur. A cela s'ajoute le problème majeur de la tomographie ionosphérique utilisant les signaux GPS : l'absence de rayons à faible incidence, qui conduit à une faible résolution verticale. Les signaux GPS utilisent des fréquences de l'ordre du gigahertz ce qui limite leur sensibilité à la couche ionosphérique la plus dense (couche
<figure><img src="image_7.png" /><figcaption>Figure 1.3: La tomographie de l’ionosphère par GPS. La densité électronique de l’inversion (gauche) et du modèle (à droite) à 21 TU présentée comme des isosurfaces. Source: (Garcia and Crespon, 2008)</figcaption></figure>

F). La tomographie de l'ionosphère à partir de mesures de radar trans-horizon en HF permet de contourner ces 2 problèmes et étendre la sensibilité tomographique aux couches inférieures.

## 1.6 La tomographie de l'ionosphère par radar trans-horizon

### 1.6.1 Le tracé de rayons

La première étape pour développer une méthode de tomographie de l'ionosphère nécessite une modélisation du temps de propagation des ondes électromagnétiques (EM) dans le plasma ionosphérique. Ceci est réalisé par une modélisation de la propagation des ondes EM dans l'ionosphère en utilisant la théorie des rayons. Le code TDR (Tracé De Rayon), utilisé ici, modélise la propagation des ondes EM dans l'ionosphère sur la base de la géométrie optique. Cela signifie qu'il modélise les ondes EM comme des rayons se propageant dans un plasma de densité locale Ne, qui affecte l'indice de réfraction. La propagation est faite dans un système tridimensionnel cartésien centré sur une terre ellipsoïdale entourée par une ionosphère tridimensionnelle obtenue par le modèle empirique NeQuick. Le modèle est donc totalement tridimensionnelle et prend un compte les hétérogénités latérales de l'ionosphère.

Nous utilisons l'équation Eikonal, qui décrit la propagation des rayons dans un milieu et nous la résolvons par la méthode de Runge-Kutta (Press et al., 1992) au
quatrième ordre. Le résultat fournit le chemin suivi par les rayons dans l'ionosphère et donc, le temps de propagation.

<figure><img src="image_8.png" /><figcaption>Figure 1.4: Exemple de Tracé de rayon dans l'ionosphère générée pour Octobre à midi en utilisant le modèle NeQuick. Pour chaque fréquence (6 MHz, 8 MHz, 10 MHz, 12 MHz, 14 MHz and 16 MHz), l'angle d'élévation varie entre 10° and 60°.</figcaption></figure>

### 1.6.2 La théorie de la tomographie de l'ionosphère

Comme brièvement évoque auparavant, le temps de propagation des ondes EM est influencé par la densité d'électrons Ne, qui affect l'indice de réfraction $n(\vec{r})$ selon la formule suivante :

$$n(\vec{r}) = \sqrt{1 - \frac{e^2 N_e(\vec{r})}{4\pi^2 \epsilon_0 m_e f_e^2}} = \sqrt{1 - \frac{80.6 N_e(\vec{r})}{f_e^2}} = \sqrt{1 - \frac{f_p^2}{f_e^2}} \quad (1.1)$$

Toute variation du modèle ionosphérique introduit une variation de la vitesse de propagation des ondes et donc une modification du temps de propagation. La différence entre le temps de propagation dans un modèle ionosphérique $T_{phase}^{synth}$ et dans l'ionosphère réelle $T_{phase}^{real}$ est liée à une perturbation dans la densité électronique $\delta N_e(\vec{r})$

$$\delta T_{phase} = T_{phase}^{real} - T_{phase}^{synth} = -\frac{40.3}{c f_e^2} \int_{s_0} ds \delta N_e(\vec{r}) ds, \quad (1.2)$$

avec la vitesse de la lumière $c$, la fréquence du signal $f_e^2$ et le chemin du rayon $ds$. Le problème inverse qu'il convient de résoudre est donc de calculer la perturbation dans
la densité électronique $\delta N_e(\vec{r})$ à partir d'une différence dans le temps de propagation mesuré, $\delta T_{phase}$.

Afin de résoudre ce problème inverse, trois méthodes d'inversion ont été développées et appliquées: deux méthodes linéaires (inversion par les moindres carrés (Section 4.5.2)) avec et sans contraintes (Section 4.5.3)) et une méthode non-linéaire (méthode du gradient conjugué, Section 4.6.5). Les méthodes sont explicitement développés avec l'ensemble des détails analytiques dans le Chapitre 4 et également dans (Roy et al., 2014) contenu dans l'appendice de ce manuscrit.

Par la suite ont détaille que les méthodes linéaires inversés par les moindres carrés.

### 1.6.3 L'inversion par les moindres carrés

Dans l'inversion par les moindres carrés, l'ionosphère a été paramétrée par $N$ blocs homogènes, ne se recouvrant pas en distance et en altitude. Les problèmes inverses sont souvent des problèmes mal conditionnés parce que le nombre de paramètres à estimer est supérieur au nombre d'observations. Par conséquent, la solution n'est pas unique et de petites variations dans les données peuvent introduire de grands changements dans la solution. Donc, la solution doit être régularisée, en prenant en compte l'écart dans les données et l'écart au modèle initial $\mathbf{m}^0$. Pour trouver une solution stable, satisfaisant les conditions ci-dessus, Menke (1989), propose une inversion par les moindres carrés en utilisant un paramètre de régularisation $\lambda$.

$$ \frac{\delta \mathbf{m}}{\mathbf{m}^0} = (\mathbf{A}^T \cdot \mathbf{A} + \lambda \cdot \mathbf{I})^{-1} \cdot \mathbf{A}^T \cdot \frac{\delta \mathbf{T}}{\mathbf{T}^0}, \quad (1.3) $$

où $\frac{\delta \mathbf{m}}{\mathbf{m}^0}$ est un vecteur avec $N$ perturbations dans la densité électronique, le vecteur $\frac{\delta \mathbf{T}}{\mathbf{T}^0}$ contient $M$ perturbations dans le temps de propagation et $\mathbf{A}$ représente une matrice géométrique avec $M \times N$ segments du trajet des rayons $ds_{ij}$ du rayon $j$ dans le bloc $i$. $\mathbf{I}$ est la matrice d'identité.

### 1.6.4 Les méthodes “Velocity” et “Velocity&Ray”

Deux méthodes linéaires ont été développées : la méthode $v$ et la méthode $v&r$. La première prend un compte seulement l'effet de la densité électronique sur la vitesse des ondes électromagnétique, c'est à dire, on suppose que le chemin des rayons dans l'ionosphère réelle reste le même que dans le modèle initial. Cependant, cette hypothèse est fausse car le défi majeur dans le développement de la tomographie de l'ionosphère par radar trans-horizon est liée à l'indétermination sur la position des points de rétrodiffusion au sol. Ceux-ci peuvent changer pour un angle d'élévation constant car ils dépendent de la variation de la densité d'électrons dans l'ionosphère.
La déviation du trajet des rayons, introduit par la variation du point de rétrodiffusion, peut ajouter une variation supplémentaire dans le temps de propagation. Pour considérer cet effet, la méthode *v&r* a été mise au point. Cette méthode tient compte de la déviation du trajet des rayons en imposant une perturbation d'électrons connu, aléatoire $dm^*$ dans les cases traversées par les rayons et en calculant l'effet sur le temps de propagation $\delta T^*$. La méthode se base donc sur l'hypothèse démontrée que le rapport $\frac{\delta T^*}{dm^*}$ est constant et indépendant de la perturbation imposée ; cela permet donc le calcul du kernel $K_{ji} = \frac{\delta T^*}{dm^*}$. L'inversion par les moindres carrés décrit dans la section 1.6.3 (l'équation 1.3) est appliqué à la matrice **A** pour la méthode *v* et la matrice $\mathbf{A} + \mathbf{K} = \mathbf{M}$ pour la méthode *v&r*. Le meilleur paramètre de régularisation $\lambda$ est choisi par le critère de la *L*-curve, après avoir déterminé les bornes de calcul de la *L*-curve à l'aide d'une courbe d'erreur. Les résultats sont également résumés par Roy et al. (2014).

## 1.7 Résultats obtenus

Pour valider les méthodes mises au point (*v* et *v&r*), elles ont été appliquées à un ensemble de tests sur signaux synthétiques (perturbations localisées ou en damier). Un cas idéalisé et irréel ou la variation de la densité électronique n'affecte pas de tout la propagation des rayons (appelée ici *frozen ray*) a aussi été calculé pour quantifier les limites de résolution envisageable par la géométrie et la paramétrisation de notre problème. Les résultats de l'inversion de *frozen ray* représentent les meilleures résultats que les méthodes *v* et *v&r* peuvent atteindre. En général, ces tests montrent les meilleures performances de la méthode *v&r* et justifie la prise en compte de la déviation du trajet des rayons dans l'inversion. Seule la méthode *v&r* est capable de trouver la position correcte des anomalies. Pour une itération, la méthode *v* n'est pas suffisante, car elle ne reconstruit pas bien les modèles cibles (Figure 1.5). La quantification de performance montre également que la méthode *v&r* reproduit à 60% le modèle cible, là où la méthode *v* est limitée à 41%. (Roy et al., 2014) démontrent aussi que le méthode itérative améliore les résultats à 70% et 60% respectivement.
<figure><img src="image_9.png" /><figcaption>Figure 1.5: Résultats d’inversion après la première itération pour un test de damier et une perturbation localisée de 0.1 % pour différents paramètres de régularisation en utilisant l’approche des *rayons figés*, la *v-* et la *v&r*-méthode. Le résultat de l’inversion pour le meilleur paramètre de régularisation est présenté dans chaque colonne dans la seconde figure à partir du haut.</figcaption></figure>

1.7.1 Tests de résolutions

Les méthodes *v* et *v&r* ont été appliquées à un ensemble de tests sur signaux synthétiques de tailles et d'amplitudes différentes pour analyser la résolution accessible à
la tomographie de l'ionosphère (Chapitre 5). Dans ces essais, l'influence de la grille, les dimensions de la cellule, ainsi que la plus petite amplitude de la perturbation ont été analysées (Chapitre 5).

Les principales conclusions de ces tests de résolution sont que les deux méthodes d'inversion ($v$ et $v&r$) sont indépendantes de la grille. Ceci a été validé en utilisant le test du damier ou des perturbations localisées. Pour ces deux types de perturbations, les deux méthodes peuvent identifier des perturbations de la densité électronique, même pour des petites grilles de 25 km × 20 km en distance et en altitude. En plus, deux petites perturbations localisées de 165 km × 20 km peuvent être identifiées et séparées par les deux méthodes (Figure 1.6).

Pour la première itération, la méthode $v&r$ donne de meilleurs résultats pour des perturbations en damier et pour une perturbation localisée. Elle permet d'identifier la position exacte de l'anomalie maximale et peut reconstruire plus de détails concernant la perturbation en damier à grande distance. La méthode $v$ localise mal la perturbation, ne peut pas reproduire la forme et dans le cas de l'échiquier, le résultat est fortement amorti.

<figure><img src="image_10.png" /><figcaption>Figure 1.6: Résultats d'inversion pour trois modèles cibles et différentes grilles en utilisant la méthode $v$ et la méthode $v&r$.</figcaption></figure>

Cependant, dans les deux résultats de l'inversion ($v$ et $v&r$), des anomalies néga- tives à grande échelle sont visibles pour une perturbation cible localisée. Ces anoma- lies ne sont ni présentes dans le modèle cible ni dans l'inversion *frozen ray*. (Figure 1.5). Ces anomalies peuvent être expliquées par des effets non linéaires dans la dévi-
ation des chemins des rayons. A temps de propagation constant, une surévaluation positive de la perturbation entraîne l'apparition d'une perturbation négative à une autre position.

Tous les essais illustrent la dépendance de la solution à la couverture des rayons et à la zone de sensibilité maximale dans le milieu. Cette zone correspond aux endroits où la fréquence du plasma (liée à la densité électronique locale) est proche de la fréquence du signal. Le test en damier démontre clairement ces dépendances : les perturbations ne sont pas bien résolues ni à grandes distances, ni à basse altitude, ni dans les zones peu couvertes par les rayons (Figure 1.5).

## 1.7.2 L'approche itérative

Les résultats d'inversion après la première itération pour une perturbation localisée ont montré que pour la méthode *v&r* l'amplitude de la solution est affaiblie, et pour la méthode *v*, la perturbation est mal localisée. Afin d'améliorer les résultats de l'inversion, une technique itérative a été développée pour les 2 méthodes. Les résultats préliminaires pour une approche itérative publiés par Roy et al. (2014) ont montré une décroissance des écarts en temps de propagation et des écarts par rapport au modèle cible avec les itérations, qui permet de reproduire le modèle cible à 70% avec la méthode *v&r* et 60% avec la méthode *v*. Cependant, (Roy et al., 2014) ont aussi mis en évidence une perte des performances de la solution après une nombre critique d'itérations. Cela a conduit à la conclusion que l'approche itérative améliore généralement les résultats de l'inversion, car les écarts sont réduits par rapport à la première itération, mais le système devient instable après un nombre critique d'itérations.

L'instabilité de l'inversion, est plus importante dans la méthode *v&r* que dans la méthode *v*. Ceci peut s'expliquer par le bruit numérique induit par la discrétisation des cellules, et par la non-linéarité de la déviation du trajet des rayons (Chapitre 6).

Pour explorer davantage ce problème et réduire l'instabilité, on a introduit un gain de boucle appelé coefficient de feedback (Chapitre 6). Les effets du feedback et de la régularisation ont été testés pour les deux méthodes. Le premier paramètre règle l'amplitude, la seconde amortit la solution ré-injectée.

Ces deux paramètres ralentissent la convergence dans les deux méthodes d'inversion, mais ils ne réduisent pas l'instabilité. Cependant, la méthode *v&r* atteint l'optimum généralement plus rapidement que la méthode *v*. Ceci peut s'expliquer par le fait que la méthode *v&r* tient compte de la déviation du trajet des rayons dès la première itération, alors que la méthode *v* prend en compte cette déviation avec le tracé de rayon dans le modèle mis à jour après la première itération. L'analyse des temps de calcul a montré que pour une itération de la méthode *v&r*, 15 itérations avec la
<figure><img src="image_11.png" /><figcaption>Figure 1.7: Ecarts en temps de propagation et des écarts par rapport au modèle cible avec les itérations pour la méthode *v* (haut) et la méthode *v&r* (bas). Roy et al. (2014)</figcaption></figure>

méthode *v* peuvent être calculées dans le même temps après parallélisation du code.
Cependant, bien souvent la méthode *v* n'a pas encore atteint un optimum après 30
itérations.

1.7.3 L'inversion de données réelles

Les méthodes *v* et *v&r* sont appliquées avec succès à des données réelles obtenues à
l'aide du radar transhorizon NOveau Système TRAnshorizon Décamétrique Appli-
quant les Méthodes Utilisées en Studio (Nostradamus) mis en œuvre par l'ONERA.
L'inversion produit une carte de la densité électronique dans l'azimut choisi qui per-
met de calculer un profil vertical de la densité d'électrons à toute distance dans la
zone atteinte par le radar (Figure 1.8).

Electron density perturbation obtained by real data inversion. The data as
collected at Over-the-horizon (OTH) radar Nostradamus March 14, 2006, at 1855
UT in four Il s'agit d'une grande amélioration par rapport aux techniques classiques
d'inversion des ionogrammes de rétrodiffusion. En effet, ces techniques ne donnent
qu'un profil vertical au point milieu du trajet des rayons. Pour la première fois,
la distribution réelle de la densité électronique dans l'ionosphère a été obtenue par
l'inversion de données d'un radar transhorizon.

Nous avons montré le fort potentiel de cette méthode de tomographie pour re-
construire la distribution du plasma sur le continent européen. Dans un test prélim-
inaire, les profils verticaux de la densité électronique obtenus par inversion étaient
comparés aux profils fournis par le sondeur vertical de Chilton. La possibilité de
combiner la tomographie par radar transhorizon avec la tomographie par GPS a été
démontrée dans un deuxième test préliminaire sur des données synthétiques (Figure
<figure><img src="image_12.png" /><figcaption>Figure 1.8: Perturbation dans la densité électronique sur l'Europe dans 4 azimut obtenue par l'inversion des données réelles (haut). Profile vertical de la densité d'électrons à 500 km en distance du radar obtenue par l'inversion des données réelles mesuré par Nostradamus en mars 2006 à 247° azimut (bas). La perturbation de la densité d'électrons $\delta\text{Ne}/\text{Ne}$ obtenue par l'inversion est montré dans les images à droit. La perturbation de la densité électronique d'un point dans le tracé de gauche correspond au bloc à partir de l'altitude du point dans l'image à droite.</figcaption></figure>

1.9, et Chapitre 7)
## 1.8 Conclusion et Perspectives

Dans cette thèse, une méthode de tomographie pour le radar transhorizon monostatique a été développée pour la première fois. Cela a eu pour but de contourner les limitations de la tomographie par GPS, liées à la haute fréquence utilisée et l'absence des chemins horizontaux des rayons. Deux méthodes ont été développées, la méthode *v* et la méthode *v&r*. La première prend en compte seulement l'effet de la densité électronique sur la vitesse des ondes électromagnétique. La méthode *v&r* considère aussi la déviation du trajet des rayons, introduit par la variation de la densité du plasma. Cela impose l'indétermination du point de rétrodiffusion, ajoutant ainsi une variation supplémentaire dans le temps de propagation.

Les méthodes ont été validées par des données simulées. Les tests montrent que les deux méthodes peuvent trouver les modèles cibles, mais que la méthode *v&r* peut mieux reconstruire leur position et la forme. Les tests de résolutions effectués ont revelés les possibilités de cette méthode de la tomographie par radar transhorizon et ses limites. Plus important, ils ont montrés l'indépendance des résultats à la grille utilisé.

Une approche itérative améliore les résultats d'inversion, parce que des écarts en temps de propagation et des écarts par rapport au modèle décroissent avec les itérations (Chapitre 6).

La ré-augmentation d'écart par rapport au modèle cible après un nombre critique d'itérations peut être explique par le bruit introduit par la discrétisation du modèle ou par les effets non linéaires du déviation du trajet des rayons.

Après validation aux données simulées, la méthode *v* et la méthode *v&r* ont été appliqué avec succès aux données réelles obtenues à l'aide du radar transhorizon Nostradamus mis en œuvre par l'ONERA. L'inversion produit une carte de la densité électronique dans l'azimut choisi qui permet de calculer un profil vertical de la densité d'électrons à toute distance dans la zone atteinte par le radar, couvrant ainsi l'ensemble de l'Europe.

Les profiles doivent être validé par comparaison avec les profiles des sondeurs verticaux. Cela était fait dans un test préliminaire avec le sondeur à Chilton, que se trouve dans un azimut et en distance sondée par Nostradamus. La comparaison montre que l'écart entre le profile du modèle NeQuick et celui de Chilton peut se réduire avec les perturbations d'électron reconstruites par la méthode *v&r*. Ce test préliminaire extrêmement encourageant peut davantage d'être améliorer en exploitant plus en détail la méthode d'inversion utilisée par l'équipe de sondeur vertical de Chilton. L'utilisation des données brutes des sondeurs pourrait nous permettre de les inverser suivant les mêmes hypothèses que nos méthodes, voir les intégrer dans une inversion conjointe sondeur vertical - radar transhorizon.
La méthode de la tomographie de l'ionosphère par radar trans-horizon peut être effectivement combinée avec d'autres techniques, non seulement les sondeurs verticaux, mais surtout les mesures de TEC obtenues par les stations GPS au sol, ou embarquées sur satellites (occultation). La possibilité d'une inversion conjoint avec le GPS a été démontre dans une test préliminaire (Chapitre 7). L'avantage d'une inversion conjointe est que les deux méthodes peuvent se compléter: La tomographie avec le GPS permet uniquement une bonne reconstruction de la densité d'électrons dans la région F2 (~300 km), et avec le radar transhorizon, une excellent sondage des altitude inférieure à 300 km est possible.

<figure><img src="image_13.png" /><figcaption>Figure 1.9: L'inversion conjoint des données simulées du radar trans-horizon et du GPS (haut) pour un test de damier.</figcaption></figure>
# Chapter 2

## Earth's ionosphere and ionospheric sounding

### 2.1 Introduction

The ionosphere is defined as the upper part of the atmosphere extending above the mesosphere from 60 km altitude up to about 1000 km, with its maximum of ionization around 300 km (Figure 2.1). This part is electrically neutral, but solar radiation produces electrically charged atoms and molecules. The ionosphere is a dynamic and highly variable environment that depends on the activity of the Sun and the effect of geomagnetic storms occurring there.

Although the ionosphere makes up less than 1% of the mass of the atmosphere above 100 km, its great importance results from its influence on electromagnetic (EM) waves. All satellite signals that cross the ionosphere (telecommunication, altimetry, radar, GPS, etc.) are affected by its presence. Hence, the exact knowledge of the structure of the ionosphere remains an important scientific and technological topic. The ionosphere has to be known precisely in Science to understand processes such as, for instance, spatial and temporal evolution of the plasma and to investigate all scales of wave activity ranging from planetary waves and tides down to small-scale Traveling Ionospheric Disturbances (TID)s generated at the ground. Apart from that, in the field of military and defense, where targets detected by radar need to be localized exactly as well as in satellite systems, the effect of the ionosphere on the measurements has to be known precisely. Consequently, the fine knowledge of the plasma density of the ionosphere is extremely important to detect perturbations, TIDs as well as to better understand the ionospheric behavior. In this chapter, the characteristics of the ionosphere as well as the methods of ionospheric sounding will be presented in details.
## 2.2 Discovery of the ionosphere: A brief history

The discovery of the ionosphere goes back to the first suggestion made by the physicist Carl F. Gauss in 1839 (Schunk and Nagy, 2009), who speculated that the upper atmosphere might contain ionized regions to explain observed variations of the magnetic field at the surface of the Earth. Important steps in the discovery of the ionosphere were the proof of the existence of radio waves by Hertz in 1887 and the first transatlantic communication established in 1901 by Guglielmo Marconi with a radio transmission at a frequency of 300 kHz.

In 1902, Kennelly and Heaviside independently explained this achievement by suggesting the presence of a permanent electrically conducting layer high in the atmosphere (Gillmor, 1982). In the following years, theories were developed on how the conducting properties of such a layer can be produced, but it was not until 1924 that Appleton was able to prove the existence of such a layer in the atmosphere. This was achieved using a BBC transmitter at Bournemouth to broadcast a continuous signal to Oxford. By measuring the difference between the direct wave traveling along the ground and the signal reflected in the atmosphere, Appleton and his student Barnett were able to determine the height of the reflecting layer at an altitude of 100 km (Gillmor, 1982). On December 12th 1924, Appleton called this layer the E layer. The reason for his choice is not clear, but one possible explanation is that in describing the reflected radio waves, he used the letter E to represent the electric field of the waves.

In 1925, Breit and Tuve in the USA were the first who measured the reflecting layer directly using the first radar-like device. It transmitted a very short radio pulse and the time it took to be reflected determined (Schunk and Nagy, 2009). This method quickly became a standard technique for measuring the height of the reflecting layers and led to extensive investigations of the atmosphere. Finally, this technique resulted in the discovery of two other ionospheric layers, which were named in alphabetical order D and F layer, based on their location below and above the E layer, respectively. The F region is usually subdivided into F1 and F2 layers.

In 1926, the Scottish physicist Robert Watson-Watt first proposed the term ionosphere (Hagfors and Schlegel, 2001). In the following years, there was a rapid progress in understanding the ionospheric behavior. Appleton developed the equation that describes the complex index of refraction and the polarization for a plane wave propagating in a plasma, taking into account the magnetic field, absorbing effects and collisions of the electrons, known as the Appleton-Hartree equation. A few years later, Sydney Chapman presented his theory of the formation of an ionized layer due to the action of solar Ultraviolet (UV) radiation (Chapman, 1931).
## 2.3 Characteristics of the ionosphere

The Earth's atmosphere is composed primarily of the gases molecular nitrogen ($N_2$) (78 %), dioxygen ($O_2$) (21 %) and Argon (1 %) that make up more than 99 % of the air. All charged atoms in the ionosphere are produced either by photoionization, impact ionization of neutral atoms and molecules or indirectly by ionic-chemical reactions. The rate of ionization depends on neutral densities that decrease with height and incoming solar radiation that increases with height.

The most important quantities needed to understand the ionosphere are electron and ion concentrations as a function of altitude. They are solutions of the continuity equation describing the conservation of mass. Since ions and electrons are continuously created by ionization and also disappear in recombination producing neutral molecules, the continuity equation for these particles is (Rees, 1989)

$$ \frac{\partial N_s}{\partial t} + \nabla (n_s \mathbf{u}_s) = P_s + L_s, \qquad (2.1) $$

where $n_s$ is the number density, $P_s$ the production rate per unit volume, and $L_s$ the loss rate per unit volume for ion species $s$. The second term on the left side, $\nabla (n_s \mathbf{u}_s)$, is a transport term and represents the flux divergence that depends on the velocity $\mathbf{u}_s$. In other words, it describes the difference between the flux of ions entering and leaving a given altitude.

The formation of the ionosphere mostly depends on the Sun radiation and is generated by photoionization. The Sun emits enormous quantities of radiation of all wavelengths that impacts atoms (X) or molecules (XY) in the atmosphere leading to free electrons and positively charged ions, e.g.,

$$ X + h\nu \rightarrow X^+ + e \qquad (2.2) $$

$$ XY + h\nu \rightarrow XY^+ + e. \qquad (2.3) $$

By contrast, the amount of ions decreases by recombination. The most important recombination processes are direct and dissociative recombination described by

$$ X^+ + e \rightarrow X + h\nu \qquad (2.4) $$

$$ XY^+ + e \rightarrow X^* + Y^*. \qquad (2.5) $$

The first process leads to the production of a photon, the second one to decomposition of a molecule forming two atoms in an excited state. The recombination rate depends on ion densities, which are decreasing with height. Therefore, recombination is more important in lower ionospheric layers and dissociative recombination of
electrons with molecular ions is the most important loss mechanism for electrons in the E-region.

In addition, collisions become frequent at low altitude and negative ions may form by attachment of an electron:

$$XY + e \rightarrow XY^{-} \qquad (2.6)$$

Negative ions are generated in the D region (60 to 90 km altitude), where ion chemistry is complicated. They may undergo backward reaction under detachment of the electron:

$$XY^{-} + h\nu \rightarrow XY + e \qquad (2.7)$$

At high latitudes, electric fields, particle precipitation, and field-aligned currents play an important role in ionization. Figure 2.1 illustrates ion densities in the atmosphere. $O^+$ predominates at altitudes around 200 km and $H^+$ starts to increase above 300 km (F-region). In the E-region (~150 km), the $NO^+$ and $O_2^+$ ions prevail.

In general, the ionosphere is described by layers. This is an incorrect description and could be misleading, since ionization occurs over the whole ionosphere, but other physical phenomena, including transport and diffusion, play an important role. The level of ionization varies with altitude, depending on the type of radiation, different types of ion recombination and various transport processes. The maximum ionization may be considered as different layers or, more correctly, regions called D, E and F. At lower altitude, a C region can also be present, but its level of ionization is too low for having any effect on radio signals, therefore it is rarely mentioned. A general overview of ionospheric 'layers' during day and night is shown in Figure 2.2.

## D region

The D region reaching from 60 to 90 km is the lowest region within the ionosphere that affects radio communication signals to any degree. It is generated mainly by H-Lyman-$\alpha$ radiation ($\lambda$ = 121 nm). This radiation corresponds to a UV emission line of the hydrogen atom and is emitted when the electron returns from the second to the lowest energy level. The energy of this radiation is sufficiently high to ionize nitric oxide (NO), which is found in small amounts.

Due to the air density and the resulting high recombination rate of ions, the layer is weakly ionized and disappears within a few minutes at night, when the incident radiation for ionization declines. On the other hand, the collision frequency of electrons and other particles during the day is very high with ca. $10^6$ collisions per second, so High Frequency (HF) radio waves are not reflected in the D layer, but suffer loss of energy leading to their attenuation and decrease of intensity.
<figure><img src="image_14.png" /><figcaption>Figure 2.1: Daytime ion densities in the atmosphere based on mass spectrometer measurements. Source: Luhmann (1995).</figcaption></figure>

## E region

The E layer located above the D layer mainly consists of ionized nitric oxide (NO⁺) and ionized oxygen (O⁺) atoms produced by the action of X-ray and extreme ultraviolet (EUV) solar radiation on oxygen atoms. The ionization reaches a maximum at noon. It diminishes during night as electrons recombine with molecular ions such as ionized molecular oxygen (O₂⁺) and NO⁺ and the ionization source is no longer present. The recombination occurs slower than in the D-region and recombination with metallic ions such as Na⁺ is very inefficient. Regarding the propagation of radio waves, this layer can only reflect radio waves with frequencies lower than about 10 MHz, higher frequencies are slightly attenuated.

## Sporadic E layer

In certain cases, especially during the summer months, thin clouds of intense ionization, the so called sporadic E layer, may sporadically form between 90 km and 130 km altitude. They may last for just a few minutes or up to several hours and reflect signals with frequencies up to 50 MHz and higher (Hagfors and Schlegel, 2001). The layer densities reach up to one order of magnitude greater than the background densities, and this layers contain primarily metallic ions (e.g., Fe⁺, Mg⁺), originating from meteoric sources (Schunk and Nagy, 2009).

In general, sporadic E layers are very narrow (0.6 – 2 km), but they can be found at all latitudes and multiple layers can occur simultaneously, separated by
6–10 km. The mechanisms leading to the formation of sporadic E layers are not well understood. Several phenomena are assumed to give rise to their formation, e.g., electrical storms, auroral activity, and upper atmosphere winds, but the involvement of gravity waves is also discussed (Budden, 1985, Chapter 1).

## F region

The F layer, also known as Appleton layer, contains the highest concentration of free electrons and ions in the atmosphere. Therefore, it is greatly affected by the position of the Sun and varies in the course of a day, falling at night as the radiation from the Sun disappears. During the day, the F layer splits into two sub-layers F1 and F2. After sunset these two layers merge again into one single F layer. Unlike other ionospheric layers, the height of the F layer as well as its electron density is highly variable due to large daily, seasonal, and sunspot-cycle variations leading in combination to a generally highly variable behaviour. During day, the electron density can vary between $10^{11} m^{-3}$ to $10^{12} m^{-3}$ (Figure 2.4).

The F layer acts as a "reflector" of signals in the HF component of the radio spectrum enabling ground based worldwide radio communications to be established. It is the most important region associated with HF signal propagation. As the density of gases at this altitude is much lower, recombination of ions and electrons takes place more rarely, with nearly about a quarter of the rate occuring in the E region.

<figure><img src="image_15.png" /><figcaption>Figure 2.2: Ionospheric layers during day and night. Source: Encyclopaedia Britannica</figcaption></figure>

## Influence of the Sun

Since the Sun is the most important source of ionization and the Earth axis of rotation is inclined with respect to its orbital plane, the ionosphere undergoes diurnal,
seasonal and solar cycle variations. There are two factors influencing the ionosphere, namely the solar zenith angle and the solar radiation. Whenever they change, the ionosphere will change, too. Seasonal and diurnal variations of the ionosphere are related to a solar zenith angle change, while its solar cycle variation corresponds to a change in solar EUV and X-ray radiation fluxes. At solar maximum, solar EUV fluxes are greater than those at solar minimum and, as a consequence, these conditions lead to higher electron densities.

The ionosphere undergoes a diurnal variation due to the Earth rotation. At sunrise, the electron density begins to increase rapidly due to photoionization. It rises slowly further during the day and starts decaying afterwards when the photoionization source disappears. Electron density of the ionosphere can be visualized

<figure><img src="image_16.png" /><figcaption>Figure 2.3: Global Total Electron Content (TEC) maps for day 112 in 2012 (April 22) at 2 UT, 8 UT, 12 UT and 20 UT. This map was generated on a daily basis with data from about 200 GPS/GLONASS sites. Source: University of Bern</figcaption></figure>

by various sensing techniques that will be explained in detail in Section 2.6. All of them take advantage on the refractive properties of the ionosphere on HF radio waves. Ionosondes, for instance, are special radars that transmit vertically short EM pulses, which are reflected in the ionosphere. The electron density above the
sounder can be estimated from the time it takes the EM pulse to return. Ionospheric sounding techniques based on monostatic OTH radars use the signal backscattered from the ground to estimate a vertical profile of electron density at the midpoint of the oblique ray path. More recently, techniques based on Global Positioning System (GPS) were developed that enable imaging of the ionosphere.

Contrary to the previous mentioned techniques, they are based on the inversion of the TEC defined as the integral of electron density along the ray-path between a satellite and a receiver at the ground, which is computed from the phase difference between two signals emitted by the GPS satellite and measured at the receiver at the ground (Mannucci et al., 1998). Because of the very dense coverage of satellites and ground stations, with this technique, it is possible to visualize the ionosphere globally, and consequently the global variation linked to the Sun, too.

Figure 2.3 shows the daily variations of the global TEC. TEC is usually visualized as global or regional two-dimensional maps showing the state of the ionosphere (Mannucci et al., 1998). In addition, the measured TEC in combination with some *a priori* information about the background ionospheric model allows estimation of the local electron density of the ionosphere using the inverse problem theory (see section 4.2). The electron density in the ionosphere can also be visualized as vertical profile as shown in Figure 2.4. It illustrates that the ionization below the F region peak reaches its maximum at noon with the solar zenith angle being smallest and decreases subsequently (Figure 2.4, right).

It is straightforward to conclude that a smaller zenith angle results in a higher electron density. Nevertheless, the electron density in the F2 layer during daytime is greater in winter than in summer, although the solar zenith angle is smaller in summer. This phenomenon, clearly evident in Figure 2.4, left, is called *seasonal anomaly* and is explained by seasonal changes in the ion concentration, where loss of ions by recombination is higher than their generation, so total ionization diminishes (Schunk and Nagy, 2009, Chapter 11). The reason for this phenomena is not clearly understood.

Solar cycles are periodic changes in activity and appearance of the Sun with an average cycle duration of about 11 years. The solar cycle modulates the solar flux as well as the solar wind that is defined as the stream of particles released by the Sun. Perturbations of the solar wind can be caused by specific solar events, particularly Coronal Mass Ejections (CME), were plasma consisting primarily of electrons and protons is ejected by the Sun. When this plasma reaches the Earth, the Interplanetary Magnetic Field (IMF) is compressed by the shock of the traveling mass of solar energetic particles and pressure is exerted on the magnetic field. This causes strong disturbance of the Earth's magnetic field resulting in disruption of communications and navigation systems, intense auroras, damage to satellites as
<figure><img src="image_17.png" /><figcaption>Figure 2.4: Vertical electron density profiles calculated using the ionospheric model NeQuick for the year 2013 (left) at noon and every two hours at March 28, 2014 (right).</figcaption></figure>

well as induced currents in power lines and pipelines that lead to power outages and
corrosion in the case of the most extreme storms.

Another typical perturbation of the ionosphere is introduced by Solar Flares. Solar Flares are intense emissions of UV and X rays by the Sun leading to strong ionization in the D region of the ionosphere which is referred to as sudden ionospheric disturbance (SID). They can occur very suddenly developing within a few seconds and cause strong ionization when the increased high-energy solar radiation reaches the upper atmosphere. The strong ionization can absorb radio waves and interrupt radio communication may be interrupted. Flares are frequent around peaks of sunspot cycle.

## 2.4 Ionospheric models

A large number of different models of the ionosphere have been developed and are used for various purposes, including scientific and practical applications such as correction of ionospheric effects on Global Navigation Satellite System (GNSS) and telecommunication channels. These models differ by their degree of complexity, calculation time and their primary purpose. At present, approximately 170 ionospheric models are in use (Schunk, 2013) and discussed in an overview (American Institute of Aeronautics and Astronautics, 1999). Basically, they can be classified into empirical, physics-based numerical, parameterized, and data assimilation models.

**Empirical models** are based on measurements collected over an extended period of time using in situ and remote methods. The data consist of easily measured parameters of each ionospheric layer such as critical frequencies ($f_oE, f_oF_1, f_oF_2$), peak heights ($h_mE, h_mF_1, h_mF_2$) and half-thicknesses ($y_mE, y_mF_1, y_mF_2$). Subsequently, the collected data are averaged and fitted to simple analytical expressions or orthogonal polynomials in order to construct an electron density profile.
**Physics-based numerical models** (e.g., Sami2 is Another Model of the Ionosphere (SAMI2)) (Huba et al., 2000) describe the distribution of electron density together with ion temperatures and drifts. Hence, they allow to study ionospheric behavior over time, including variation of electron and ion densities and temperatures with altitude, latitude, longitude, and solar cycle as well as the behavior under geomagnetic conditions and at different seasons. These models are principally calculated by solving numerically conservation equations (continuity, momentum, energy, etc.) for the ions and electrons taking into account chemical and transport processes in the ionosphere.

**Parameterized models** such as the Parametrized Ionosphere Model (PIM) (Daniell et al., 1995) are empirical models fitting orthogonal functions to the output obtained from a large number of numerical simulations. It returns ionospheric parameters as well as ion composition on a global scale.

**Data assimilation models**, e.g., Global Assimilation of Ionospheric Measurements (GAIM) (Schunk et al., 2004), include different types of data from various measurements resulting in a real time ionospheric model. Data sources are, for example, (1) electron density profiles of the bottom side ionosphere from a network of ionosondes, (2) TEC between a large network of ground stations and GPS satellites, or (3) TEC from occultations between satellites. Details on ionospheric monitoring techniques will be presented in Section 2.6.

In the following, two widely used empirical ionospheric models, namely International Reference Ionosphere (IRI) and NeQuick, are discussed in detail. A visual comparison of the TEC obtained using this two models is presented in Figure 2.5. In both parts of the figure, the equatorial anomaly (the areas with higher TEC about 20° north and south of the equator) becomes evident. Empirical models benefit from their independence of the evolving theoretical understanding of processes in the ionospheric plasma, but a disadvantage is their dependence on the underlying database.

## The International Reference Ionosphere (IRI)

The IRI model (Bilitza et al., 2014) is the standard reference model for the ionosphere. It is a global empirical model based on worldwide collected ground and space data. Data sources are the worldwide network of ionosondes and storm time ionosondes, the Incoherent Scatter Radar (ISR), incoherent scatter data obtained,
for instance, from the European Incoherent Scatter (EISCAT) radar, and satellite data such as EUV data from the Aeronomy Satellite (AEROS), measurements by the International Satellites for Ionospheric Studies (ISIS) and Alouette satellite provide the basis of the IRI model.

Established in the late 1960s by the Commitee on Space Research (COSPAR) and the Union Radio-Scientifique Internationale (URSI), IRI has been steadily improved over the years using updated data and better modeling techniques. Nowadays, a working group of 58 experts is in charge of developing and improving the model. A detailed list of improvements provided since 1968 can be found in (Bilitza et al., 2014).

For a given location, time, date and solar activity, IRI provides monthly averages of electron density and temperature, ion temperature and ion composition in an altitude range from about 60 km to about 2000 km. Thus, IRI can describe monthly varying electron densities, but no day-to-day variability that requires real-time data and an update or assimilation technique combining IRI with these data. This is briefly discussed in (Bilitza et al., 2014).

To obtain the electron density with IRI, the vertical profile is described by 7 subsections: the D region, the E-bottomside, the E-valley, the intermediate region between E and F1 layer, the F1 layer, the F2-bottomside, and F2-topside (Bilitza, 1990). The topside and the bottom side of the electron density profiles are normalized to the F2 peak density and heights. The latest version of the model, IRI-2012, is also able to describe storm effects in the auroral E-region and includes auroral boundaries that allow a better representation of density and temperature features at these boundaries (Bilitza et al., 2014).

Ionospheric models such as IRI rely on solar indices involving daily and seasonal variations as well as the impact of solar activity on ionospheric conditions. In general, these solar indices are the sunspot number *R* (number of dark spots on the solar disc) and the solar radio flux at 10.7 cm wavelength (F10.7). Both can be observed from the ground and long data records exist. IRI uses the 12-months running mean of the sunspot number observed at the Zurich observatory (*R*<sub>12</sub>) and the IG12 index (Bilitza, 1990). *R*<sub>12</sub> is a smoothed value over the values obtained for six months before and six months after the month for which *R*<sub>12</sub> is computed. IG12 is based on F peak plasma frequencies measured by 13 ionosonde stations and on the linear regression with solar activity taken from the International Telecommunication Union-Radiocommunication (ITU-R) model.
## NeQuick model

The “quick calculation model” NeQuick (Radicella and Leitinger, 2001) is an empirical model because it is based on a model introduced by Di Giovanni and Radicella (1990) taking into account the physical properties of the ionospheric layers. It has been developed by the International Centre for Theoretical Physics (ICTP) in Trieste, Italy, in collaboration with the University of Graz in Austria. NeQuick is a three-dimensional and time-dependent ionospheric electron density model, which provides electron densities in the ionosphere as a function of position and time. The input parameters of this model are the position (longitude, latitude, and height), the period (month and UT), and the solar activity (given by monthly-mean sunspot number $R_{12}$ or 10.7 cm solar radio flux).

The electron density distribution is reproduced analytically up to the F2 layer peak using five semi-Epstein layers. The model uses the peaks of the E, F1, and F2 layers as anchor points modeled from the ionosonde parameters ($f_oE, f_oF_1, f_oF_2$). However, the critical frequency for the F2 layer, $f_oF_2$, is modeled by the Comité Consultatif International pour la Radio (CCIR) maps, the critical frequency for the E layer, $f_oE$, is a formulation modified by John Titheridge, and the critical frequency for the F1 layer, $f_o, F_1$, is assumed to be proportional to $f_oE$ in daytime and 0 during night (Leitinger et al., 2005). The topside is represented by another semi-Epstein layer, with a height-dependent, empirical thickness parameter.

As for the IRI model, several efforts have been made to improve the analytical formulation of the NeQuick model. Leitinger et al. (2005) improved the bottom side description of the model, and Coïsson et al. (2006) made major changes in the topside formulation leading to NeQuick 2 (Nava et al., 2008) that is now adopted by ITU-R recommendation. The NeQuick model is also used to correct the ionospheric delay for the Galileo satellite system (Arbesser-Rastburg, 2006).
<figure><img src="image_18.png" /><figcaption>Figure 2.5: Global ionospheric TEC maps for 26th February 2010 at 12 UT modeled by the empirical models NeQuick and IRI. Although they are different in the way they model the ionosphere, they both distinguish the equatorial anomaly. Source: *Najman and Kos (2014)*</figcaption></figure>
## 2.5 Propagation of electromagnetic waves in the ionosphere

The main topic of this thesis is related to radio signals of high frequency radars, therefore details of the physics regarding the propagation of radio waves in a plasma, known as magneto-ionic theory, are discussed here. The cornerstone of this theory is the formula for the refractive index of an ionized medium in a magnetic field, known as Appleton-Hartree formula (Davies, 1965, Chapter 3). The refractive index determines how much radio waves are bent when propagating through the ionosphere. The complex refractive index is given by:

$$n^2 = 1 - \frac{X}{1 - iZ - \frac{Y_T^2}{2(1-X-iZ)} \pm \left( \frac{Y_T^4}{4(1-X-iZ)^2} + Y_L^2 \right)^{\frac{1}{2}}}, \quad (2.8)$$

where the dimensionless quantities $X$, $Y_T$, $Y_L$ and $Z$ are defined as

$$X = \frac{\omega_0^2}{\omega^2}, \quad Y_T = \frac{\omega_B \sin \theta}{\omega}, \quad Y_L = \frac{\omega_B \cos \theta}{\omega}, \quad Z = \frac{\nu}{\omega} \quad (2.9)$$

with the electron plasma frequency $\omega_0$ and the electron gyro-frequency $\omega_B$

$$\omega_0 = \sqrt{\frac{Ne^2}{\epsilon_0 m}}, \quad \omega_B = \frac{Be}{m}. \quad (2.10)$$

In these equations $\epsilon_0$ denotes the vacuum permittivity, $e$ and $m_e$ are the charge and the mass of the electron, respectively, $\omega$ is the frequency of the radio signal, $\nu$ the electron collision frequency, $\theta$ the angle between the direction of propagation and the magnetic field, and $B$ is the ambient magnetic field strength. The plasma frequency is the characteristic frequency of the medium. In equilibrium, the electric fields of electrons and ions cancel each other out in a plasma. If an electron is displaced from a uniform background of ions due to thermal motion of the particles, an electric field is set up pulling the electron back to its original position. Because of their inertia, the electrons will overshoot and oscillate around the equilibrium position with a characteristic frequency, the plasma frequency.

In equation 2.8 $Z$ is the collision term describing energy loss due to collision of particles. $Y$ are the two propagation modes depending on the orientation of the magnetic field, where the subscripts $L$ and $T$ refer to the longitudinal and transverse components of the magnetic field. The $\pm$ sign in the Appleton-Hartree equation gives two separate solutions for the refractive index, since radio signals are split up due to the magnetic field of the Earth as soon as they penetrate into the ionosphere.
The EM wave propagating perpendicular to the magnetic field lines takes the name *ordinary mode* and is usually indicated with a + sign; following the same approach, the *extraordinary mode* is the EM wave propagating parallel to the magnetic field lines, and it is usually indicated with a - sign. At higher frequencies, the *ordinary* and the *extraordinary* waves often follow very similar paths, but at lower frequencies they diverge and travel along completely different paths through the ionosphere. The existence of two split waves is clearly detectable in ionograms. At high frequencies, $f > 8$ MHz, the collisional term $Z$ can be neglected and the refractive index is purely real or imaginary. Neglecting the magnetic field ($Y_L = Y_T = 0$) too and considering that the angular frequency is linked to the frequency $f$ by $\omega = 2\pi f$, equation 2.8 reduces to

$$n^2 = 1 - \frac{f_p^2}{f_e^2}. \qquad (2.11)$$

At ground level, where the electron density (and consequently the plasma frequency) is zero, the refractive index is 1. It decreases with altitude until it is zero, then the plasma frequency $f_p$ equals the frequency of the signal $f_e$. At this condition $n = 0$ the signal is reflected. If the frequency of the signal is too high, the signal is not reflected and escapes into space.

### 2.5.1 From waves to rays: Mathematical approximation

Typically, the propagation of waves in a medium is described by wave equations. A detailed description of the propagation of radio waves in the ionosphere can be found in (Budden, 1985, Chapter 4). However, here ray theory will be used to represent the wave field as rays. The well known homogeneous wave equation for EM waves is given by

$$\left(\frac{\partial^2}{\partial x} + \frac{\partial^2}{\partial y} + \frac{\partial^2}{\partial z}\right) \phi = \frac{1}{v(\vec{x})^2} \frac{\partial^2 \phi}{\partial t^2} \qquad (2.12)$$

where $v$ is the propagation velocity of the EM waves in the medium and $\phi$ the electric field. For a constant velocity, the solution of equation 2.12 could be written as a plane wave. Since the velocity varies slowly with local position, a solution where the amplitude $A$ and the velocity are functions of position is assumed.

$$\phi(\vec{x}, t) = A(\vec{x}) \cdot \exp i\omega \left( \frac{W(\vec{x})}{v_0} - t \right), \qquad (2.13)$$

with $W(\vec{x})$ as the *Eikonal* and $v_0$ a reference velocity. In an isotropic medium, the wavefronts are given by the surfaces, where $W$ is constant and the rays are orthogonal
trajectories to these surfaces. Substituting this modified solution in equation 2.12 gives two sets of equations. For the real part

$$ \nabla^2 A(\vec{x}) - \frac{\omega^2}{v_0} \nabla^2 W(\vec{x}) A(\vec{x}) = - \frac{\omega^2}{v(\vec{x})^2} A(\vec{x}) \quad (2.14) $$

and for the imaginary part

$$ i \left( \frac{2\omega}{v_0} \nabla A \nabla W(\vec{x}) + A(\vec{x}) \frac{\omega}{v_0} \nabla^2 W(\vec{x}) \right) = 0 \quad (2.15) $$

Equation 2.15 is the transport equation and can be used to compute the amplitude of propagating waves. In order to obtain information on propagation of the EM waves, only the real part is considered, for which follows

$$ \nabla^2 W(\vec{x}) - \frac{v_0^2}{v(\vec{x})^2} = \frac{v_0^2 \nabla^2 A(\vec{x})}{\omega^2 A(\vec{x})}. \quad (2.16) $$

The right hand side of this equation is a ratio of the spatial Laplacian of the amplitude to the product of amplitude and $\omega^2$. For high frequencies (like in the case of OTH radar), this term is small. In that case it follows

$$ \nabla^2 W(\vec{x}) = \frac{v_0^2}{v(\vec{x})^2} = n^2. \quad (2.17) $$

Equation 2.17 is called Eikonal equation. It means that the gradient of a wavefront at a position $\vec{x}$ is equal to the refractive index. The direction of maximum change of the wavefront defines the direction of the wave. Basically, this equation is one of several equivalent theorems in geometrical optics like Snell's law and Fermat's principle.

Further assumptions are required to satisfy the Eikonal equation: from equation 2.16 follows that the variation of the gradient of amplitude $A(\vec{x})$ for a given wavelength is smaller than the amplitude $A(\vec{x})$. This implies that

$$ \lambda_0^2 \frac{\nabla^2 A(\vec{x})}{A(\vec{x})} \ll \nabla^2 W, \quad (2.18) $$

with $\lambda_0 = c_0 \cdot \frac{2\pi}{\omega}$ as a reference wavelength. Consequently, an implicit requirement for the Eikonal equation is that the ratio of amplitude variation to amplitude has to be much smaller than the refractive index

$$ \lambda_0^2 \frac{\nabla^2 A(\vec{x})}{A(\vec{x})} \ll n^2. \quad (2.19) $$

Furthermore, the Eikonal equation (equation 2.17) implies that knowing the refrac-
tive index $n$ allows to reconstruct the direction of the ray by ray tracing. Since the *Eikonal equation* is nonlinear and therefore difficult to solve, the ray path is calculated instead of the wavefront. Based on 2.17, the equation for the ray-path is derived (Lay and Wallace, 1995) and given by

$$ \frac{d}{ds} \left( n \frac{d\vec{x}}{ds} \right) = \nabla n. \qquad (2.20) $$

This is a second order differential equation for the ray path $\vec{x}$ and means that the change of the refractive index is related to the change of ray geometry, or more precisely, the change in ray geometry is proportional to the spatial change in the refractive index $n$. The solution of this equation system is calculated numerically using the fourth order Runge-Kutta method (Press et al., 1992).

## 2.6 Ionospheric sounding

The basic instruments used to study the ionosphere and to obtain information about the electron density in the ionosphere are ionosondes, GPS and radars. All of them take advantage of the refractive properties of the ionospheric plasma on EM waves. However, the information they provide about the electron density in the ionosphere is limited to a certain region, because they use different frequencies and incident angles. In this Section, methods and instruments used to investigate the ionosphere and to obtain information about the electron density are described.

### 2.6.1 Ionosondes

Breit and Tuve developed the first prototype of an ionosonde to measure the height of the ionospheric layers (see Section 2.2), and this device has become a basic instrument in studying the ionosphere. Ionosondes are special radars that transmit vertically short EM pulses with a frequency range of 0.1 - 30 MHz. They consist of an antenna with suitable radiation pattern as well as digital control and data analysis circuits.

The EM signals emitted vertically by the ionosonde propagate in the ionosphere as long as the local plasma frequency is smaller than the emission frequency and are reflected in a certain height depending on the signal frequency. From equation (2.11) clearly follows that signals with higher frequencies penetrate deeper into the ionosphere until they are reflected. Their echos are received by the receiver and analyzed by the control system. Modern Digital Ionospheric Goniometric Ionosonde (Digisonde) measure in addition to the travel time of the signal, angles of arrival, polarization and Doppler frequency shift.
The reflection height that can be calculated from the propagation time is called virtual height of reflection, because the group velocity of the EM waves (velocity with which the envelope of a pulse propagates in a medium) is smaller than the speed of light in the vacuum due to ionospheric dispersion. The virtual reflection height is always larger than the real reflection height in the ionosphere. The propagation time of the signal is related to the virtual height of reflection, $h_v$, by

$$h_v = \frac{1}{2} \cdot c \cdot \Delta t, \qquad (2.21)$$

with c being the speed of light in the vacuum. Usually, the reflected pulse consists of two components, the ordinary and the extraordinary wave with different values of $h_v$ resulting from the influence of the Earth magnetic field on HF propagation (Section 2.5).

Figure 2.6 shows an example of a measured ionogram being the virtual height of reflection as a function of frequency at the Digisonde in Juliusruh, Germany. In an ionogram, ionospheric parameters such as the critical frequency and the maximum electron density of each ionospheric layer can be identified, but it contains no information about the vertical profile of electron density. The relation between the virtual height of reflection and the electron density given by (Reinisch and Xueqin, 1983)

$$h_v = \int_{0}^{h_r} \frac{dh}{\sqrt{1 - \frac{N(h)e^2}{\epsilon_0 m w^2}}} \qquad (2.22)$$

allows inversion of the ionogram to obtain a vertical profile of electron density. A number of methods for ionogram inversion have been proposed and developed since invention of the ionosonde.

One of the most common technique is a model-fitting method or POLynomial ANalysis (POLAN), a FORTRAN program developed by Titheridge (1985), where the graph true height versus plasma frequency is approximated by polynomials. This program solves the inversion problem by breaking up the profile into simpler sections, for which physically acceptable solutions can be found, and using extrapolation and interpolation for the rest. It calculates only a profile for the bottomside (250 km) of the ionosphere. The exponential decrease of electron density in the topside ionosphere is normally described by an α-Chapman or Epstein layer.

Automatic Real-Time Ionogram Scaler with True height (ARTIST) is the second commonly used ionogram inversion program that automatically scales digital ionograms by combining image recognition and analytical function fitting techniques. It
provides a vertical profile of electron density from data of each Digisonde that measures polarization and incidence angles (Reinisch and Xueqin, 1983). This technique has been used at the Digisonde in Germany to produce the ionogram in Figure 2.6.

<figure><img src="image_19.png" /><figcaption>Figure 2.6: Example of an ionogram (virtual height of reflection in km as a function of frequency in MHz) measured at the Digisonde in Juliusruh, Germany. The black solid line shows the vertical electron density profile obtained by inversion up to 250 km. The topside profile (dotted line) is extrapolated using a Chapman layer. The E, F1 and F2 layers can be identified, as well as the ordinary (red) and extraordinary (green) modes.</figcaption></figure>

A major drawback of this ionospheric sounding method is the fact that analysis of ionosonde data allows the reconstruction of vertical electron density profiles only for the bottom side ionosphere up to 250 km altitude for a local position on Earth. Typically, an ionosonde station obtains one ionogram recorded every 15 minutes. Figure 2.7 shows a map of the worldwide distribution of ionosondes constructed from the list provided by the National Oceanic and Atmospheric Administration (NOAA). Only 53 stations out of 420 ionosondes provide ionograms in real time to the scientific community, e.g., the one at Chilton, U.K.
<figure><img src="image_20.png" /><figcaption>Figure 2.7: Worldwide distribution of ionosondes created on the basis of the station map from NOAA as at June 2014.</figcaption></figure>

2.6.2 Inversion of backscatter ionograms

Another well known method is the inversion of backscatter ionograms obtained with OTH radar. This radar transmits a signal in an azimuthal direction and sweeps in frequency or in elevation. The signal is refracted by the ionosphere and backscattered from the ground and then analyzed to determine its energy when returning from different ranges. The results are called backscatter ionograms, which represent the amplitude of the backscattered signal as a function of range and elevation angle or frequency.

An example of a backscatter ionogram is shown in Figure 2.8. The strongest signal in the backscattered ionogram represents the leading edge, that is, for each frequency the ray with minimal group delay from the radar to a point of first contact with the ground. In other words, these are the fastest rays for each frequency. The leading edge is a function of two variables operating frequency and minimal group delay and contains information about ionospheric regions which are located thousands of kilometers away from the transmitter and are in some cases inaccessible, for instance the Arctic or the Antarctic. This turns the backscatter sounding into a powerful tool for investigation of the ionosphere. The backscattered ionograms can be inverted to obtain ionospheric parameters, but in general only the leading edge is used, since a complete inversion is extremely difficult and requires information regarding complex phenomena such as sea backscatter and absorption. This approach imposes a real limit in data exploration, because it neglects the information contained in the total backscattered signal.

The idea of exploiting data of OTH radar in order to obtain information on the
Chapter 2. Earth's ionosphere and ionospheric sounding

<figure><img src="image_21.png" /><figcaption>Figure 2.8: An elevation scan backscatter ionogram from OTH radar Nostradamus. Axis labels are not shown for confidential reasons. Source: *Benito et al. (2008)*</figcaption></figure>

electron density in the ionosphere came up in the 1970s (Hatfield, 1970; Rao, 1974) and was further developed later on (e.g., Ruelle and Landeau (1994); Landeau et al. (1997)). However, these methods use 2D ray tracing in a 1D ionosphere described by quasi-parabolic (QP) layers and invert for the three major ionospheric parameters (critical frequency $f_c$, peak height and semi-thickness for each layer) for a local position on Earth. More recently, *Benito et al. (2008)* developed an inversion method of backscatter ionograms optimized by simulated annealing, which was validated for real data of OTH radar Nostradamus.

Since previous methods were only able to provide ionospheric parameters for one local position, *Fridman and Fridman (1994)* extended the inversion of the leading edge from backscatter ionograms to reconstruct the two-dimensional electron density distribution by using additionally the vertical electron density profile measured over the sounding station. Finally, incorporation of data obtained from eight azimuthal beams and the vertical profile measured by quasi-vertical-incident (QVI) sounders (*Fridman, 1998*) allowed a 3 Dimensions (3D) reconstruction of the ionosphere.

Two results are shown in Figure 2.9 for one day in summer and winter at night-time. Their iterative method is able to produce two-dimensional snapshots of the horizontal structure of the lower F2 region, but it seems to be less effective for other ionospheric layers than the F layer. In addition, four data sets are needed as input for the inversion, that is, the leading edge data (group delay versus frequency and azimuth), the vertical profile of electron density from QVI sounders, and the root-mean-square error of the leading edge and vertical profile. Furthermore, the method
do not account for ray-path deflection, needed to be taken in account for monostatic OTH radar. In monostatic OTH radars, the endpoints of the rays (where the signal is backscattered at the ground) are unknown, so the location of the scattering point at the ground can change for a constant elevation angle and depends on the electron density variation $\delta N_e(\vec{r})$ in the ionosphere. The shift of the scattering point introduces an additional in the propagation time. This thesis focalizes on this effect and takes it into account in the inversion of OTH radar. The complete description of the method is presented in Chapter 4.

<figure><img src="image_22.png" /><figcaption>Figure 2.9: Inversion result for the plasma frequency in MHz at 235 km altitude on December 8, 1994, 2244 UT (left) and on August 5, 1994, 2257 UT at 260 km (right). Source: Fridman (1998).</figcaption></figure>

## TEC measurements by satellites

Beginning in 1960, rockets and satellites came into use and enabled to study the ionosphere from above, instead of using instruments on Earth. Shortly after the launch of the first man-made satellite Sputnik 1 by the Soviet Union in 1957, the development of the world's first navigation satellite system TRANSIT for the U.S. Navy started. The first navigation satellite of this system, which consisted of 15 navigation satellites and eight related research satellites, was launched in 1959. After more than 32 years of continuous successful service to the U.S. Navy, the TRANSIT system was replaced by the American Global Positioning System (GPS) and is now used as the Navy Ionospheric Monitoring System (NIMS).

GPS developed quickly for military purposes thereafter and was a great improvement with higher accuracy and stable atomic clocks on board to achieve precise time transfer. However, it wasn't until a civilian Korean Air airplane carrying 269 passengers, was shot down after mistakenly entering Soviet airspace that the Reagan
Administration in the US had opened up GPS for civilian applications. So aircraft, shipping, and transport could fix their positions and avoid straying into restricted foreign territory.

Since December 2012, the GPS system consists of 32 satellites arranged in orbital planes at approximately 20 200 km height. All of these satellites broadcast at the same two frequencies in the L-band, $f_1 = 1575.42$ MHz and $f_2 = 1227.60$ MHz. By receiving the transmitted signals from several of these satellites, a GPS receiver on Earth can accurately determine its position. Although initially intended for military applications, GPS is considered a dual-use technology, having significant military and civilian applications, e.g. cartography, tectonics, navigation, and ionospheric science.

An alternative to the American GPS is the Russian system Global'naya Navigatsionnaya Sputnikovaya Sistema (GLONASS). Both constellations, which are fully operational and have global coverage, are generally described as Global Navigation Satellite System (GNSS). Other GNSS are under development, for instance the European Union Galileo or the Chinese COMPASS positioning system, which is an expansion of the existing Chinese BeiDou-1 satellite system. Apart from global positioning systems, also regional networks like the Indian Regional Navigation Satellite System (IRNSS) and the Japanese Quasi-Zenith Satellite System (QZSS) are under development.

Measurements with GNSS using two or more frequencies are employed to obtain the TEC on the Line Of Sight (LOS) between transmitter and receiver, regardless of the receiver being a stationary ground station or on board of a satellite. The TEC is proportional to the phase difference introduced by the ionosphere between the two signals and is defined as the integral of electron density $N_e$ along the ray path between satellite and receiver. Usually it is visualized as global or regional two-dimensional map to show the state of ionosphere (Mannucci et al., 1998). The TEC can be inverted for the local electron density and then used in ionospheric tomography (Section 4.2).

The vertical resolution in ionospheric tomography with GNSS is usually poor due to the absence of horizontal ray paths. This problem can be addressed by including TEC occultation data. In GPS occultation measurements, the signal emitted by the GPS satellites is received by a receiver on a Low Earth Orbit (LEO) satellite, orbiting around the Earth at approximately 160 km altitude (Figure 2.10). Along its way through the ionosphere, the signal has been refracted due to the presence of free electrons with the amplitude of refraction depending on the electron density. Since both satellites move over time, this technique allows a vertical scanning of successive layers of the atmosphere. The LEO satellites Constellation Observing System for Meteorology Ionosphere and Climate (COSMIC) and CHallenging Minisatellite
<figure><img src="image_23.png" /><figcaption>Figure 2.10: Principle of GPS occultation measurement. The LEO satellite receives a signal from the GPS satellite that has been refracted along its way through the ionosphere. Source: https://www.image.ucar.edu/DAReS/DART/Research/GPS_Liu/</figcaption></figure>

Payload (CHAMP) are used in GPS radio occultation measurements.

The installation of global (e.g., EUREF Permanent Network (EPN), International GNSS Service (IGS)) and regional (e.g., GPS Earth Observation Network (GEONET) in Japan) ground-based GPS networks greatly increased the amount of available data and ray-path coverage, allowing a 3D reconstruction of the ionosphere with time evolution (3D+1). The IGS collects, archives, and distributes GPS observation data sets from an international network of more than 350 continuously operating dual-frequency GPS stations. The generated orbit and tracking raw data are sent to Operational Data Centers, where they are formatted and forwarded to regional or Global Data Centers, which make the data online available. Apart from tracking data and satellite ephemerides, the IGS provides products such as global ionospheric maps.

GEONET is a Japanese GPS station network operated by the Geospatial Information Authority of Japan and consists of approximately 1240 GNSS sites with an average spatial resolution of 25 km (Seemala et al., 2014). Its equivalent in Europe is the EPN with 247 permanently operating GNSS reference stations, which also provide data in real-time.
# Chapter 3

## Over-the-horizon radars

### 3.1 Introduction

All kinds of radar transmit EM signals that are scattered at a “target”, with a part of the signal energy returning to the emitter, where it is detected by a receiver. Comparing the energies returned and transmitted and measuring the propagation time allows determination of the target range, as well as its direction, corresponding to the radar emission direction. Hence, the acronym RADAR stands for RADAR Detection And Ranging. The velocity of an object can be computed based on the doppler shift of a signal that is scattered at an object.

The distance over which a radar works efficiently is limited by obstacles and, especially, by the curvature of the Earth. The problem of detecting a target at very long range, i.e., more than 1000 km away from the transmitter, can be solved by using the reflective properties of the ionospheric plasma located more than 100 km above the ground.

For this purpose powerful Over-the-horizon (OTH) radars using the ionospheric refraction properties on EM waves in the HF band (3-30 MHz) have been developed. There are two types of OTH radars, namely OTH sky wave and OTH surface wave radar. As the names suggests, the EM waves propagate in the ionosphere or along the ground and the sea surface, respectively. In this thesis, only data obtained from the OTH sky wave radar are used, so this kind of radar will be described more detailed.

As described in section 2.5, EM waves are refracted in the ionosphere due to the presence of free electrons. Since ionization increases with altitude up to 300 km, HF EM waves are bent towards the ground to locations beyond the geometric horizon, typically up to thousands of kilometers away from the transmitter. Therefore, this type of radar is called *Over-the-horizon radar*. Figure 3.1 illustrates the difference between a classical radar and an OTH sky wave radar.
The distance reached with an OTH radar depends on the signal frequency and the elevation angle of the radar, but for a single-path the signal may be backscattered at about 2000 km, and for multi-path at about 4000 km. A detailed description of the working principle of OTH radar is given in *Headrick and Anderson (2008)*.

OTH radars can be monostatic or bistatic. For monostatic OTH radar, transmitter and receiver are located at the same place, and signals emitted by the radar are refracted by the ionosphere and then reach the ground, where they are backscattered in all directions, especially in the incident direction. Presumably they follow the same path back to the radar, where they can be detected. In addition to information about the target, the received signal contains information about the medium it propagated through, so it can be used to study phenomena in its interior. For bistatic OTH radars emitter and receiver are located at different places, and the signal transmitted at one site is detected by a receiver several hundred kilometer away.

<figure><img src="image_24.png" /><figcaption>Figure 3.1: Schematic representation of OTH (red) and classical radar (yellow). Source: *Molinié (2011)*</figcaption></figure>

## 3.2 OTH radar worldwide

Since World War II, the request for identifying a target without being limited by the curvature of the Earth lead to the development of OTH systems all over the world. Some of these systems are described in the following sections, and an exhaustive description can be found in *Liu (2007)*.

OTH systems in the United States of America (USA) are described in detail by *Thomason (2003)*. The first experimental monostatic radar, Magnetic-Drum Radar Equipment (MADRE), was built in the late 1950s by the Naval Research Laboratory (NRL) in Chesapeake, Virginia. Data were recorded on magnetic drum devices (hence the name). Another experimental OTH radar built in the United
States is the Wide-Aperture Research Facility (WARF) in central California. This bistatic radar represents a major improvement as compared to MADRE and its large aperture results in improved azimuthal resolutions (Thomason, 2003). The first attempt of establishing an operational OTH radar was the Anglo-American project AN/FPS-95, also known as Cobra Mist. Unfortunately, after two years of testing the project was cancelled, because the radar did not perform as expected and resolutions in range and azimuth were poor. The second operational US OTH radar was Over-the-horizon-Backscatter (OTH-B) that consists of two bistatic systems located at the East and the West Coast of the US, respectively. The east system transmitter was located at Moscow (Maine) and the receiver in Columbia (Maine). The west system had its transmitter in Christmas Valley (Oregon) and its receiver at Tulelake (California). In 2002, the West coast facilities were downgraded to “cold storage” status and finally, the antenna arrays have been pulled down and removed in 2007.

The US Navy created its own system, the Relocatable Over The Horizon Radar (ROTHR), which was originally intended to keep track of ship and aircraft movements over the Atlantic. The prototype was installed on the Aleutian Island and monitored the eastern coast of Russia between 1991 and 1993. Later, it was moved to Virginia (ROTHR-VA) in order to control illegal drug trade in central America and the Caribbean. Subsequently, additional ROTHRs were installed in Texas (ROTHR-TX) and in Puerto Rico (ROTHR-PR).

In Australia, Jindalee was the first bistatic OTH radar project and started in the 1960s. Transmitter and receiver are located separately in central Australia at Alice Springs. In 2003, Jindalee was combined with two other OTH systems. The first one near Laverton in Western Australia and the other one at Longreach, Queensland, in Eastern Australia. This system is now known as Jindalee Operational Radar Network (JORN). Colegrove (2000) reported on the development of JORN in detail.

In Russia, the first OTH radar system was DUGA-2 starting in 1971. It was followed by the operational system STEEL YARD, also known as Russian Wood-speeker due to its loud and repetitive pulses in the middle of the shortwave radio bands. In 2014, the Russians announced a new system named Container, capable of reaching the Netherlands (Karpenko, 2014).

There are also OTH radars in China, Iran, and Japan, but only few details are known regarding these systems. According to a report, China may have as many as three OTH radars used in an early warning system and one OTH-B radar to provide surveillance of the South China Sea.

In France, two OTH radars are installed. One in Valensole that served for ionospheric research and operated by the University Pierre et Marie Curie (Six et al., 1996) and another one Nostradamus located south of Paris, which will be
described in detail in the following section.

## 3.3 SuperDARN

Super Dual Auroral Radar Network (SuperDARN) is a global scale network of HF and Very High Frequency (VHF) radars, constructed by engineers and scientists of a dozen different countries. Radars are located at mid-high latitudes with fields of view covering the polar regions for the study of ionospheric physics. The primary objective of SuperDARN is to provide direct continuous global-scale observations of the motion of irregularities in the ionospheric plasma density at middle (30° – 55°) and high (> 55°) latitudes (Greenwald et al., 1995).

As of 2013, this network consists of 21 continuously operating radars in the northern hemisphere and 11 radars in the southern hemisphere (Figure 3.2). It is still expanding with new radars joining the network almost yearly. The view field of the radar network covers vast regions of ocean near the poles, which play a crucial role in global climate variability. All radars have an azimuthal resolution of approximately 4° for a transmission frequency of 12 MHz (Greenwald et al., 1995).

SuperDARN has been successful in addressing a wide range of scientific questions concerning processes in the magnetosphere, ionosphere, and thermosphere. Radars measure the backscattered energy, the line-of sight Doppler velocity, and the width of the Doppler spectrum (Baker et al., 2010). The Doppler motion of ionospheric irregularities can be used to infer the strength and direction of the ionospheric electric field as well as to calculate convection maps of ionospheric plasma (plasma convection velocity) on a global scale every 1-2 minutes. Convection maps are important diagnostics of energy transfer from the solar wind into the magnetosphere and ionosphere.

Apart from that, scientific objectives of SuperDARN are the investigation of large scale waves, for instance gravity waves, the motion of ionospheric plasma, ionospheric irregularities, and high-latitude plasma structures. In addition, the roughness at the Earth surface including ocean waves and ice cover can be investigated with SuperDARN, since the scatter from the ground is most intense from water-covered areas and almost extinguished over the ice cap of central Greenland.
<figure><img src="image_25.png" /><figcaption>Figure 3.2: Fields of view of currently operational SuperDARN radars in the northern (left) and southern hemisphere (right). Source: http://superdarn.org/tiki-index.php.</figcaption></figure>

## 3.4 Nostradamus

Nostradamus, an OTH radar operated by the Office National d'Etudes et de Recherche Aérospatiale (ONERA), is located 80 km south of Paris at Dreux. In contrast to other OTH radars e.g., JORN or ROTHR, this radar is monostatic that means the antennas for emission and reception are at the same location.

This french OTH radar consists of 288 biconical antenna elements distributed randomly over the arms of a three-branch star spaced by 120°. Each antenna is 7 meter high and 6 meters wide (Figure 3.3). This star-shaped antenna arrangement allows a 360° coverage in azimuth, and the biconical form of each antenna enables the control of the transmission beam in elevation. These unique characteristics of Nostradamus permit the investigation of a very large area of more than 2000 km range (Figure 3.4) all around Europe. Only the central part of the array (96 antennas) is dedicated to transmission and reception, the entire array is used for reception allowing a great capability in receiving beam forming.

Despite the vast sounded area, the radar conserves good resolution in reception with, for example, approximately 2.35° in azimuth and 5.43° in elevation for a frequency of 11 MHz, because the number of antennas used for receiving is three times larger than that used for emission. The receiving array possesses a broad aperture allowing the scattered signals to be resolved into fine azimuth cells. In addition, timing of the received signal allows resolution of the signal into range cells. The resolution in azimuth and elevation is determined by the beam aperture at -
<figure><img src="image_26.png" /><figcaption>Figure 3.3: Example of an antenna of Nostradamus. Source: Molinié (2011)</figcaption></figure>

3 dB below the maximum. This angle characterizes the angle range of the antenna, into which at least half of the maximal power is radiated. Azimuth and range of the returned signal can be determined with a certain resolution depending on the frequency of the signal, since the antenna beam width is frequency depending. In general, resolution in elevation is better for higher frequencies.

The returned signal can therefore be separated into azimuth-range cells, as shown in Figure 3.4). The resulting range-azimuth resolution cell pattern (Figure 3.4, right) is then treated as a search plane for targets, which would manifest as local maxima of received signal power in a cell relative to the surrounding cells. Local maxima are declared as detections. To separate moving targets from the clutter echo a Doppler processing is used. Tracking the location of these detections over time provides target trajectories. Since the amount of data was very large at the time of construction, the antennas are assembled in sub-arrays in order to reduce the data quantity in the computer. The signal from each sub-array is digitalized and merged together giving a receiving beam.

Nostradamus has two different operating modes - the radar and the sounding mode. Due to its unique antenna arrangement, Nostradamus allows to do backscatter sounding by frequency sweep (BSS) and backscatter sounding by elevation scanning (ELS). In the BSS method, a signal is transmitted at a fixed azimuth and elevation angle and scanned in frequencies between 6 and 28 MHz. A measurement for one frequency and one elevation angle takes around 3.5 s, and a complete measurement (all frequencies and elevation angles) requires approximately 21 min. After digitalization of the received signals, cross-correlation with the emitted signal and a
<figure><img src="image_27.png" /><figcaption>Figure 3.4: Simulations for the capacities of OTH radar Nostradamus: Reached distances ranges (left) and range-azimuth resolution cell pattern example for the red square in the left figure (right). Source: Molinié (2011).</figcaption></figure>

Doppler processing, a graph elevation angle as function of group path is obtained. In the **ELS** operating mode, the radar scans in elevation, but at fixed frequency and azimuth. The radar can perform a panoramic sounding, where the azimuth varies from 0° to 360° for a given frequency and elevation angle, as exemplified in Figure 3.5. In this radar mode, the first objective is the determination of range and distance of targets.

<figure><img src="image_28.png" /><figcaption>Figure 3.5: Variability of the ionosphere: Panoramic sounding for one frequency and one elevation angle. Colors represent the amplitude of the backscatter echo. Axis labels are not shown due to confidential reasons. Source: Bazin et al. (2006)</figcaption></figure>
## 3.5 Scientific applications of OTH radar

Although applications of OTH radar are mainly military and defense-orientated (tracking boats and planes, detection of targets), these radars have scientific relevance, too, as already demonstrated by SuperDARN. As mentioned above, the propagation of HF radio waves in the ionosphere depends on the properties of the ionospheric plasma, therefore, OTH radar can be used to study the ionosphere and phenomena occurring inside. Apart from ionospheric science, OTH radars have applications in several different scientific fields.

One of these fields is oceanography, where this radar is very useful to map surface wind directions and ocean currents. Winds over the ocean generate waves on the water that can be identified by HF OTH radar. The Doppler frequency spectrum produced by echoes from water waves, allows determination of the direction of the waves generated by the wind and hence the direction of the wind itself. In addition, the strength of the waves (indicating the state of the sea, or roughness) can be ascertained (Young et al., 1997).

According to the same principle, it is possible to map ocean currents. This method is based on the fact that HF signals are primarily backscattered from *Bragg* resonant surface waves. These are waves having wavelengths which are exactly half of the wavelength of the radio signal. The echo from the sea surface consists of two sharp *Bragg* lines systematically centered on either side of the zero Doppler if no currents are present. The Doppler position corresponds to the velocity of the waves. Otherwise, in the case of underlying surface currents, the *Bragg* lines are no longer symmetric to zero Doppler, since the current changes the wave velocity (Georges et al., 1998).

Also in astrophysics OTH radar is useful. Thomas and Netherway (1989) showed the potential to use OTH radar for the observation of meteors. In ionospheric science, Occhipinti et al. (2010) and Coïsson et al. (2011) showed that perturbations in the ionospheric plasma related to Acoustic Gravity Waves (AGW) and Internal Gravity Wave (IGW) excited by tsunamis or earthquakes associated Rayleigh waves can be detected using backscatter echoes of monostatic OTH radar. In this thesis the possibility of OTH radar to produce a 3D image of electron density in the ionosphere is discussed.
# Chapter 4

## Ionospheric tomography

### 4.1 Tomography

The term tomography is derived from the ancient Greek words το΄μος (slice) and γράφω (to write). It denotes a technique that allows to reconstruct the volume of an object (a human body in the case of medical imaging or a geologic structure in geophysics) from a series of measurements taken from outside the object. It is a non-invasive technique enabling the visualization of internal structures without destroying the object. These measurements result in a reconstruction of certain properties within the object depending on the type of information provided by the sensors (sound pressure, attenuation of a light beam, variation in velocity or polarization of seismic waves).

From a mathematical point of view, tomography consists of two stages. At first, it requires the development of a forward model describing accurately the physical phenomena that are measured. Then, in a second step, the model or 3D distribution of the properties within the object is determined, based on the forward model.

A simple example explaining tomography is a Sudoku. The objective is to fill in a 9 × 9 grid with digits so that each column and each row contains all numbers from 1 to 9. The forward model is to calculate the sum of each row and each column. Knowing the sums of each row and column, the inverse model is to complete the grid.

In medical tomography physicians use X-rays (λ = 0.01-10 nm) sent in all directions through the body of a patient in order to have a look at organs and bones without surgery. X-rays passing through the patient are attenuated and its intensity is reduced. Repeating this procedure for different angles gives information on how the inside looks like. Computers then combine these images into a 3D picture of the body.

Tomography has many other applications, which are quite different from those
in medicine. An interesting example comes from archaeology, where tomography was used to determine the cause of death of the Pharao Tutankhamun (Hawass, 2005). In seismology, seismic waves propagating through the Earth provide, for instance, information about its 3D velocity structure (e.g., Wawerzinek et al. (2013), Appendix A).

In analogy to the examples given above, it is possible to monitor the state of the ionosphere using tomography. Radio waves propagating through the ionosphere are refracted due to the presence of free electrons, and the measured time delay with respect to an ionospheric model can be used to obtain information on the electron density distribution. Most of ionospheric tomography methods use signals emitted by GPS satellites and received on the ground. In the following section, this method as well as its limitations are discussed in detail, followed by a description of the new developed method of ionospheric tomography by OTH radar.

## 4.2 Ionospheric tomography and the role of GPS

Although initially intended for military applications, GPS is a dual-use technology, meaning that it has significant military and civilian applications. Since GPS data became available for scientific applications, it can be used for ionospheric tomography. Figure 4.1 illustrates the principle of GPS tomography. The GPS satellite emits two signals with the frequencies $f_1 = 1575.42$ MHz and $f_2 = 1227.60$ MHz that traverse the ionosphere and are recorded at the ground stations.

As described above and in Section 3.1, EM waves are perturbed when traveling through the ionosphere. The ionospheric refraction causes amongst others phase delays of the signals, which are measured at the ground GPS receivers. The differential phase shift encountered by dual-frequency satellite signals for each satellite location and receiver pair is proportional to the relative TEC along the corresponding ray-path. The phase difference $\Delta\phi$ between the two signals is correlated to the integral of electron density $N_e$ along the ray path between satellite and receiver by

$$ \Delta\phi \sim \left(\frac{1}{f_1^2} - \frac{1}{f_2^2}\right) \int N_e ds. \qquad (4.1) $$

The sum of electron density $N_e$ along the ray-path is called Slant Total Electron Content (STEC)

$$ \text{STEC} = \int N_e ds, \qquad (4.2) $$

where $ds$ is the ray-path. The measured TEC along with some a priori information about the background ionospheric model allows estimation of the local electron density of the ionosphere. Despite of its indirectness, this method has the capability
<figure><img src="image_29.png" /><figcaption>Figure 4.1: Principle of GPS ionospheric tomography. The GPS satellite orbiting around the Earth is emitting signals that are measured at several receivers at the ground. The measured STEC can be inverted for the local electron density perturbation in the red area of the ionosphere. Source: http://gnss.be/ionosphere_tutorial.php</figcaption></figure>

to provide a complete multidimensional image. This is an advantage as compared to other measuring techniques explained before.

First simulations of the feasibility of GPS topography were carried out by Austen et al. (1988) using TEC values computed for a realistic satellite-receiver configuration and applying the simultaneous iterative reconstruction technique (SIRT) (Gilbert, 1972) algorithm for image reconstruction. He successfully showed that numerical tomography techniques can be used to produce two-dimensional vertical cross-sections of the electron density in the ionosphere. Several other numerical tests have been carried out by others trying different algorithm for image reconstruction (Raymund et al., 1990; Raymund, 1994), until the first experimental measurements were reported.

Andreeva et al. (1990) published the first experimental results using TEC data collected at three receivers in Russia and reconstructed the first two-dimensional vertical cross-section of the ionosphere by satellite radio tomography. Since then, GPS-based computerized ionospheric tomography has been widely used to investigate the temporal and spatial variations of ionospheric structures (Hansen et al., 1997; Hajj et al., 1994; Hernández-Pajares et al., 1998; Bust et al., 2004; Garcia and Crespon, 2008; Ma et al., 2005; Yizengaw et al., 2005; Wen et al., 2007). In addition, the installation of global (e.g., IGS) and regional (e.g., GEONET in Japan)
ground-based GPS networks enabled ionospheric tomography worldwide. In Europe (Mitchell and Spencer, 2003), Russia (Kunitsyn et al., 1994), Antartica (Heaton et al., 1996), China (Wen et al., 2007) as well as at the polar caps (Pokhotelov et al., 2011).

Although the amount of available data is huge due to a multitude of satellites and receivers, ionospheric tomography is an ill-posed problem. In general, the reconstruction of images from measurements is an inverse problem that is often ill-posed because measurements are both an incomplete sampling of information and often corrupted by noise. As a consequence, the solution is often unstable (because it is very sensitive to small data changes) and not unique (more than one reconstruction fits the observed data).

The literature of the last years reveals a lot of different reconstruction algorithms tested for ionospheric tomography to address these problems of instability and non-uniqueness. Kersley et al. (1993); Heaton et al. (1995); Mitchell et al. (1995); Pryse et al. (1995); Vasicek and Kronschnabl (1995) applied the iterative Multiplicative Algebraic Reconstruction Technique (MART) (Gordon et al., 1970) algorithm that attempts to minimize differences between measured and calculated TEC values by modifying the background ionosphere until the differences are acceptably small.

The Algebraic Reconstruction Technique (ART) incorporates some prior information on each pixel in an ionospheric grid and was used by Kunitsyn et al. (1994, 2001). However, the solution obtained with this algorithm is, considerably limited, because it is very sensitive to the initial ionospheric model, used for those cells in the tomographic grid, which are not hit by any ray. To cope with this problem, Wen et al. (2010) recently proposed the constrained algebraic reconstruction technique (CART) algorithm, where cells not hit by any ray extract information from their neighbors. They validated their algorithm with numerical simulations and finally applied it for the reconstruction of the electron density distribution over China. Kunitake et al. (1995) successfully used a modified version of singular value decomposition to reconstruct the TEC over Japan in a magnetically disturbed period.

## Validation of tomographic images

Tomographic images need to be verified to be sure that the obtained image corresponds to the reality. Pryse and Kersley (1992) reported the first observation regarding the verification of a tomographic image by another instrument. They used TEC measurements collected at two receivers in Scandinavia and verified their tomographic image obtained with a simple SIRT algorithm by comparison with measurements of the EISCAT radar. After this first successful implementation of the EISCAT radar in radio tomography, it served also in other experimental results for
validation for tomographic images of electron density in Scandinavia (Nygren et al.,
1997; Walker et al., 1997).

Other authors compared peak electron densities and the height of the layer of
maximum electron density from tomographic images with ionosonde data. Fehmers
et al. (1998) show in their ionospheric tomography experiment over midlatitude
Europe that the tomography undervalues the maximum electron density $N_{e,max}$ by
20 % and overvalues the height of the layer of maximum electron density between
80 – 200 km.

Dear and Mitchell (2007) observed that the peak electron density in simulation
studies is underestimated by 15 % and up to 40 % for real GPS data by their to-
mographic reconstruction with Multi-Instrument Data Analysis System (MIDAS)
compared to ionosonde data. They explain the underestimation by a poor represen-
tation of the vertical profile in the simplified set of basis functions they used in the
reconstruction algorithm to constrain the result. They show that incorporation of
the F2 peak height from two ionosondes and a range of scale heights as input in the
reconstruction algorithm improves the result giving a better match with the peak
electron density of the ionosondes.

Limitations

Although GPS is a powerful tool for studying the ionosphere, theoretical limita-
tions of ionospheric tomography using a LEO satellite-to-Earth configuration have
been discussed in detail by Yeh and Raymund (1991) and Na and Sutton (1994).
Apart from spatial and temporal limitations associated with the experimental ar-
rangements, data sets are often incomplete because of the finite sampling interval in
time and a limited view angle of each receiver. Furthermore, an invariant ionosphere
during the time of measurement was assumed. The major problem of ionospheric
tomography using GPS is the absence of horizontal ray-paths resulting in low ver-
tical resolution. Finally, the high frequency of the GPS signal limits the sensitivity
of the EM waves to the maximum of electron density in the ionosphere that is the
F-region.

Three approaches are described to overcome the limitations in vertical resolu-
tion due to the absence of horizontal ray-paths. The first approach is to incorporate
extra information from other experiments, for instance ionosonde data, which give
information on the lower vertical profile. Heaton et al. (1995) tested the incorpo-
ration of scaled ionograms into the imaging and found improvement in the vertical
profiles. Kersley et al. (1993) showed promising results with the incorporation of
ionosonde data into the reconstruction algorithm. Markkanen et al. (1995) applied
a Bayesian approach to simulated results and incorporated peak heights as a priori
information in the reconstruction algorithm.

The second approach is to make use of the general knowledge regarding the shape
of ionospheric profiles to constrain the reconstruction results and fill the information
gap. The most widely used solution is that proposed by *Fremouw et al. (1992)*, who
applied a set of vertical orthonormal functions, created from ionospheric models
to image the vertical profile. However, existing ionospheric models are partially
incorrect, because they fail to predict the strong day-to-day variations. Therefore,
reconstructions constrained with vertical profiles from ionospheric models are often
inaccurate.

*Fehmers et al. (1998)* proposed a model-independent algorithm compensating the lack of horizontal ray paths with information that does not depend on a specific model. They impose that the electron density cannot be negative at high and low altitudes, and assumed it to be smooth and vertically stratified. These constraints are incorporated in the optimization problem, and the mathematical form of the problems allows additionally incorporation of ionosonde data as constraints in the optimization problem. Their tests show moderate success with an error in the layer-height estimation in the order of 90 km.

The third approach is to combine measurements from different sources in the inversion, i.e., from ionosondes, radar or occultation measurements. For instance, *Rius et al. (1997)* showed that combination of GPS TEC and occultation data improved the vertical resolution, whereas the use of ground data alone is insufficient for vertical reconstruction of electron density.

*Fridman and Nickisch (2001)* developed a method for simultaneous inversion of vertical and backscatter ionograms from OTH radar with TEC measurements collected by receivers operated in the Caribbean to obtain a 3D, smooth ionospheric model over a fixed geographical region. *Fridman et al. (2006)* and *Fridman et al. (2009)* developed a 3D real-time reconstruction method of the ionosphere (GPSII), that is able to use absolute and relative TEC data from ground and space GPS receivers including occultation, vertical TEC from altimeters, and in situ measurements on LEO satellites such as CHAMP as well as electron density profiles from vertical sounders. Finally, *Fridman et al. (2012)* upgraded their GPSII method enabling the incorporation of oblique backscatter sounding data. The combined algorithm produces a dynamic model of electron density for a fixed geographical region. This model is consistent with backscatter leading edge data, vertical sounding data as well as with absolute and relative TEC data from a number of GPS/LEO receivers.

With the installation of global and regional GPS networks the amount of data
increased and the angular as well as the ray path coverage improved, making the 3D
reconstruction of the ionosphere with a time evolution (3D+1) possible. Over Scan-
dinavia, van de Kamp (2013) imaged the ionosphere in 3D for the entire December 2006 using TEC measurements performed at the GPS network GEOTRIM in Finland and together with EISCAT radar data as input to constrain the vertical profile. Their inversion results show that the inversion can better resolve vertical profiles of irregular structures when using the profiles from EISCAT instead of Chapman profiles in the inversion process.

## Mapping of disturbances in the ionosphere

Mitchell et al. (1995) have shown the ability of tomography to image large-scale ionisation depletions, which are known as troughs and generally found on the night-side auroral mid-latitude boundary. Since then, tomographic imaging has been used by several authors to image troughs in the ionosphere. Thus, Heaton et al. (1996) imaged a trough over Antarctica for the first time and Kersley et al. (1997) used ionospheric tomography to map a trough over the U.K and Scandinavia.

However, several papers based on experimental results showing tomographic images of Traveling Ionospheric Disturbances (TID) have been published (e.g., Cook and Close (1995); Markkanen et al. (1995); Pryse et al. (1995). TIDs are large-scale irregularities in electron density, which travel in the form of a wave over large distances without significant change of their shape. Large-scale TIDs have periods of the order of one hour, wavelengths of about 1000 km, and horizontal speeds greater than 250 m s⁻¹ (Schunk and Nagy, 2009, Chapter 11).

They are gravity waves generated usually by temperature change, interaction with the wind with the topography, geomagnetic storms, tsunamis, tropical storms, etc. Although various authors pointed out that tomographic methods cannot resolve TID propagating in the ionosphere, Nygren et al. (1997) showed in their tomographic experiment at four receiving stations in Scandinavia that the tomography is unable to resolve wave fronts lying northward of the receiver chain and explained this with rays being perpendicular to the wave fronts. Saksman et al. (1997) demonstrated mathematically and by simulation how a wave-like TID may be partly invisible to transmitters on the ground due to the experimental setup.

## 4.3 Summary

Different methods of ionospheric tomography using GPS were described in the previous sections and their potential and limitations were discussed. Due to high availability of stations and data this method is very powerful and enables a worldwide imaging of the ionospheric plasma in 3D with time evolution and the detection of perturbations in the ionosphere. Nevertheless, ionospheric tomography by GPS is an ill-posed problem due to spatial and temporal limitations associated with the
experimental arrangements, incomplete data sets and limited view angle of each
receiver that has been addressed in several papers in the past years.

Strong limitations of GPS tomography are the poor vertical resolution due to
the absence of horizontal ray-path and its sensitivity to the region of maximum
ionization in the ionosphere because of the high frequency of the GPS signal. In the
next sections a new ionospheric tomography method based on OTH radar will be
presented to overcome some limitations.

4.4 Ray tracing

The first step in developing a tomography method for OTH radar requires a modeling
of the propagation of EM waves in the ionospheric plasma. There are a number
of studies for modeling the propagation of rays in the plasma, all of them rely on
different hypothesis. Most of these studies assume an isotropic ionosphere neglecting
the Earth's magnetic field and losses in the ionosphere, others are limited to 2D
modeling and based on the assumption that the azimuth does not change along
the ray-path. There are two main types of ray tracing in the ionosphere, called
numerical and analytical ray tracing.

**Analytical** ray tracing mostly applies to a spherically symmetrical ionosphere, whose profiles are described by QP or quasi-cubic (QC) segments. As its name suggests, analytical ray tracing uses explicit equations to define the ionosphere and to determine ray parameters such as ground range, reflection height, phase, and group path. Consequently, analytical ray tracing can only be applied to some ionospheric models, but it is much faster than numerical ray tracing. This makes the method particularly valuable in HF applications, where a large number of rays has to be traced. On the other hand, it is difficult to include the effects of the 3D Earth's magnetic field, horizontal gradients, and the ellipticity of the Earth in analytical ray tracing. Thus, analytic ray tracing is limited to simple and unrealistic ionospheric models. Nevertheless, *Norman and Cannon (1997)* developed an analytical ray tracing code that takes into account horizontal gradients by segmenting the ionosphere horizontally as well as vertically.

**Numerical** ray tracing usually integrates the position and the ray direction for
each ray point from the equations, which describe the propagation of waves in a
medium. The most widely used ray tracing technique in ionospheric science is that
described by *Jones and Stephenson (1975)*. It allows to trace rays in several different
ionospheric models taking into account the Earth's magnetic field, electron collision,
and absorption. *Tsai et al. (2010)* developed a 3D numerical ray tracing technique
on a phenomenological ionospheric model for Taiwan, considering the Earth's mag-
netic field and horizontal gradients in order to simulate synthetic ionograms. More
recently, Azzarone et al. (2012) developed a MATLAB software, IONORT, which allows ionospheric ray tracing in a geocentric spherical coordinate system, taking into account the geomagnetic field and particle-electron collisions. Other numerical ray-tracing programs have been developed by different groups, but they are not necessarily freely available.

Although a lot of methods exists, none of them has been used for this thesis, but the code Tracé De Rayons (TDR) developed by Occhipinti (2006) has been implemented and was improved and developed further during this research work. This code traces rays in a 3D heterogeneous ionospheric model starting at a given point on the ground or into the space in the WGS84 coordinate system (National Imagery and Mapping Agency, 2000). In this work the starting point is always considered as the geographical position of the OTH radar Nostradamus (lat: 48.64°N, lon: 1.08°E). The code neglects the Earth's magnetic field, since its application was tailored for Nostradamus, but the resolution of this radar is not sufficient for resolving the two waves (ordinary and extraordinary). In addition, the radar cannot determine the polarization of the waves, making it impossible to distinguish them. Moreover, the tomography is designed for real time application, so ray tracing should run as fast as possible.

## 4.5 Theory of tomography method

### 4.5.1 v-method

As described in Section 2.5.1, EM waves can be approximated as rays, and their propagation in a medium is described by the Eikonal equation

$$ \frac{d}{ds} \left( n \frac{d\vec{x}}{ds} \right) = \nabla n \qquad (4.3) $$

that relates the change in ray geometry $\vec{x}$ along an arc length $s$ to the variation in the refractive index $n$. For a stationary, isotropic, and horizontal stratified ionosphere, the refractive index depends only on the frequency $f_e$ of the emitted EM signal and on the electron density $N_e(\vec{r})$ (Davies, 1990)

$$ n(\vec{r}) = \sqrt{1 - \frac{e^2 N_e(\vec{r})}{4\pi^2 \epsilon_0 m_e f_e^2}} = \sqrt{1 - \frac{80.6 N_e(\vec{r})}{f_e^2}} = \sqrt{1 - \frac{f_p^2}{f_e^2}}, \qquad (4.4) $$

where $\vec{r}$ is the location vector, $\epsilon_0$ the vacuum permittivity, $e$ and $m_e$ are the charge and the mass of the electron, respectively. The electron density $N_e$ is provided by the semi-empirical ionospheric model NeQuick. This model has been chosen for this
study, because there are only slight differences between the IRI and NeQuick models in the description of the ionospheric bottom side and ray tracing should run as fast as possible. IRI accounts for the magnetic field and therefore runs slower than the quick-run model NeQuick.

The solution of equation 4.3 is calculated numerically using the fourth order Runge-Kutta method (Press et al., 1992), where the step size h in the Runge-Kutta method was adapted to the local elevation angle φ along the ray with

$$h = \max \left( 1000 \cdot \left| \frac{\sin \phi}{\sin 1^{\circ}} \right|, 10 \right), \qquad (4.5)$$

so that it varies between 1000 m in distance for local elevation angles larger than 1° and 10 m at the reflection point. This variable step size ensures a better sensitivity at the reflection points. This is of importance, since rays are most sensitive to the medium in the region around the reflection point, where the plasma frequency $f_p$ and the signal frequency $f_e$ (equation 4.4) are equal. The sensitivity of the rays to the medium is discussed in detail in Section 4.6.2. At their ends, rays are interpolated to altitude zero in order to avoid an artificial time delay introduced by different ray-path lengths.

Figure 4.2 shows traced rays in the ionospheric model NeQuick for six different frequencies and varying elevation angles. The ionosphere was generated for October at 12 UT and a solar flux of 198.1 solar flux units (SFU). Different ray types are clearly visible. Lower angle rays, reaching larger distances up to 2500 km with increasing frequency, and high angle rays that penetrate deeper into the ionosphere with higher frequencies. Additionally, guided rays exists, which get trapped in the ionosphere and move almost parallel to the Earth's surface (e.g., ray for 12 MHz with elevation angle 13° in Figure 4.2). These rays are called *Pedersen Rays* and may be generated by large gradients in electron density.

Knowing the ray-path in the ionospheric model by solving equation 4.3 allows to calculate the propagation time of an EM wave from the radar through the ionosphere and to the ground, where it is backscattered.

$$T_{phase} = \frac{1}{c} \int_{s(n)} n(\vec{r}) ds, \qquad (4.6)$$

where c is the speed of light and $n(\vec{r})$ the refractive index of the medium. By replacing the refractive index in equation (4.6) by equation (4.4) and linearizing by first order Taylor series expansion, the propagation time of an EM wave can be
<figure><img src="image_30.png" /><figcaption>Figure 4.2: Example of rays traced in azimuth 89° in the ionospheric model NeQuick for October at 12 UT and a solar flux of 198.1 solar flux units. For each frequency (6 MHz, 8 MHz, 10 MHz, 12 MHz, 14 MHz and 16 MHz), the elevation angle varies between 10° and 60°.</figcaption></figure>

separated into two parts (Roy et al., 2014)

$$ 
\begin{aligned}
T_{phase} &= \frac{1}{c} \int_{s(n)} \sqrt{1 - \frac{80.6 N_e(\vec{r})}{f_e^2}} ds \\
&\approx \underbrace{\frac{1}{c} \int_{s(n)} ds}_{\text{vacuum}} - \underbrace{\frac{40.3}{c f_e^2} \int_{s(n)} N_e(\vec{r}) ds}_{\text{ionosphere}},
\end{aligned} 
\quad (4.7) 
$$

where the first integral in the second line describes the propagation in vacuo and the second one the delay introduced by the ionosphere.

Figure 4.3 illustrates the limit of linearization comparing the first and the second line of equation (4.7). The difference between exact and linearized propagation time is mostly negligible and is in order of few percent for high elevation angles. Another interesting feature evident in Figure 4.3 is the boundary between E and F layer that separates the down-left corner (E region) from the rest of the figure (F region). The lack of values in the right upper corner (white or navy-blue cells) is associated with rays that have not been reflected in the ionosphere that means they traverse the ionosphere and escape into space.

Equation (4.7) allows to calculate by ray tracing the propagation time $T_{phase}^{synth}$ of any EM wave with a given frequency $f_e$ in any given ionospheric model with electron density $N_e^0(\vec{r})$. The electron density in the real ionosphere can be described
<figure><img src="image_31.png" /><figcaption>Figure 4.3: Validation of the linearization of the refractive index. Exact propagation time calculated with the first line of equation (4.7) (left), the linearized propagation time calculated with the second line of equation (4.7) (middle), and their difference in percent (right) for 1071 traced rays with elevation angles 10° – 60° and frequencies 6–16 MHz. Source: Roy et al. (2014)</figcaption></figure>

by the electron density of an a priori ionospheric model $N_e^0(\vec{r})$ plus a perturbation $\delta N_e(\vec{r})$ with respect to the model, i.e., $N_e^{real} = N_e^0(\vec{r}) + \delta N_e(\vec{r})$. The ray-path is also perturbed in the real ionosphere compared to the ionospheric model. Assuming that $s_0$ is the ray-path in the a priori model $N e^0$, the ray-path in the real ionosphere can be described as $s' = s_0 + \delta s$. The propagation time in the real ionosphere is given by the term

$$T_{phase}^{real} = \frac{1}{c} \int_{s_0+\delta s} ds - \frac{40.3}{cf_e^2} \int_{s_0+\delta s} (N_e^0(\vec{r}) + \delta N_e(\vec{r})) ds, \quad (4.8)$$

which can be separated into

$$T_{phase}^{real} = \frac{1}{c} \int_{s_0} ds - \frac{40.3}{cf_e^2} \int_{s_0} (N_e^0(\vec{r}) + \delta N_e(\vec{r})) ds \\ + \frac{1}{c} \int_{\delta s} ds - \frac{40.3}{cf_e^2} \int_{\delta s} (N_e^0(\vec{r}) + \delta N_e(\vec{r})) ds. \quad (4.9)$$

The integrals in the second line are neglected if starting and endpoints of the ray-paths are fixed (e.g., epicenter and seismometers in the case of seismic tomography).
Since rays are traced between two fixed points and Fermat’s principle states that for fixed endpoints the traveltime along a ray-path is stationary with respect to perturbations in the path. The linearized synthetic propagation in an a priori ionospheric model along the unperturbed ray-path $s_0$ is given by

$$T_{phase}^{synth} = \frac{1}{c} \int_{s_0} ds - \frac{40.3}{cf_e^2} \int_{s_0} N_e^0(\vec{r}) ds. \qquad (4.10)$$

The terms in the second line of equation (4.9) describe the ray-path deflection. Neglecting them, the difference $\delta T_{phase}$ between real and synthetic propagation time (eq. 4.10) in an a priori ionospheric model takes the form Roy et al. (2014)

$$\delta T_{phase} = T_{phase}^{real} - T_{phase}^{synth} = -\frac{40.3}{cf_e^2} \int_{s_0} \delta N_e(\vec{r}) ds. \qquad (4.11)$$

This approach is based on the hypothesis that the ray-path $s(n)$ in the real ionosphere and in the a priori model are the same. That means that the electron density variation $\delta N_e$ only introduces variation in the speed of the EM waves and not in the ray-path. Therefore, this approach is named v-method.

### 4.5.2 Damped least squares inversion

The inverse problem to be solved is to reconstruct the electron density perturbation $\delta N_e(\vec{r})$, which is linked to the time delay between real and synthetic propagation time $\delta T_{phase}$ (equation 4.11). In order to solve numerically equation 4.11, the model is parametrized by $N$ homogeneous, non-overlapping blocks indexed $i$ (Figure 4.4). The electron density perturbation $\delta N_e(\vec{r})$ is expressed by a linear combination of electron density perturbation $\delta m_i$ in block $i$ and $N$ known basis functions $B_i(\vec{r})$ by

$$\delta N_e(\vec{r}) = \sum_{i=1}^{N} \delta m_i \cdot B_i(\vec{r}), \qquad (4.12)$$

where $B_i(\vec{r})$ is defined as

$$B_i(\vec{r}) = \begin{cases} 1 & \text{if } \vec{r} \text{ in block } i \\ 0 & \text{otherwise.} \end{cases}$$

Substituting equation (4.12) into equation (4.11), equation (4.11) becomes for jth-
<figure><img src="image_32.png" /><figcaption>Figure 4.4: Schematic representation of the parametrization of j-rays and i-cells as well as of the ray-path deflection induced by a localized perturbation (e.g., the gray cell). This produces a perturbed ray (red) as compared to the unperturbed ray (black).</figcaption></figure>

measurement of travel time perturbation

$$ \delta T_j = - \frac{40.3}{f_{ej}^2 c} \sum_{i=1}^{N} \delta m_i ds_{ij}, \qquad (4.13) $$

where $ds_{ij}$ is the length of ray-path segments within block $i$ of ray $j$ and $\delta m_i$ is the electron density perturbation in block $i$ (Figure 4.4). In this study, relative values of perturbations have been used, to facilitate the comparison of the obtained perturbation with the a priori ionospheric model. For relative values, equation (4.13) becomes

$$ \frac{\delta T_j}{T_j^0} = - \frac{40.3}{f_{ej}^2 c T_j^0} \sum_{i=1}^{N} m_i^0 \frac{\delta m_i}{m_i^0} ds_{ij}. \qquad (4.14) $$

$T_j^0$ ($T_j^0 = T_{phase}^{synth}$) is the synthetic propagation time in the a priori ionospheric model $m_i^0$ and $i$ is the block index. Introducing a matrix $\mathbf{A}$ of size $M \times N$ with $M$ being the number of travel time measurements and $N$ the number of basic functions (that is the number of blocks in the parametrization of the ionospheric model), equation (4.14) can be rewritten as

$$ \frac{\delta T_j}{T_j^0} = \sum_{i=1}^{N} \frac{\delta m_i}{m_i^0} A_{ji} \qquad (4.15) $$
or in tensor notation as

$$ \frac{\delta \mathbf{T}}{\mathbf{T}^0} = \mathbf{A} \cdot \frac{\delta \mathbf{m}}{\mathbf{m}^0}. \qquad (4.16) $$

$\delta \mathbf{m}$ is the model vector with $N$ unknown electron density perturbations $\delta m_i$. The vector $\delta \mathbf{T}$ contains $M$ observed travel time perturbations $\delta T_j$, and $\mathbf{A}$ is the geometric matrix containing $M \times N$ ray-path segments $ds_{ij}$ of ray $j$ in block $i$ multiplied by a coefficient

$$ A_{ji} = - \frac{40.3}{f_{ej}^2 c T_j^0} m_i^0 ds_{ij}. \qquad (4.17) $$

In order to find relative perturbations $\frac{\delta \mathbf{m}}{\mathbf{m}^0}$, the inverse of the matrix $\mathbf{A}$ has to be calculated according to the equation

$$ \frac{\delta \mathbf{m}}{\mathbf{m}^0} = \mathbf{A}^{-1} \cdot \frac{\delta \mathbf{T}}{\mathbf{T}^0}. \qquad (4.18) $$

Since there is no exact solution of this problem, the best thing to do is to estimate that model parameters giving the best approximate solution. According to Menke (1989), the best approximate solution $\frac{\delta \mathbf{m}}{\mathbf{m}^0}$ of an inverse problem has to minimize the $L_2$ norm of data and model, nominally

$$ \left\| \frac{\delta \mathbf{T}}{\mathbf{T}^0} - \mathbf{A} \cdot \frac{\delta \mathbf{m}}{\mathbf{m}^0} \right\|^2 = \min \qquad \left\| \frac{\delta \mathbf{m}}{\mathbf{m}^0} \right\|^2 = \min. \qquad (4.19) $$

The first condition imposes the best fit to the data, the second one minimizes the discrepancy from the a priori model. The Euclidean length of the solution quantifies the simplicity of the solution. A solution is defined to be simple if it is small when measured under the $L_2$-norm. This is a prior assumption to the solution of the problem. Then, in least-squares sense the solution of equation (4.18) is

$$ \frac{\delta \mathbf{m}}{\mathbf{m}^0} = (\mathbf{A}^T \cdot \mathbf{A})^{-1} \cdot \mathbf{A}^T \cdot \frac{\delta \mathbf{T}}{\mathbf{T}^0} \Leftrightarrow \left\| \frac{\delta \mathbf{T}}{\mathbf{T}^0} - \mathbf{A} \cdot \frac{\delta \mathbf{m}}{\mathbf{m}^0} \right\|^2 = \min. $$

$$ \frac{\delta \mathbf{m}}{\mathbf{m}^0} = \mathbf{A}^T \cdot (\mathbf{A} \cdot \mathbf{A}^T)^{-1} \cdot \frac{\delta \mathbf{T}}{\mathbf{T}^0} \Leftrightarrow \left\| \frac{\delta \mathbf{m}}{\mathbf{m}^0} \right\|^2 = \min. $$

The minimum norm solution in the second line can be obtained using Lagrange multipliers and its derivation is described in detail by (Menke, 1989, Chapter 3).

It is known that inverse problems in geophysics generally present a number $N$ of parameters (here the vector $\delta \mathbf{m}$) larger than the number $M$ of observations (here the vector $\delta \mathbf{T}$). These problems are called under-determined inverse problems. In practice, inverse problems are never purely under-determined, but mixed-determined,
meaning that some model parameters are over-determined and others are under-determined. For instance, in this tomography problem there may be cells passed through by several rays and others which are missed entirely (Figure 4.4). Consequently, the matrices $A^T \cdot A$ and $A \cdot A^T$ cannot be inverted, because they are almost singular. Even in the case of formally existing inverse matrices, these are often ill-conditioned that means small changes in the data vector ($\delta T$) lead to large changes in the model estimation ($\delta m$).

To find a stable solution ($\frac{\delta m}{m^0}$) balancing the sensitivity to the data as well as coherence with the a priori model, Menke (1989) suggests a damped least squares solution

$$ \frac{\delta m}{m^0} = (A^T \cdot A + \lambda \cdot I)^{-1} \cdot A^T \cdot \frac{\delta T}{T^0}, \quad (4.20) $$

minimizing the cost function

$$ F = \left\| \frac{\delta T}{T^0} - A \cdot \frac{\delta m}{m^0} \right\|^2 + \lambda \left\| \frac{\delta m}{m^0} \right\|^2. \quad (4.21) $$

$\lambda$ is a regularization parameter and $I$ is the identity matrix.

Simultaneous minimization of both terms is not possible, but the parameter $\lambda$ controls the emphasis that is put on the conflicting requirement. If $\lambda$ is set to zero, the misfit to the data is minimized, but no prior assumption to the solution is made to rule out the under-determined model parameters. The choice of this regularization parameter $\lambda$ is a very sensitive issue and will be explained in Section 4.6.3.

### 4.5.3 Damped least squares inversion with constraints

The damped least squares inversion described before imposes the best fit to the data and minimizes the discrepancy from the a priori model, given by equations 4.19. However, this inverse problem is an under-determined problem, because the number of rays is smaller than the number of pixels ($M < N$), e.g., 1056 rays and 2000 cells. In particular, the number of cells not hit by any ray ~ 1600 is large as compared to the number of cells crossed by rays ~ 400. As a result, the ray geometry defines a singular matrix, which makes the reconstruction of the electron density in the probed media difficult.

For those pixels without any traversing rays them, or cells with just a few number of rays, additionally constraints that impose smoothness to the solution can be included in the algorithm. These constraints should minimize the difference between the solution in one cell and its eight (or less in the case of border cells) adjacent
cells (Figure 4.5). For each of this eight constraints a matrix can be constructed.

<figure><img src="image_33.png" /><figcaption>Figure 4.5: Sketch for creating eight constraint matrices minimizing the difference of the solution in cell m,n and its eight adjacent cells.</figcaption></figure>

For instance, a constraint matrix minimizing the difference between the solutions in cell m, n and cell m + 1, n can be defined as

$$
\mathbf{R}_1 = \begin{pmatrix}
1 & -1 & 0 & 0 & 0 & \cdots & 0 \\
0 & 1 & -1 & 0 & 0 & \cdots & 0 \\
0 & 0 & 1 & -1 & 0 & \cdots & 0 \\
\vdots & \vdots & 0 & 1 & -1 & \cdots & 0 \\
\vdots & \vdots & \vdots & \vdots & \ddots & \ddots & 0 \\
0 & 0 & 0 & 0 & 0 & -1 & 1
\end{pmatrix} .
$$

Attention has to be paid at the boundaries in altitude and distance and the arrange-
ment of the grid. Similarly, seven additional matrices **R** can be created with each
matrix being different. Then, the cost function (equation 4.21) becomes

$$
F = \underbrace{\left\| \frac{\delta \mathbf{T}}{\mathbf{T}^0} - \mathbf{A} \cdot \frac{\delta \mathbf{m}}{\mathbf{m}^0} \right\|^2}_{\text{data misfit}} + \underbrace{\lambda \left\| \frac{\delta \mathbf{m}}{\mathbf{m}^0} \right\|^2}_{\text{solution norm}} + \underbrace{\lambda_1 \left\| \mathbf{R}_1 \frac{\delta \mathbf{m}}{\mathbf{m}^0} \right\|^2 + \cdots + \lambda_8 \left\| \mathbf{R}_8 \frac{\delta \mathbf{m}}{\mathbf{m}^0} \right\|^2}_{\text{difference to the 8 adjacent cells}} . \quad (4.22)
$$

Minimizing the cost function (that means $\frac{\partial F}{\partial m} = 0$) yields

$$
\frac{\partial F}{\partial m} = 0 = 2\mathbf{A}^T \mathbf{A} \frac{\delta \mathbf{m}}{\mathbf{m}^0} - 2\mathbf{A}^T \frac{\delta \mathbf{T}}{\mathbf{T}^0} + 2\lambda \frac{\delta \mathbf{m}}{\mathbf{m}^0} + 2\lambda_1 \mathbf{R}_1^T \mathbf{R}_1 \frac{\delta \mathbf{m}}{\mathbf{m}^0} + \dots + 2\lambda_8 \mathbf{R}_8^T \mathbf{R}_8 \frac{\delta \mathbf{m}}{\mathbf{m}^0}. \quad (4.23)
$$

The best solution of this inverse problem is

$$
\frac{\delta \mathbf{m}}{\mathbf{m}^0} = (\mathbf{A}^T \mathbf{A} + \lambda \mathbf{I} + \lambda_1 [\mathbf{R}_1^T \mathbf{R}_1 + \cdots \mathbf{R}_8^T \mathbf{R}_8])^{-1} \cdot \mathbf{A}^T \cdot \frac{\delta \mathbf{T}}{\mathbf{T}^0}, \quad (4.24)
$$

assuming that $\lambda_1 = \lambda_2 = \cdots = \lambda_8$. There are two regularization parameters to choose, namely $\lambda$ and $\lambda_1$.
### 4.5.4 Consideration of the ray-path deflection

In monostatic OTH radars, the endpoints of the rays (where the signal is backscattered by the ground) are unknown, so the location of the scattering point at the ground can change for a constant elevation angle and depends on the electron density variation $\delta N_e(\vec{r})$ in the ionosphere (Figure 4.4). Thus, the ray-paths in the real ionosphere and in the ionospheric model are not equal, $s(n) \neq s(n_0)$, introducing an additional shift $\delta T_{phase}^{ray}$ in the propagation time. Snieder and Spencer (1993) showed that both effects (velocity variation and ray-path deflection) can be combined in a single perturbation theory. At first-order, the two effects are simply additive (Snieder and Spencer, 1993) and the perturbation in the propagation time can be expressed as

$$ \delta T_{phase}^{real} = \delta T_{phase}^{velocity} + \delta T_{phase}^{ray}. \quad (4.25) $$

The deduction of an equation for the propagation time that takes into account both effects, the velocity variation and the change of the ray-path, requires that the integrals in the second line of equation (4.9) are not neglected. Replacing $\delta s$ in equation 4.9 with $\delta s = s' - s_0$, gives the following term for the propagation time in the real ionosphere (Roy et al., 2014)

$$ \begin{aligned} T_{phase}^{real} = & \frac{1}{c} \int_{s'} ds - \frac{1}{c} \int_{s_0} ds \\ & - \frac{40.3}{c f_e^2} \left[ \int_s (N_e^0(\vec{r}) + \delta N_e(\vec{r})) ds - \int_{s_0} (N_e^0(\vec{r}) + \delta N_e(\vec{r})) ds \right] (4.26) \\ & + \frac{1}{c} \int_{s_0} ds - \frac{40.3}{c f_e^2} \int_{s_0} (N_e^0(\vec{r}) + \delta N_e(\vec{r})) ds. \end{aligned} $$

That is

$$ T_{phase}^{real} = \frac{1}{c} \int_{s'} ds - \frac{40.3}{c f_e^2} \int_{s'} (N_e^0(\vec{r}) + \delta N_e(\vec{r})) ds. \quad (4.27) $$

This equation describes both effects, velocity variation and ray-path deflection, on the propagation time of EM waves in the real ionosphere. It is obvious that this equation can be simply derived from equation 4.7 by replacing $N_e(\vec{r})$ with $N_e^0(\vec{r}) + \delta N_e(\vec{r})$ and $s(n)$ by $s'$. For the difference in the propagation time $\delta T_{phase}^{real}$ between real (eq. 4.27) and synthetic propagation time (eq. 4.10) in the a priori ionospheric
model now follows

$$
\begin{align}
\delta T_{\text{phase}}^{\text{real}} &= T_{\text{phase}}^{\text{real}} - T_{\text{phase}}^{\text{synth}} \nonumber \\
&= \frac{1}{c} \int_{s'} ds - \frac{1}{c} \int_{s_0} ds \nonumber \\
&\quad - \frac{40.3}{c f_e^2} \int_{s'} (N_e^0(\vec{r}) + \delta N_e(\vec{r})) ds + \frac{40.3}{c f_e^2} \int_{s_0} N_e^0(\vec{r}) ds. \tag{4.28}
\end{align}
$$

This form of equation 4.25 allows to describe the difference between data and syn-
thetics accounting for both effects (velocity and ray-path deflection). Consequently,
this approach is called the **v&r-method**. Now it follows that the two terms of
equation 4.25, $\delta T_{phase}^{velocity}$ and $\delta T_{phase}^{ray}$, are given as

$$
\begin{equation}
\begin{split}
\delta T_{\text{phase}}^{\text{ray}} = {}& \frac{1}{c} \int_{s'} ds - \frac{1}{c} \int_{s_0} ds \\
& - \frac{40.3}{c f_e^2} \left[ \int_{s'} (N_e^0(\vec{r}) + \delta N_e(\vec{r})) \, ds - \int_{s_0} (N_e^0(\vec{r}) + \delta N_e(\vec{r})) \, ds \right]
\end{split}
\tag{4.29}
\end{equation}
$$

and

$$
\delta T_{\text{phase}}^{\text{velocity}} = - \frac{40.3}{c f_e^2} \int_{s_0} \delta N_e(\vec{r}) \, ds. \quad (4.30)
$$

$\delta T_{phase}^{ray}$ depends on the unknown electron density perturbation $\delta N_e$, that is the solution of the inverse problem mentioned above. To set up an inverse problem, which allow the determination of $\delta N_e$ based on observations of $\delta T_{phase}^{real}$, $\delta T_{phase}^{ray}$ needs to be expressed in terms of known $\delta N_e$. The study of *Snieder and Spencer (1993)* also indicates that sensitivity kernels $K(s)$ can be defined so that

$$
\delta T_{\text{phase}}^{\text{ray}} = \int_{s_0} K(\vec{r}) \cdot \delta N_e(\vec{r}) \, ds, \quad (4.31)
$$

where $K(\vec{r})$ is the data kernel, which in the *v*-method described above is just a delta
function along the unperturbed ray-path $s_0$. Here, the kernel contains the Fréchet
derivatives $\partial T / \partial m$, where $\partial T$ is a perturbation in the propagation time caused
by a perturbation in the model *m*. If the relation between the model *m* and the
propagation time *T* is linear, the sensitivity function can be computed numerically
*K*(\vec{r}) by

1. assuming a Dirac delta perturbation $\delta N_e = \delta(s - s_i)$, where $i$ denotes the cell,

2. running the Tracé De Rayons (TDR) code mentioned before to obtain the travel time delay due to velocity variation, $\delta T_j^{\text{velocity}}$ for ray j (i.e., equation 4.30) and the total time delay $\delta T_{\text{phase}}^j$, (equation 4.28)
3. calculating $\delta T_{ray}$ by means of equation 4.25.

The resulting value can be substituted on the left-hand side of equation 4.31, while the right-hand side collapses to $K(s_i)$ by the properties of the Dirac delta. Iterating the procedure over the location of all cells $s_i$ allows to compute $K$ throughout the volume of interest. In practice, the Dirac delta is replaced by a localized perturbation of an arbitrary but finite known value $\delta N_e^*$ and $K$ is obtained by dividing the resulting value of $\delta T_{phase}^{ray}$ by $\delta N_e^*$, so the kernel is given by

$$K_{ji} = \frac{\delta T_j^*}{\delta N e_i^*}.$$

Once kernels have been computed, equation 4.25 can be replaced by

$$\delta T_{phase,j}^{real} = -\frac{40.3}{f_{ej}^2 c} \sum_{i=1}^{N} \delta m_i ds_{ij} + \sum_{i=1}^{N} K_{ji} \delta m_i \quad (4.32)$$

or in a tensor formalism as

$$\delta \mathbf{T} = (\mathbf{A} + \mathbf{A}') \cdot \delta \mathbf{m}, \quad (4.33)$$

where (for absolute values)

$$A_{ji} = -\frac{40.3}{f_{ej}^2 c} ds_{ij} \quad (4.34)$$

and

$$\mathbf{A}' = K_{ji} = \frac{\delta T_{i,j}^*}{\delta N e_i^*} \quad (4.35)$$

The damped least squares inversion described in Section 4.5.2 is applied to the matrix $\mathbf{A}$ in the case of the v-method, and to the matrix $\mathbf{A} + \mathbf{A}' = \mathbf{M}$ in the case of the v&r-method.

In order to proof the linearity of the problem and to choose the amplitude of the localized, imposed electron density perturbation $\delta N e_i^*$ in the calculation of the ray-path deflection, the ratio $\frac{\delta T^*}{dm^*}$ was traced as a function of the imposed perturbation $dm^*$ for one ray (i.e., 11 MHz, 35° take-off angle) perturbing all crossed cells one after another. Figure 4.6 shows the linear relationship between the propagation time difference of the perturbed and the not-perturbed ray $\delta T^*$ and the imposed perturbation $dm^*$. Furthermore, the figure proves that the ratio $\frac{\delta T^*}{dm^*}$ is independent of the imposed electron density $\delta N e_i^*$ for most of the cells. Only for cells at altitude around 80 km the ratio varies slightly due to the strong variability of the a priori ionospheric model. In the following, the imposed electron density $\delta N e_i^*$ was set to 1%.
<figure><img src="image_34.png" /><figcaption>Figure 4.6: Proof of the linearity of the inverse problem. The propagation time difference of one ray traced in the *a priori* ionospheric model and one ray perturbed by a localized, imposed electron density of abritray amplitude $dm^*$ (left) and the ratio $\frac{\delta T}{dm^*}$ as a function of imposed electron density perturbation $dm^*$. A column of points represents all cells crossed by the ray, and their color corresponds to the altitude.</figcaption></figure>

## 4.6 Inversion results for synthetics

In order to validate and compare the two methods described in the previous sections, the method was tested for synthetic data generated by ray tracing TDR in the continuous *a priori* ionospheric model NeQuick (Radicella and Leitinger, 2001). The data were simulated by tracing 1071 rays with elevation angles between 10° and 60° and in a frequency range of 6–16 MHz in the ionospheric model plus an additional, known perturbation $\delta Ne^{target}$. $\delta T$ is then the difference in propagation time between the ray tracing in the perturbed ionospheric model, and the ray tracing in the *a priori* ionospheric model, i.e., equation (4.11). The solution of the inverse problem, the obtained perturbation $\delta m$, has to correspond to $\delta Ne^{target}$. In the following sections, results obtained for a checkerboard perturbation (Lévêque et al., 1993) and for a localized perturbation of 0.1 % will be presented and discussed. The background ionosphere $Ne^{apriori}$ was generated for October at 12:00 UT with a solar flux of 198.1 SFU.

The ray-path deflection in the $v&r$-method was calculated by adding 1% of electron density perturbation cell by cell, since the ratio of travel time perturbation to imposed electron density perturbation $\frac{\delta T}{dm^*}$ is independent of $dm^*$ (see Section 4.5.4).

### 4.6.1 Effect of ray-path deflection

The theoretical limit of the developed method was investigated by calculation of a vector of traveltime perturbations $\delta T^{frozen}$ satisfying exactly the hypothesis that the electron density perturbation $\delta Ne^{target}$ modifies only the velocity of EM waves,
i.e., rays are frozen in the a priori model configuration. This theoretical data set represents the idealized case of no ray-path perturbation. This resembles to classic seismic tomography (where both endpoints of the ray are known). It should be inverted to explore the resolution linked to the data coverage, as well as the numerical noise introduced by the discretization of the models (number and dimensions of the cells).

<figure><img src="image_35.png" /><figcaption>Figure 4.7: Inversion results for a checkerboard benchmark test using analytically (left) and numerically (right) calculated frozen rays. The second plot from the top of each column shows the solution for the best regularization parameter, chosen as described in Section 4.6.3.</figcaption></figure>

Two tests were realized to explore the theoretical limit using numerical and analytical frozen rays. In the case of analytical frozen rays, the vector $\delta T^{frozen}$ of travel time perturbations was calculated analytically to calculate the propagation time, according to

$$ \left( \frac{\delta T}{T} \right)_{ana}^{frozen} = A_{s0} \cdot \left( \frac{\delta Ne}{N_e^0} \right)^{target}, \qquad (4.36) $$

where $A_{s0}$ contains the ray-path in the unperturbed model $m_0$ and is given by equation 4.17. Next using equation 4.20, $\delta T^{frozen}$ was inverted for the relative electron density perturbation $\delta Ne/N_e^0$. This is, of course, just a mathematical operation, but proves the functionality of the code and allows to estimate the noise level introduced by ray tracing and the discretization of the model by cells.

Numerical frozen rays were traced in the a priori ionospheric model and a perturbed propagation time along the unperturbed ray-path $s_0$ was calculated, using
the equation

$$T_{phase} = \frac{1}{c} \int_{s_0} ds - \frac{40.3}{f_{ej}^2 c} \int_{s_0} (N_e^0 + \delta N_e(\vec{r})) ds. \quad (4.37)$$

Then for the relative time delay $\frac{\delta T}{T}_{num}$ follows

$$\left(\frac{\delta T}{T_{phase}^{synth}}\right)_{num}^{frozen} = \frac{T_{phase} - T_{phase}^{synth}}{T_{phase}^{synth}} = -\frac{40.3}{f_{ej}^2 \cdot c \cdot T_{phase}^{synth}} \int_{s_0} \frac{\delta N_e(\vec{r})}{N_e^0} ds, \quad (4.38)$$

were $T_{phase}^{synth}$ is given by equation 4.10. The inversion is then performed by means of equation 4.20. Figures 4.7 and 4.8 show the inversion for three different regulariza-

<figure><img src="image_36.png" /><figcaption>Figure 4.8: Inversion results for a localized perturbation using analytically (left) and numerically (right) calculated frozen rays. The second plot from the top of each column shows the solution for the best regularization parameter, chosen as described in Section 4.6.3.</figcaption></figure>

tion parameters $\lambda$, using analytical and numerical frozen rays for a checkerboard and a localized perturbation. The inversion with numerically calculated frozen rays with the best regularization parameter (i.e., $\lambda = 6.3 \cdot 10^{-6}$ for the checkerboard and $\lambda = 5.6 \cdot 10^{-5}$ for the localized perturbation) reproduces well the target model, although some noise is present in the solution. At distances larger than 1500 km, the targets are not resolved due to the sparse ray coverage (Figure 4.9). Both tests show that the inversions do not depend significantly on the regularization parameter. This
result is also observed for the inversion using the analytically calculated travel time perturbation $\delta T_{ana}^{frozen}$, but the solution is less noisy and more damped. It reflects the noise introduced by tracing rays in a continuous ionospheric model like NeQuick instead of using a discretized model with the parametrization of homogenous blocks. The mean noise is around 12%.

<figure><img src="image_37.png" /><figcaption>Figure 4.9: Ratio of plasma to signal frequency along the rays, calculated using equation 4.40 in order to show the sensitivity of the rays to the medium. Rays are most sensitive to the medium where the ratio is approximately 1.</figcaption></figure>

The limited complexity of the solution can be explained by the varying sensitivity of the rays to the medium. The electron density in the ionosphere is strongly varying with altitude, from zero at around 80 km to $10^{11}-10^{12}e/m^3$ at around 300 km (Figure 2.4). Consequently, EM waves emitted by the radar are particularly sensitive to the zone where the rays are reflected that is where the plasma frequency $f_p$ approaches the emission frequency $f_e$. A more exact equation for the ratio of plasma to emission frequency at the reflection point can be derived from Snell's spherical law,

$$R \cdot n(0) \cdot \cos \phi_0 = (R + h) \cdot n(h) \cdot \cos \phi_h, \quad (4.39)$$

where R is the Earth radius, $n(0)$ and $n(h)$ the refractive index at the Earth surface and at the altitude $h$, respectively, $\phi_0$ the elevation angle of the ray, and $\phi_h$ the local elevation angle along the ray. Assuming that $n(0)$ is 1 and using equation 4.4 for
the refractive index, the ratio of plasma to emission frequency is given by

$$ \frac{f_p}{f_e \cdot \sqrt{1 - \frac{\cos^2 \phi_0 \cdot R^2}{(R+h)^2}}} = 1. \qquad (4.40) $$

Figure 4.9 shows this ratio along the traced rays color coded. The ratio equals 1 between 200 and 300 km altitude at around 500 km horizontal distance and at 1000 km distance. By contrast, the corner around 200 km of altitude and close to the radar (around 100 km distance) is not well reproduced. Indeed, this zone lacks of reflected rays and is only crossed by nearly vertically rays. Consequently, these regions are not well resolved in the inversion.

### 4.6.2 v-method vs. v&r-method

The inversion results $\frac{\delta m}{m}$ for the checkerboard benchmark test and a localized perturbation are illustrated in Figured 4.11 and 4.10, respectively, where in the solutions obtained with the v-method and the v&r-method are directly compared for three different regularization parameters $\lambda$. The model is parametrized by homogenous blocks of 25 km x 20 km, covering an area of 2500 km from the geographical coordinates of the radar Nostradamus and reaching up to 400 km altitude. As expected, the solution models depend significantly on the regularization parameters. This effect was less serious in the case of the idealized frozen ray (Figures 4.7 and 4.8).

With ray-path deflections (endpoint perturbations) taken into account in the inversion, resolution deteriorates and the results are affected more profoundly by the choice of $\lambda$ (eq. 4.20). The correct location (500 km distance, 200 km altitude) of maximum anomaly of the localized perturbation is only reconstructed by the v&r-method provided that an adequate value is assigned to $\lambda$; the v-method identifies a high $\delta Ne/Ne$ anomaly in the common area of 200–500 km in horizontal distance and 100–200 km altitude, but slightly mislocates it and does not reproduce its shape.

The better performance of the v&r- compared to the v-method is evident in Figure 4.12, where the v&r L-curve (Tikhonov, 1963) has a more pronounced corner than that from the v- method. Compared to the ideal frozen-ray case (Figure 4.8), both the v- and v&r-methods occasionally introduce large-scale negative $\delta Ne/Ne$ anomalies that do not correspond to any feature of the target model, where the sign of perturbation is always positive. These artifacts can be explained with nonlinear effects of ray-path deflections, which (endpoints not fixed) can, in principle, result in faster propagation times even if the velocity perturbation is negative. While the v&r-method allows identification of more details of the target in the case of the checkerboard perturbation (Figure 4.11), the v-method identifies only the negative
<figure><img src="image_38.png" /><figcaption>Figure 4.10: Inversion results for the *v*-method (left) and *v&r*-method (right) for a localized perturbation of 0.1 % using three different regularization parameters. The inversions for the best regularization parameter are shown in the second plot from top.</figcaption></figure>

and positive perturbation at 500 km distance, but cannot reproduce their shape and amplitude. This result is grid independent, as demonstrated in Figures 5.2 in Section 5, where the pixel dimensions were enhanced in distance and altitude, respectively. The plots show the inversion for the best value of regularization, chosen as described in the following Section. Independent of the grid, the *v&r*-method is able to reproduce more details of the perturbations, even for larger distances. Resolution of the tomography method for different targets and grid sizes will be discussed more detailed in Section 5. The bad reconstruction at the edges and for larger distances is due to the sparse ray coverage, as demonstrated in Section 4.6.1, this is also due to the sparse ray coverage, not to the ray-path deflection.

### 4.6.3 Dependence of the solution on regularization

As Roy et al. (2014) showed, the regularization parameter $\lambda$ strongly controls the noise level in the solution. This is also clearly visible in Figure 4.11. In Section 4.5.2 is discussed that the chosen regularization parameter $\lambda$ has to minimize the misfit to the data $dT$ as well as the solution norm quantifying the discrepancy with the *a priori* model (i.e., equation 4.19).

To explore the influence of $\lambda$, the inverse problem (eq. 4.20) is solved for a set of regularization parameters, and the misfit is analyzed as function of the normalized solu-
<figure><img src="image_39.png" /><figcaption>Figure 4.11: Inversion results for the *v*-method (left) and *v&r*-method (right) for a checkerboard perturbation of 0.1 % using different regularization parameters. The inversions for the best regularization parameter are shown in the second plot from the top.</figcaption></figure>

tion norm. Following Tikhonov (1963) the obtained plot (e.g., Figure 4.12) is called trade-off curve, also commonly known as *L*-curve, because of its shape. The best compromise between the minimum misfit and the minimum solution norm is usually obtained geometrically at the maximum curvature of the *L*-curve (Tikhonov, 1963). Examples for *L*-curves are given in Figure 4.12 for the localized and the checkerboard benchmark test using the *v&r* and the *v*-method. The best compromises are $\lambda_{best} = 0.5$ and $\lambda_{best} = 0.63$, respectively. Smaller regularization parameters for both tests and both inversion methods add noise in the inversion, so the solution norm is large, but the misfit is small. On the other hand, larger regularization parameters strongly damp the solution, therefore the solution norm is small, but the misfit is large.

It is important to notice that the best regularization parameter, chosen at the maximum curvature of the *L*-curve, strongly depends on the range of the explored parameters. Indeed, the minimum value $\lambda_{min}$ of the explored set of $\lambda$ strongly controls the resulting $\lambda_{best}$. This is clearly visible in Figure 4.13, showing the inversion results obtained for different ranges of regularization parameters ($\lambda = [\lambda_{min}, 5000]$), with $\lambda_{min} = [10^{-5}, 0.1]$) for the checkerboard perturbation. For each of these four ranges, a *L*-curve has been calculated, and the best value of regularization has been determined geometrically at the maximum curvature of the corresponding *L*-
<figure><img src="image_40.png" /><figcaption>Figure 4.12: *L*-curves for the localized perturbation (a) and the checkerboard benchmark test (b) and for using the *v*- and *v&r*-methods for inversion. Inversion results for the magenta points are shown in Figures 4.11 and 4.10. Diamonds represent the best regularization parameter. Misfit and solution norm are given by equation 4.19.</figcaption></figure>

curve as described before. Obviously, the choice of a too small value of $\lambda_{min}$ leads to a $\lambda_{best}$ value that imposes very much noise on the result of the inversion, e.g., $\lambda = 0.02$. At the opposite extreme, a larger value of $\lambda_{min}$ results in a $\lambda_{best}$ that reduces considerably the amplitude of the solved pattern, because the solution is damped too much.

In order to find the best range of regularization parameters for the calculation of the *L*-curve, an error curve has been calculated: For each $\lambda$-range [$\lambda_{min}, \lambda_{max}$], with $\lambda_{min} = [10^{-5}, 10^3]$ and $\lambda_{max} = 5000$, the *L*-curve has been calculated, and the best value of regularization has been chosen from this curve. Then, the error i.e., $\sum_{i=1}^{cells} \left( \left( \frac{\delta Ne}{Ne} \right)_i^{target} - \left( \frac{\delta Ne}{Ne} \right)_i \right)^2$, between the target model and the solution has been determined. The $\lambda_{min}$ minimizing the error curve is $\lambda_{min} = 10^{-3}$ for the checkerboard perturbations. The corresponding $\lambda_{best}$ is 0.5 and the inversion result using this regularization parameter is shown in Figures 4.11 in the second plot from top in the right column.
<figure><img src="image_41.png" /><figcaption>Figure 4.13: a) Inversion results for the checkerboard benchmark test using different ranges of regularization parameters for the calculation of the *L*-curve shown in b). b) *L*-curve for different *λ*-ranges for the checkerboard benchmark test. Diamonds represent the best value of regularization for each curve for which the inversion is shown in a). c) Corresponding error curve for the explored range of regularization parameters. Colored points represent the errors for inversions using the best regularization parameters from the *L*-curves. The red line is the sum over the target models, i.e., $\sum_{i=1}^{\text{cells}} (\frac{\delta N_e}{N_e})_i^2$. The errors converge to the red line when the solution is totally damped. Source: Roy et al. (2014)</figcaption></figure>

### 4.6.4 Damped least squares inversion with constraints - Results

In order to reduce the noise in the inversions showed by Roy et al. (2014) in particular at low altitudes around 100 km, the damped least squares inversion with constraints (Section 4.5.3) was applied to the *v&r*-method. This procedure is not expected to change radically the inversion result, but it smoothes the solution. As the previous inversions show, results obtained with the *v*-method are already smooth in the first
iteration.

The major problem is to choose the two regularization parameters for the mini-
mum solution norm, $\lambda$ and $\lambda_1$, for the constraint matrices. A useful starting point
for at least getting 'on the map' is using the trace of the matrices as regularization
parameters. That means, the regularization parameters are chosen as

$$
\lambda = \frac{\operatorname{trace}(\mathbf{M}^{\mathrm{T}} \cdot \mathbf{M})}{\operatorname{trace}(\mathbf{I})} \qquad \lambda_1 = \frac{\operatorname{trace}(\mathbf{M}^{\mathrm{T}} \cdot \mathbf{M})}{\operatorname{trace}(\mathbf{R})}. \qquad (4.41)
$$

I is the identity matrix, **M** denotes the matrix **A** + **A**' used in the *v&r*-method, and
**R** is

$$
\mathbf{R} = \mathbf{H}^{\mathrm{T}}\mathbf{H}, \tag{4.42}
$$

where **H** is the sum of the eight constraint matrices minimizing the difference to the
eight adjacent cells

$$
\mathbf{H} = \mathbf{R}_1 + \mathbf{R}_2 + \mathbf{R}_3 + \mathbf{R}_4 + \mathbf{R}_5 + \mathbf{R}_6 + \mathbf{R}_7 + \mathbf{R}_8. \quad (4.43)
$$

Unfortunately, in the present case, this choose of regularization means $\lambda$ is around

<figure><img src="image_42.png" /><figcaption>Figure 4.14: Errorcurves for the determination of the best combination of parameter $\alpha_0$ and $\alpha_1$ for a checkerboard perturbation and the localized perturbation for different daytimes. First, $\alpha_1$ was modified for a fixed $\alpha_0$. The red points mark the minimum errors, and the red lines are the sum over the cells of the target model, $\sum_i^{\text{cells}} (\frac{\delta N_e}{N_e^0})_i^2$.</figcaption></figure>
100 times larger than $\lambda_1$ and consequently, the weight given to the constraints would be negligible small. Therefore, the regularization parameters were fixed to the trace of the matrix as described in equation 4.41, but multiplied by two factors, $\alpha_0$ and $\alpha_1$. Thus, the best inversion result is given by

$$ \frac{\delta \mathbf{m}}{\mathbf{m}^0} = (\mathbf{A}^\mathrm{T} \mathbf{A} + \alpha_0 \lambda \mathbf{I} + \alpha_1 \lambda_1 \mathbf{R}) \cdot \mathbf{A}^\mathrm{T} \cdot \frac{\delta \mathbf{T}}{\mathbf{T}^0}. \quad (4.44) $$

The best combination of the parameters $\alpha$ are determined by calculation of an error curve keeping one parameter fixed and alternating the other. For each pair of these two parameters, the inversion was calculated with equation 4.44 and the difference between target and solution was plotted as a function of $\alpha$. The best combination of these two parameters is defined by the minimum difference to the target model. The explored range for the two parameters is 0 to 100 with a sampling interval of 2.

Figure 4.14 shows these error curves as function of the parameters $\alpha$ for a localized perturbation at different daytimes (8h, 12h, and 18h), and the checkerboard perturbation at 12h. For the localized perturbation at 12h and 8h, as well as for the checkerboard test at 12h, the error curves return the same value for $\alpha_0$ and slightly different values for $\alpha_1$. The difference in the second parameter is probably due to the sampling interval of 2. These results indicate the independence of the parameters $\alpha$ from perturbation and daytime.

However, at night, (18h), $\alpha_0$ is 2 as for the other daytimes, but $\alpha_1 = 2$ is much smaller. At night, rays are reflected to a lesser extent, due to low ionization, therefore, the tomographic matrix has changed dramatically as compared to 12h.

Figure 4.15 and Figure 4.16 compare the inversion results obtained with the damped least squares inversion without and with constraints for a checkerboard test and a localized perturbation at 12h. In the first case, the regularization parameter was chosen by the L-curve criterion as described in Section 4.6.3, and in the second case by the error curve for the $\alpha$ parameters in Figure 4.14. The effect of smoothing is clearly visible for both target models. The noisy solution at altitude 100 km disappears in the smoothed solution for the localized perturbation; and for the checkerboard perturbation, the shape of the negative perturbation at altitude 100 km and at the borders are better reconstructed.
<figure><img src="image_43.png" /><figcaption>Figure 4.15: Damped least squares inversion vs. damped least squares with constraints for a localized perturbation. From bottom to top: Target model, damped least squares inversion with the best regularization parameter $\lambda_{best}$ chosen from the L-curve, and damped least squares inversion with constraints.</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>Figure 4.16: Damped least squares inversion vs. damped least squares with constraints for a checkerboard perturbation. From bottom to top: Target model, damped least squares inversion with the best regularization parameter $\lambda_{best}$ chosen from the L-curve, and damped least squares inversion with constraints.</figcaption></figure>
4.6.5 Conjugate Gradient algorithm

The results published by Roy et al. (2014) that were described above, are based on the solution of equation 4.18 by a damped least squares inversion. To take into account more the nonlinearity of the problem, introduced by the ray-path deflec- tion, another method was applied, namely the Conjugate Gradient Method (CGM) (Hestenes and Stiefel, 1952). This method is an iterative algorithm for the numeri- cal solution of sparse systems of linear equations, like equation 4.16. In analogy to in the damped least square inversion, the minimum of a cost function F is searched, minimizing the error to the data and the discrepancy from the a priori model,

$$
\begin{equation}
\begin{aligned}
F(\mathbf{m}) = & \underbrace{(\mathbf{T} - \mathbf{T}^{\text{obs}})^T \cdot \mathbf{W}_{\text{d}} \cdot (\mathbf{T} - \mathbf{T}_{\text{obs}})}_{\text{data misfit}} \\
& + \underbrace{(\mathbf{m} - \mathbf{m}_{\text{prior}})^T \cdot \mathbf{W}_{\mathbf{m}} \cdot (\mathbf{m} - \mathbf{m}_{\text{prior}})}_{\text{model misfit}}.
\end{aligned}
\tag{4.45}
\end{equation}
$$

This equation a cost function similar to that in equation 4.21 except that the cost function is expressed on absolute values and the damped least squares solution was transformed into a more generally regularized least-squares solution. The latter avoids the choice of an optimum regularization parameter λ, a problem, which has no simple solution.

In general, the weight matrices **W**<sub>m</sub> and **W**<sub>d</sub> can be anything. In Bayesian statistics one treats the inverse problem from a statistical point of view combining a priori information about the data and the model with the data that are actually measured. The weight matrices reflect true physical *a priori* information for data and model. In such a Bayesian approach, the weight matrices are given by the covariance matrices for data and model,

$$
\mathbf{W}_{\mathrm{m}}=\mathbf{C}_{\mathrm{m}}^{-1} \quad \mathbf{W}_{\mathrm{d}}=\mathbf{C}_{\mathrm{d}}^{-1}
\qquad(4.46)
$$

The CGM is similar to the method of steepest descent, but with a small difference. The method of steepest descent starts at an arbitrary point, choosing the direction, in which *F* decreases most quickly, which is the direction opposite to the gradient of *F*. After a series of steps *x*<sub>(1)</sub>, *x*<sub>(2)</sub> ... the errors are sufficiently small and the solution is the desired one. This method often takes steps in the same direction as earlier steps (Figure 4.17), undoing success achieved previously. It would be better to move in none interfering directions, exhausting the error in one direction and ensure that no additional error is introduced in that orientation in the next step.

Therefore, the directions in the CGM are *conjugate* (perpendicular in some sense)
to all previous ones. The first iteration of CGM is identical to the first iteration of
Steepest Descent, meaning that one moves in the direction opposite to the gradient. The gradient of the cost function (equation 4.45) is given by (Tarantola, 2005, Chapter 6),

$$g_0 = C_M \cdot \frac{\partial F}{\partial \mathbf{m}} = C_M \cdot \mathbf{G} \cdot C_D^{-1} \cdot (\delta \mathbf{T} - \mathbf{A}, \delta \mathbf{m}) + \delta \mathbf{m}, (4.47)$$

where $\mathbf{G}$ is a matrix that contains the partial derivatives $\frac{\partial \mathbf{T}}{\partial \mathbf{m}}$.

However, an important and critical point in the CGM is the choice of a step length, which quantifies how far one moves in one direction. Generally, this is done by performing a line search that means one searches a step length that minimizes the cost function in the search direction. In this work, a modified algorithm of the conjugate gradient was applied, following Tape et al. (2007). Their work uses a polynomial line search. The CGM for $k$ iterations can be summarized by the

<figure><img src="image_45.png" /><figcaption>Figure 4.17: Conjugate gradient vs. Steepest Descent Method. Solid grey lines are contours of the quadratic form with constant value $F(x)$. Starting at the arbitrary point $x_0$, the conjugate gradient method (green) arrives at the solution $x^*$ after two iterations, the steepest descent method (purple) after five iterations. Source: http://visiblegeology.com/rowan/presentations/2013/MSc2PhD/#/4/2</figcaption></figure>

following steps.

1. For an initial arbitrary starting model $\mathbf{m}_0$, calculate the misfit $F(\mathbf{m}_0)$ (equation 4.45) and the steepest ascent vector $g_0$ (equation 4.47).

2. Calculate the search direction $\mathbf{p}$. In the first iteration the search direction is set
to minus the gradient. In subsequent iterations conjugate search directions are calculated based on the previous search directions, in order to avoid searching in directions that have been searched before and spoiling the minimization along another :

$$
\mathbf{p}_k = \begin{cases} -\mathbf{g}_k & k=1, \\ -\mathbf{g}_k + \beta_k \mathbf{p}_{k-1} & k>1 \end{cases}
$$

The parameter $\beta_k$ is given for instance by the Polak-Ribíre formula (Tarantola, 2005, Chapter 6).

$$
\beta_k = \frac{(\mathbf{g}_k - \mathbf{g}_{k-1}) \cdot \mathbf{C}_M^{-1} \cdot \mathbf{B}_0 \cdot \mathbf{g}_k}{\mathbf{g}_{k-1} \cdot \mathbf{C}_M^{-1} \cdot \mathbf{B}_0 \cdot \mathbf{g}_{k-1}}, \quad (4.48)
$$

where $\mathbf{B}_0$ is a preconditioning operator. The simplest choice for $\mathbf{B}_0$ is the identity matrix. There exist other formulations of the parameter $\beta$, too.

3. Perform a polynomial line search (Tape et al., 2007) to choose a step length $v_k$ that minimizes $F(\mathbf{m} + v_k \cdot \mathbf{p})$ that is the misfit in the search direction $\mathbf{p}$. The choice of this step size, which determines how far one should go in the search direction, is a key decision in this algorithm. To find a step size

(a) Calculate at first a trial step, or test parameter $v_t$ by interpolating the misfit function F with a quadratic polynomial

(b) Update the model in the search direction $\mathbf{p}$ to obtain a test model $\mathbf{m}_t = \mathbf{m} + v_t \cdot \mathbf{p}$ and calculate the misfit for this test model $F(\mathbf{m}_t)$

(c) Interpolate the misfit function F between the current misfit $F(\mathbf{m}_t)$ and the misfit for the test model $(\mathbf{m}_t)$ by a quadratic polynomial and choose the step size $v_k$ that gives the minimum of the polynomial.

4. Update the model $\mathbf{m}_{k+1} = \mathbf{m}_k + v_k \cdot \mathbf{p}$

5. If the misfit $F(\mathbf{m}_{k+1})$ is suitable small, $\mathbf{m}_{k+1}$ is the desired model, otherwise restart from 2.

4.6.6 Results

The above described CGM was applied to synthetic data generated by ray tracing
TDR in the continuous *a priori* ionospheric model NeQuick (Radicella and Leitinger,
2001) for October at 12 UT with a solar flux of 198.1 SFU. The data was simulated
by tracing 1071 rays with elevation angles between 10° and 60° and in a frequency
range of 6–16 MHz in the ionospheric model plus an additional, known perturbation
$\delta Ne^{target}$. Again, the solution of the inverse problem, the obtained perturbation $\delta m$, has to correspond to $\delta Ne^{target}$.

The CGM was first validated with the analytical frozen rays (Section 4.6.1), satisfying the hypothesis that the electron density perturbation $\delta Ne^{target}$ modifies only the velocity of EM waves, i.e., rays are *frozen* in the *a priori* model configuration. This theoretical data set represents the idealized case of no ray-path perturbation.

In order to choose the data and model covariance matrix, the idealized CGM was applied to a localized perturbation of 5%. The model variance signifies the allowed deviation of the solution from the *a priori* ionospheric model, the data variance signifies the precision in the data. A small model variance does not allow the solution to differ much from the background ionospheric model. The model covariance was modified from 1% to 10% and the data covariance was fixed to $8.3 \cdot 10^{-6}$s.

Figure 4.18 and 4.19 illustrate the evolution of data and model misfit as a function of iterations.

<figure><img src="image_46.png" /><figcaption>Figure 4.18: Model misfit as a function of iterations for different model covariances and a fixed data covariance of $8.3 \cdot 10^{-6}$s. The CGM was applied to a localized perturbation of 5% in the idealized case of frozen rays.</figcaption></figure>

Smaller model variances do not allow the solution to differ much from the *a priori* ionospheric model, therefore, data and model misfit are quickly converging for 1% and not much changing with iterations. The larger the model variance, the more the model can change, and the longer it takes to converge. After 16 iterations, the data misfit is smallest for the largest covariance.

The CGM method was now applied to two synthetic, known electron density
<figure><img src="image_47.png" /><figcaption>Figure 4.19: Model misfit as a function of iterations for different model covariances and a fixed data covariance of 8.3 · 10⁻⁶s. The CGM was applied to a localized perturbation of 5% in the idealized case of frozen rays.</figcaption></figure>

perturbations that are a localized perturbation and a checkerboard perturbation of 5%. The covariance matrix for the model and data, $C_M$ and $C_M$ are diagonal matrices with 10% model variance and 8.3·10⁻⁶s data variance, respectively. Figure 4.20 shows the evolution of the solution with iterations for the localized perturbation. Already after four iterations, the perturbation is at the right horizontal and vertical place and with further iterations, the surrounding noise is reduced. The data misfit (equation 4.45) is reduced with iterations and converges, as shown in Figure 4.21. The model misfit (equation 4.45) increases, because the solution is a difference to the *a priori* ionospheric model. After iteration seven, it is constant.

Similar results are obtained for the checkerboard perturbation (Figure 4.22). After six iterations all negative and positive perturbations up to 500 km are well recovered. In larger distances the perturbations are poorly reproduced due to the sparse ray coverage (Figure 4.9). For the checkerboard perturbation, the data misfit decreases and converges, the model misfit increases and converges to a constant value, too (Figure 4.23).

These preliminary tests show that the Conjugate Gradient method is working in the idealized case of the frozen rays. In a next step, this method has to be applied to synthetics considering that the ray-path in the *a priori* ionospheric model is not the same as in the perturbed ionosphere (Section 4.5.4).
<figure><img src="image_48.png" /><figcaption>Figure 4.20: Solution of the Conjugate Gradient Method in the idealized case of the frozen rays for a localized perturbation.</figcaption></figure>

<figure><img src="image_49.png" /><figcaption>Figure 4.21: Data and model misfit (equation 4.45) as a function of iterations for the Conjugate Gradient Method for a localized perturbation.</figcaption></figure>
<figure><img src="image_50.png" /><figcaption>Figure 4.22: Solution of the Conjugate Gradient Method in the idealized case of the frozen rays for a checkerboard perturbation.</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>Figure 4.23: Data and model misfit (equation 4.45) as a function of iterations for the Conjugate Gradient Method applied to a checkerboard perturbation.</figcaption></figure>
# Chapter 5

## Resolution tests

In order to evaluate the resolution of the tomography method, tests were performed for several checkerboard tests and localized perturbations of different size and amplitude. In this chapter, inversion results obtained for a selection of target models are discussed. This selection includes large scale perturbations of several hundred kilometers as well as very small perturbations of only few kilometers. Localized perturbations were created using gaussian functions and their size was estimated at the 3 dB level of the maximum amplitude.

### 5.1 Effect of cell dimension on inversion results

In order to gain insight into the possible resolution using this tomography method, a first test was performed with cell dimensions being increased in altitude for a fixed distance of 50 km and vice versa (for a fixed altitude of 20 km cell dimensions are enlarged in distance). All of these test sets use the same amplitude of electron density perturbation of 0.1%.

In the second test, the amplitude of the electron density perturbation was scaled up for a fixed grid, to find out the order of still resolvable perturbations. The *v* and the *v&r*-method were used for inversion of both test sets, and the best value of regularization was chosen as described in Section 4.6.3. In summary, these tests show the independence of the inversion results of the grid and demonstrate the possibility of this method to reconstruct small scale anomalies and separate them spatially.

#### 5.1.1 Localized perturbation

Two resolution tests of localized perturbations with dimension 470×71 km and 235×47 km were performed. The results are illustrated in Figure 5.1 and 5.2 respectively.
<figure><img src="image_52.png" /><figcaption>Figure 5.1: Inversion results for a localized perturbation of 0.1% using the *v*-method and the *v&r*-method. The cell dimensions are 50 km in distance and varied in altitude from 10 km to 25 km and 40 km in the top panel. In the bottom panel, cell dimensions are 20 km in altitude and varied from 25 km to 50 km and 100 km in distance.</figcaption></figure>

In both cases, the correct location (500 km distance, 180 km altitude) of maximum anomaly in the target model is only reconstructed by the *v&r*-method. The *v*-method identifies a high $\delta\text{Ne}/\text{Ne}$ anomaly in the general area of 200–500 km in
<figure><img src="image_53.png" /><figcaption>Figure 5.2: Inversion results for a localized perturbation of 0.1% using the *v* and the *v&r*-method. The cell dimensions are fixed to 50 km in distance and increased in altitude from 10 km to 25 km and 40 km in the top panel. In the bottom panel cell dimensions are fixed to 20 km in altitude and enlarged from 25 km to 50 km and 100 km in distance.</figcaption></figure>

horizontal distance but slightly mislocates it at 100–200 km in altitude and does not reproduce its shape. Both *v*- and *v&r*-method introduce large-scale low $\delta Ne/Ne$ anomalies that do not correspond to any feature of the target model, where the
sign of the perturbation is always positive. As discussed in section 4.6.2, these are
artifacts to the nonlinear effects of ray-path deflections, which (endpoints not being
fixed) can in principle result in a faster propagation time even if the velocity per-
turbation is negative. These are more clearly visible for larger cell dimensions. The
grid sizes does not influence the result.

**5.1.2 Checkerboard tests**

The inversion results for an unrealistic checkerboard test with alternately changed perturbation (positive and negative) are illustrated in Figure 5.3 and Figure 5.4, where the size of the perturbations is the same, but the second checkerboard test was shifted downwards by 60 km.

The *v*-method only recognizes a positive and negative perturbation at around 200 km and 300 km altitude at 500 km horizontal distance. The two negative perturbations at the top and at the right boundary cannot be separated by the *v*-method. The shape of the positive anomaly is not well reconstructed and mislocated. While *v&r*-method cannot separate the two negative perturbations at the top neither, more details of the perturbation, also in larger horizontal distances are recognized in the inversion. These resolution tests show that the grid size is not influencing the inversion results. Results are not worse or better with a smaller or larger grid size.

The test with the downwards shifted checkerboard test confirms that the *v&r*-method performs quite better than the *v*-method, because all perturbations are recognized, although not clearly separated by each other. In the *v*-method the missing of details in the inversion result is obvious, and only two perturbations at between 200 km and 300 km altitude are visible.

In the unrealistic case of the radical checkerboard test (Figure 5.5) the sign of
the perturbation is alternating from cell to cell. While the *v*-method can partly
reproduce the pattern of alternating positive and negative perturbation between
200 and 300 km altitude, the *v&r*-method finds large anomalies at 100 km altitude
that do not correspond to the pattern in the target and mask the perturbations at
larger altitudes.
<figure><img src="image_54.png" /><figcaption>Figure 5.3: Inversion results for a checkerboard perturbation of 0.1% using the *v*- and the *v&r*-method. The cell dimensions are fixed to 50 km in distance and enlarged in altitude from 10 km, to 25 km and 40 km in the top panel. In the bottom panel, cell dimensions are fixed to 20 km in altitude and enlarged from 25 km to 50 km and 100 km in distance.</figcaption></figure>
<figure><img src="image_55.png" /><figcaption>Figure 5.4: Inversion results for the checkerboard perturbation in Figure 5.3, but shifted downward at around 60 km. The *v*- and *v&r*-method were used for inversions. The cell dimensions are fixed to 50 km in distance and enlarged in altitude from 10 km to 25 km and 40 km in the top panel. In the bottom panel, cell dimensions are fixed to 20 km in altitude and enlarged from 25 km to 50 km and 100 km in distance.</figcaption></figure>
<figure><img src="image_56.png" /><figcaption>Figure 5.5: Inversion results for a checkerboard perturbation of 0.1% using the *v* and the *v&r*-method. The cell dimensions are fixed to 50 km in distance and enlarged in altitude from 10 km, to 25 km and 40 km in the top panel. In the bottom panel, cell dimensions are fixed to 20 km in altitude and enlarged from 25 km to 50 km and 100 km in distance. The sign of the perturbations is alternating from cell to cell (radically).</figcaption></figure>
5.1.3 Two localized perturbations (165 × 47 km)

The tomography method is also able to recognizes and separate smaller scale anomalies, e.g., Figure 5.6. These tests represents a realistic case of a depletion and an enhancement of electron density in the ionosphere. The two localized (positive and negative) perturbations are located at 300 and 500 km respectively, in 250 km altitude. The v&r-method identifies both perturbations and is able to separate them, independent of the grid. Both maxima anomalies are reproduced well. In the v-method the negative anomaly is always correctly located, but the positive localized perturbations is smeared and not well recovered.
<figure><img src="image_57.png" /><figcaption>Figure 5.6: Inversion results for two localized perturbation of 0.1% using the *v*- and the *v&r*-method. The cell dimensions are fixed to 50 km in distance and increased in altitude from 10 km, to 25 km and 40 km in the top panel. In the bottom panel, cell dimensions are fixed to 20 km in altitude and enlarged from 25 km to 50 km and 100 km in distance.</figcaption></figure>
## 5.1.4 Cell perturbations

Supposing the unrealistic case of a perturbation in only one cell (Figure 5.7), the *v*-method performs much better than the *v&r*-method, because only the *v*-method is able to reproduce the correct location of the perturbed cell, even for the smallest cells of 25 × 20 km in distance and altitude. By contrast to the checkerboard tests and the localized perturbations, in this tests, the target model changes when the dimension of the cells are enlarged.

The *v&r*-method does not find the cell perturbation in any case, even for the largest pixel of 100 km horizontal distance, and the inversion result is strongly damped and blurred over the model space, so that no maximum anomaly can be found.
<figure><img src="image_58.png" /><figcaption>Figure 5.7: Inversion results for a localized perturbation in one cell of 0.1% using the *v* and the *v&r*-method. The cell dimensions are fixed to 50 km in distance and increased in altitude from 10 km, to 25 km and 40 km in the top panel. In the bottom panel, cell dimensions are fixed to 20 km in altitude and enlarged from 25 km to 50 km and 100 km in distance. The size of perturbations is changing with the grid.</figcaption></figure>
## 5.2 Which amplitude of perturbation is still resolvable ?

In the second test set, cell dimensions are fixed to 25 km × 20 km and only the perturbation amplitude is increased. Figure 5.8 shows the inversion results for the $v$ and $v&r$-method for the same localized perturbation as in Figure 5.2. The results indicate that small perturbations of electron density down to $10^{-3}\%$ and up to 20% can still be detected by both inversion methods. Perturbation of $10^{-4}\%$ are too small to be resolved by either of the inversion methods. While these small perturbations are still somehow visible in the $v$-method, the $v&r$-method detects nothing.

Independent of the amplitude, the $v$-method misallocate the perturbation, and in both method these large scale negative perturbations discussed previously are visible. Both inversion methods cannot solve localized perturbation larger than 50%, neither the amplitude nor the shape of the perturbation are resolved.
<figure><img src="image_59.png" /><figcaption>Figure 5.8: Inversion results for a localized perturbation of 0.1% using the *v* and the *v&r*-method. The cell dimensions are fixed to 25 km × 20 km and the amplitude of the perturbation is increased.</figcaption></figure>

## 5.3 Conclusions

The resolved area in all tests is limited by the uneven ray coverage, corresponding to a single OTH radar. Due to the strong variation of the ionospheric background