The Possible Seismo-Ionospheric Perturbations Recorded by the China-Seismo-Electromagnetic Satellite

: Driven by the objective of earthquake disaster prevention and mitigation, China launched the Zhangheng mission to build a stereoscopic earthquake monitoring system from the lithosphere to space. This report brieﬂy presents the possible seismic ionospheric disturbances recorded by the ﬁrst probe of the Zhangheng mission, which is known as the China-Seismo-Electromagnetic Satellite (CSES). The routine data preprocessing and seismo-ionospheric information analysis methods are brieﬂy introduced. The possible seismo-ionospheric disturbances that appeared during the strong shallow earthquakes (with a magnitude over 7 and a depth shallower than 30 km) are analyzed by using CSES and other multi-source data. Investigating seismo-ionospheric mechanisms requires multidisciplinary knowledge involving geophysics, atmosphere/ionosphere physics, geo-chemistry/atmospheric chemistry, etc. We state that the results from the CSES scientiﬁc application center are preliminary, calling for international scientists to contribute to the seismo-ionospheric perturbation phenomena, which is one of the most challenging scientiﬁc problems.


Introduction
Strong earthquakes (EQs) are one of the most destructive natural hazards, claiming countless deaths and economic losses in human history. The diversity and scale of EQ activity are daunting; however, the earthquake monitoring capabilities of human society still lack development. Promoting earthquake disaster prevention and reduction capabilities is a common issue faced by all countries worldwide.
The abnormal geophysical and geochemical variation around the EQ preparation zone is interpreted as the seismic precursors [1][2][3]. Electromagnetic signals are regarded as the most promising candidate for short-term earthquake prediction [4][5][6]. With the development of space technology, by the early 1980s, scientists found abnormal electromagnetic emissions, plasma parameter irregularities, and energetic particle precipitations over seismic fault zones from electromagnetic satellite observations [3,[7][8][9]. For example, Larkina et al. [7] reported an abnormal electromagnetic field intensity enhancement over the seismic zones from the Intercosmos 19 satellite data. Pulinets and Ouzounov [10] systemically summarized previous studies and concluded that the seismic signals can propagate upward from the lithosphere to the atmosphere and ionosphere, in particular circumstances, even up to the inner magnetosphere.
These new scientific findings bring hope for a breakthrough in earthquake science. Several countries have started building dedicated electromagnetism satellites to monitor these earthquake precursors from space. In the early 21st century, France launched the DEMETER (Detection of Electromagnetic Emissions Transmitted from Earthquake Regions) satellite mission, which successfully operated from December 2004 to December 2010, for studying ionospheric perturbations caused by earthquakes, volcanic eruptions, and human activities [11]. Since then, increasing attention has been given to the seismo-ionospheric disturbances phenomena. Much evidence on possible seismo-ionospheric precursor signatures has been accumulated after extensive studies have been carried out.
Abnormal electromagnetic emissions in a broad frequency range were observed over the epicenter zones or their conjugated points during the impending earthquake time [11][12][13][14]. Nemec et al. [14] statistically evaluated the abnormal electromagnetic field before M 4.8+ shallow-depth EQs by using 3.5 years of DEMETER data and confirmed the existence of a very small but statistically significant decrease in the electromagnetic wave intensity before the mainshock time. The ionospheric parameter irregularity at different altitudes over the EQ preparation zones is most commonly reported [9,[15][16][17]. For example, Zhao et al. [15] reported an unusual large-scale enhancement of electron density during the 2008 great Wenchuan earthquake. Marchetti et al. [18] statistically analyzed the variation in the magnetic field and electron density recorded by Swarm satellites and reported anomalies starting about 40 days before the beginning of the seismic sequence, coinciding and following surface and atmospheric alterations. Energetic particle precipitation was also observed during the seismic time [11].
Driven by the objective of earthquake disaster prevention and mitigation, China has been devoted to building a stereoscopic earthquake monitoring system from ground to space, planning to launch a series of Low-Earth-Orbit (LEO) electromagnetism microsatellites in the following decades step by step [19]. This mission is called the Zhangheng mission, named after the ancient Chinese scientist Zhangheng who invented the world's first seismo-scope in the second-century CE.
The first electromagnetism microsatellite of the Zhangheng mission, known as the China Seismo-Electromagnetic Satellite (CSES or ZH-1), was successfully launched into a sun-synchronous circular orbit at an altitude of 507 km in the upper ionosphere [19]. TheCSES flies 15.2 orbits around the Earth per day at the local time around 02:00 a.m. (nightside) and 2:00 p.m. (dayside). It has a 5 days revisiting period for the same area. Until now, the CSES has been steadily operating in orbit for over four years. Its identical successor satellite will be launched to the same orbit space around February 2023. This paper briefly reports how the CSES scientific application center analyzes seismic activities using satellite data.

Satellite and Payloads
The CSES is specially designed to monitor the ionospheric perturbations associated with seismic activities in the lithosphere. It carries eight types of scientific payloads, including 13 individual instruments to acquire the global geomagnetic field, electromagnetic field/waves, plasma parameters, and energetic particles data. The onboard configuration was introduced by Shen et al. [19]. It also benefits space physics-related studies (e.g., space weather, ionospheric irregularities, and plasma environment). There are two working modes for payload operation: survey and burst mode. The former is operated along the whole orbit trajectory with lower sampling rates; the latter is only triggered above the global main seismic belts with a higher sampling rate. Here, we briefly introduce the scientific payloads as follows. Details of the scientific payloads are listed in Table 1. The geomagnetic field is detected by the high-precision magnetometer (HPM), which consists of two tri-axial fluxgate magnetometers (FGMs) [20] and one coupled dark-state magnetometer (CDSM) [21]. The FGMs provide the magnetic field vector data in the frequency range from DC (Direct Current) to 15 Hz, and the CDSM serves as a reference to the FGM by providing the scalar value of the total magnetic field. A tri-axial search coil magnetometer (SCM) detects the alternating magnetic field in the frequency range from 10 Hz to 20 kHz [22], providing data in the three detection frequency bands: ULF (Ultra-Low-Frequency, 10-200 Hz), ELF (Extremely-Low-Frequency, 200 Hz-2.2 kHz), and VLF (Very-Low-Frequency, 1.8 kHz-20 kHz). An electric field detector (EFD) measures the spatial electric field in the frequency band from DC to 3.5 MHz [23], providing data in the four detection frequency bands: ULF (DC-16 Hz), ELF (6 Hz-2.2 kHz), VLF (1.8 kHz-20 kHz), and HF (High-Frequency, 18 kHz-3.5 MHz). One advantage of the CSES satellite is that it provides high-resolution six-component electromagnetic field waveform data in ULF/ELF bands continuously along the whole orbit trajectory in the survey mode, allowing us to examine the wave propagation parameters at any time of interest. The waveform data of VLF/HF is downlinked only during burst mode due to the limitation of telemetry capability.
For the ionospheric plasma parameter detection, a Langmuir probe (LAP) and a set of plasma analyzer packages (PAPs) [24] provide the in situ ionospheric plasma parameters (electron and ion); a GNSS occultation receiver (GOR), and a tri-band beacon (TBB) transmitter together with about 15 ground-based receivers [25] provide the ionospheric profile parameters, such as the TEC (total electron content), electron density profile, and ionospheric scintillation index. For the energetic particles, there are a set of high-energy-particle packages (HEPPs) [26], and a high-energy-particle detector (HEPD) [27] provided by the Italian Space Agency, jointly detecting the energetic particles in an energy level range from 0.1 MeV to 200 MeV. The first six-month in-orbit commission tests demonstrate that the scientific instruments generally perform well, except for the PAP, EFD, and TBB, which have some defects in different aspects. PAP was contaminated after four months in orbit, leading to lower absolute values than expected; further evaluation suggests that the relative values of ion densities can be used in scientific applications. The EFD is heavily interfered with by the satellite-internal communication over the equatorial region and has a high noise level in the HF band. The middle frequency band (400 MHz) of the TBB (needs the ground-based receivers to realize the coherent beacon system) malfunctioned after launch, so the data quality of the TBB is still under evaluation. We suggest that extra caution is required when using the PAP, EFD, and TBB for scientific research.
The scientific data are packed into HDF (Hierarchical Data Format) with the descending (orbit path from north pole to south pole) and ascending (from south to north) half-orbit data files. There are five levels of data in total based on the working principles of each payload, which are described as follows: Level 0: The data are reconstructed after a series of preprocessing, including frame synchronization, de-randomization, decoding, and de-formatting on the raw telemetry data.
Level 1: The data obtained after general error elimination, physical unit conversion, format conversion, etc. based on the Level-0 data. Level 2: The calibrated physical values correspond to Level 1 data with orbit information after coordination system transformation and necessary data inversion.
Level 2A: Only the EFD in the ULF band and GOR deliver the Level 2A data. For the EFD, Level 2A is generated after eliminating the Vs × B effect in the ULF band, where vs. is the satellite's velocity and B is the geomagnetic field. For GOR data, it is the TEC values obtained after conversion with precision orbit determination information on Level 2.
Level-3: Time sequential data along satellite orbits generated after resampling, necessary spectral analysis based on Level-2 data.
Level-4: The global interpolation maps of physical values from Level-2 data. The Level 2 data products that can be directly used for scientific application are accessible via the website to the international scientific community (http://leos.ac.cn accessed on 26 December 2021).

Data Preprocessing Methods
The ionosphere is easily disturbed by various natural and artificial sources, such as the geomagnetic storms that mainly originate from the Sun or the inner magnetosphere, the intense lightning activity generated from the atmosphere, or the artificial VLF radio waves emitted by powerful terrestrial transmitters. Therefore, before searching for the seismic information from the satellite data, it is necessary to eliminate any possibilities from all other non-earthquake sources. Several data preprocessings are routinely operated to avoid non-EQ disturbance sources.
Step 01: Data cleaning, through which the orbit data affected by satellite control operations or the other artificial experiments are directly excluded, according to the health condition data of the platform and payloads.
Step 02: Space-weather condition check. Compared to the disturbance induced by the solar wind, geomagnetic storms, or lightning activities, the influence from earthquakes is relatively feeble. Under the disturbed space-weather conditions (Dst ≤ −30 nT or Kp ≥ 3), we usually do not analyze the seismo-ionospheric disturbances. That does not mean it is impossible to extract seismic signals under disturbed space-weather conditions. As the physical processes of the seismo-ionospheric phenomena are still poorly understood, selecting the quiet space weather condition can avoid many complicated factors. Additionally, as most previous studies have carried out, we mainly focus on the nighttime data to avoid the complicated ionospheric ionization effects in the daytime.
Step 03: The EQ influential area selection. According to the empirical equation of the earthquake preparation zone put forward by Dobrovolsky et al. [28], the influential zone of an earthquake is roughly computed by the equation R = 10 0.43M , where R is the diameter of the influential earthquake zone in the km scale, and M is the magnitude of the EQ. For convenience, considering the projection feature of a satellite orbit on the ground, we choose a square area, that is, the epicenter ± 5 • area for a magnitude 6 EQ and epicenter ± 10 • for a magnitude 7 EQ.
Step 04: Frequency band selection is essential for preprocessing electromagnetic data. The electromagnetic field or waves can be easily interrupted by parasitic signals or other interference sources, which cause uncertainties or errors. For example, below 200 Hz, the interferences originating from instruments or satellite platforms, such as electrostatic turbulence from solar panels or magnetic torque, must be carefully excluded [12]. Selecting a frequency range free of artificial interferences is essential before searching for seismicrelated emissions. According to the calibration tests in orbit, at the frequency of 27.5 Hz and its harmonics (from 2 to 5 times and 8-12 times), the SCM data are influenced by the magnetic torque that generates a magnetic moment to maintain the satellite platform position. Over the equatorial area, the electric field data in ULF/ELF/VLF bands are interfered with by the internal communication signals of the satellite platform.

Results
The CSES has been steadily operating in the topside ionosphere for over four years. As the CSES launched on 2 February 2018, there have been 38 strong shallow EQs (M 7+, depth shallower than 100 km) worldwide (http://www.ceic.ac.cn/ accessed on 26 December 2021) until 1 December 2021; the epicenter distribution is shown in Figure 1. Previous studies [6,11,14,29] have demonstrated that the seismo-ionospheric perturbation phenomena predominately appear during the shallow strong EQs, so we mainly focus on the strong EQs with a depth shallower than 30 km in this report. The possible seismic ionospheric disturbances are listed in Table 2, and the specific analysis methods are introduced in the following. All the listed information represents possible seismic anomalies; they still need more evidence before confirming them. We share these preliminary results with the public, calling for global efforts for further study with profound methods.

Single Orbit Analysis
The single orbit analysis is the first step to locate any perturbations recorded by satellites over the EQ influential zone. The physical values (i.e., the power spectral density or amplitude of the electromagnetic field, the density/temperature of plasma parameters, or the flux of energetic particles) recorded by those orbits passing over the EQ influential zone are processed separately as Level 3 data product. Figure 2 shows an example of the Level 3 data product routinely produced in the CSES scientific application center, which shows the electron density (Ne) and temperature (Te) recorded by orbit No. 13962_1 on 8 August, ten days before the 18 August 2020, M 7.0 Sumatra EQ (see Table 2). The red curves show the observations from the current Remote Sens. 2022, 14, 905 6 of 18 orbit; the black (average values) and blue (upper and lower boundary) curves present a background trend that is calculated by the revisiting orbits in the previous 30 days (6 orbits in total). It can be seen that the latest Ne data (red) recorded by the current orbit No.13962_1 abnormally increased during the 18 August 2020, M 7.0 Sumatra EQ, significantly exceeding out of the upper (blue) boundary. Considering the space weather conditions (see Figure 3), the enhancement of electron density on 8 August is possibly related to the space weather conditions based on the rough analysis from the Level 3 data product. It is noted that Te data need further validation, not recommended for seismic information analysis directly by using Level 3 data.  The magnetic field enhancement at the frequency of 155 Hz nearest the epicenter 7 and 3 days before the mainshock. The electron/ion disturbed 7, 6, 5, and 2 days before the mainshock.
Loyalty Is-29 August The electron density increased; the PSD values of the electromagnetic field at the ELF fre-   The electron density significantly increased 10 days before the mainshock. The electromagnetic field intensity in the ULF/ELF band increased on 4 days, and one day before the mainshock. The energetic particle flux from 100 to 200 keV increased 4 and 3 days before the mainshock.
The third panel shows the information of EQs that occurred along the current orbit trajectory within one month; the y-axis represents the direct distance between the epicenters and the satellite footpoint; the square, circle, and triangle denote the EQs that occurred before, during, and after the time of the current orbit, respectively; and the EQ magnitude indicates the size of the square, circle, and triangle. The red and green color gradient denote the specific time of EQs that occurred before and after the passing time of the satellite, respectively. The bottom-left figure shows the space weather index (Kp and Dst) within one month (before the date of the current orbit). The bottom-right global map shows the orbit project and the epicenter. Alternatively, the orbits passing through the epicenter areas are plotted as scatter values of different physical values along the orbit traces, as Figure 2 of Zhima et al. [13], which shows the abnormal enhancement of the power spectral density values in the frequency range from 300 to 800 Hz of the magnetic field. The third panel shows the information of EQs that occurred along the current orbit trajectory within one month; the y-axis represents the direct distance between the epicenters and the satellite footpoint; the square, circle, and triangle denote the EQs that occurred before, during, and after the time of the current orbit, respectively; and the EQ magnitude indicates the size of the square, circle, and triangle. The red and green color gradient denote the specific time of EQs that occurred before and after the passing time of the satellite, respectively. The bottom-left figure shows the space weather index (Kp and Dst) within one month (before the date of the current orbit). The bottom-right global map shows the orbit project and the epicenter.
Alternatively, the orbits passing through the epicenter areas are plotted as scatter values of different physical values along the orbit traces, as Figure 2 of Zhima et al. [13], which shows the abnormal enhancement of the power spectral density values in the frequency range from 300 to 800 Hz of the magnetic field.

Multi-Orbits Analysis
After the rough analysis using the Level 3 data products, the second step is to compare the relative long-term observations with the multi-orbits and their revisiting orbits. The physical values over the EQ influential zone from each orbit are extracted and rebuilt as time sequence data. The running quartile method [13,30] is applied to examine the longterm trend. The running medians and the interquartile ranges (i.e., IQR, the difference between the third and first quartiles) were computed using the previous six revisiting orbits (30 days) and the current one (7 orbits in total). The running median values of the time sequence are used to build the background trend, while the median values of the current orbit are regarded as the current observation. The upper and lower bounds are computed by the running median values of physical values ± n*IQR values; n is set from 1 to 5, depending on the different physical quantities. The disturbance amplitude (the difference values between the current observation and the background trend) is normalized as percentages, and if current observation values are within the upper and lower thresholds, the disturbance amplitude is then set to zero. Figure 3 shows an example of the revisiting orbit analysis, which shows the Ne and Te variation from 25 July to 24 August over the M 7.0 Sumatra epicenter (18 August 2020).

Multi-Orbits Analysis
After the rough analysis using the Level 3 data products, the second step is to compare the relative long-term observations with the multi-orbits and their revisiting orbits. The physical values over the EQ influential zone from each orbit are extracted and rebuilt as time sequence data. The running quartile method [13,30] is applied to examine the long-term trend. The running medians and the interquartile ranges (i.e., IQR, the difference between the third and first quartiles) were computed using the previous six revisiting orbits (30 days) and the current one (7 orbits in total). The running median values of the time sequence are used to build the background trend, while the median values of the current orbit are regarded as the current observation. The upper and lower bounds are computed by the running median values of physical values ± n*IQR values; n is set from 1 to 5, depending on the different physical quantities. The disturbance amplitude (the difference values between the current observation and the background trend) is normalized as percentages, and if current observation values are within the upper and lower thresholds, the disturbance amplitude is then set to zero. Figure 3 shows an example of the revisiting orbit analysis, which shows the Ne and Te variation from 25 July to 24 August over the M 7.0 Sumatra epicenter (18 August 2020). From top to bottom are the space weather index (Kp and Dst), the variation in Ne, the disturbance amplitude of Ne, the variation in Te, and the disturbance amplitude of Te. The physical values during the current orbit are plotted as blue curves, and the background median values are plotted as black curves; the upper and lower boundaries are plotted as pink curves. It can be seen that an abnormal increase in Ne appeared on 8 August (10 days); according to the Kp and Dst index, we suggest that it is mainly related to the space weather conditions. The slight increase of 13 August (5 days before the mainshock) appeared under quiet space weather conditions, suggesting that it is possibly linked to the EQ. Figure 4 shows another example of the electric field, which shows the PSD (power spectral density) values of the Ex, Ey, and Ez components in the frequency range of 1000 Hz ± 100 Hz during the 18 August 2020, M 7.0 Sumatra EQ, in the same format as Figure 3. It can be seen that there was a slight enhancement in the three components of the electric field under the quiet space weather conditions on 13 August, 6 days before the mainshock.
linked to the EQ. Figure 4 shows another example of the electric field, which shows the PSD (power spectral density) values of the Ex, Ey, and Ez components in the frequency range of 1000 Hz ± 100 Hz during the 18 August 2020, M 7.0 Sumatra EQ, in the same format as Figure  3. It can be seen that there was a slight enhancement in the three components of the electric field under the quiet space weather conditions on 13 August, 6 days before the mainshock. Alternatively, they are plotted in the format of Figure 5, which shows the ULF wave activities in the frequency band from DC to 15 Hz recorded by the revisiting orbits during two Alaska EQs, which occurred on 29 July 2021 (M 8.0) and August 14 (M 7.0) (see Table  2). It is noted that only orbits with abnormal ULF emissions and their revisiting orbits are plotted in Figure 5. Observations on July 28, and August 4 and 7 are removed as there is an apparent AE index enhancement due to the disturbed space weather condition. The variation in the magnetic field intensity along the latitude (−65° to 65°) obtained by the revisiting orbits shows that on 19 July (10 days) and 27 July (2 days) before the 29 July M Alternatively, they are plotted in the format of Figure 5, which shows the ULF wave activities in the frequency band from DC to 15 Hz recorded by the revisiting orbits during two Alaska EQs, which occurred on 29 July 2021 (M 8.0) and August 14 (M 7.0) (see Table 2). It is noted that only orbits with abnormal ULF emissions and their revisiting orbits are plotted in Figure 5. Observations on July 28, and August 4 and 7 are removed as there is an apparent AE index enhancement due to the disturbed space weather condition. The variation in the magnetic field intensity along the latitude (−65 • to 65 • ) obtained by the revisiting orbits shows that on 19 July (10 days) and 27 July (2 days) before the 29 July M 8.1 mainshock, the eastward magnetic field in a frequency range from DC to 15 Hz appeared with apparent disturbances. On 10 August (4 days), 12 Auguest (2 days), and 14 August (on the mainshock day) before the 14 August M 7.0 EQ, there were still similar abnormal ULF wave activities. On the other days (some revisiting orbits shown in Figure 4), the strange ULF wave activities were not discovered in the non-EQ activity times over the Alaska epicenter area. It is interesting that the similar but weak ULF emissions also occurred in the southern conjugate point of the epicenter, which is physically possible. The electromagnetic disturbances can propagate along the magnetic field line to the conjugate point of the counterpart hemisphere [10]. Thus, we suggest that the possible abnormal emissions in the ULF band were emitted during the Alaska EQs. However, more evidence is still needed to confirm whether they are linked to the seismic precursors. 4), the strange ULF wave activities were not discovered in the non-EQ activity times over the Alaska epicenter area. It is interesting that the similar but weak ULF emissions also occurred in the southern conjugate point of the epicenter, which is physically possible. The electromagnetic disturbances can propagate along the magnetic field line to the conjugate point of the counterpart hemisphere [10]. Thus, we suggest that the possible abnormal emissions in the ULF band were emitted during the Alaska EQs. However, more evidence is still needed to confirm whether they are linked to the seismic precursors.

Background Map
The longer-term observations (months or yearly data) under quiet space weather conditions are selected to build a background map. The global background map is built by the median or mean values for specific physical values routinely produced as Level 4 data products for every five days or one year or different seasons (dependent on the types of physical values). The global background map helps depict large-scale ionospheric structures or the electromagnetic field intensity distributions [31,32].
A residual map is built for the specific EQ influential area by the difference values between the background map and the current EQ days (every five days within one month). For example, Figure 6 shows the residual map about the magnetic field in the frequency range from 60 Hz to 70 Hz during the 29 July 2021 and 14 August, M 7.0 Alaska EQs. It is roughly seen from the residual maps that for these two EQs, the PSD values of the magnetic field became enhanced 10 to 5 days before the mainshock (see the map of 20 to 29 July, 4 to 13 August); then, they dropped to a lower level after the mainshock (see the maps of 30 July to 3 August and 14 to 18 August). Admittedly, the residual map always looks disturbed; it is difficult to obtain reliable seismic information most of the time, as shown in Figure 6.

Background Map
The longer-term observations (months or yearly data) under quiet space weather conditions are selected to build a background map. The global background map is built by the median or mean values for specific physical values routinely produced as Level 4 data products for every five days or one year or different seasons (dependent on the types of physical values). The global background map helps depict large-scale ionospheric structures or the electromagnetic field intensity distributions [31,32].
A residual map is built for the specific EQ influential area by the difference values between the background map and the current EQ days (every five days within one month). For example, Figure 6 shows the residual map about the magnetic field in the frequency range from 60 Hz to 70 Hz during the 29 July 2021 and 14 August, M 7.0 Alaska EQs. It is roughly seen from the residual maps that for these two EQs, the PSD values of the magnetic field became enhanced 10 to 5 days before the mainshock (see the map of 20 to 29 July, 4 to 13 August); then, they dropped to a lower level after the mainshock (see the maps of 30 July to 3 August and 14 to 18 August). Admittedly, the residual map always looks disturbed; it is difficult to obtain reliable seismic information most of the time, as shown in Figure 6.

Multi-Parameter Comparisons
The lithosphere-atmosphere-ionosphere is an open and coupling system; the electromagnetic emissions and the chemistry gases emitted by the stressed rocks in the EQ preparation zones lead to the variation in geophysical and geochemical field, thus directly impacting the ground-based instruments. The seismic signals can also couple into the atmosphere, breaking the chemical reaction balance in there, leading to the formation of local electric currents, air ionization, atmospheric gas condensation, and plasma parameter irregularity [10].
Therefore, in addition to CSES data, the CSES scientific application center also collects multi-source data. At present, the multi-source observations from infrared/hyperspectral remoting sensing satellites, ground-based electromagnetic field instruments, the vertical/oblique ionosonde, and ground-based GNSS receivers etc. can be collected for the first time to comparatively explore the seismic signals with the CSES data. Remote Sens. 2022, 14, x FOR PEER REVIEW 13 of 18

Multi-Parameter Comparisons
The lithosphere-atmosphere-ionosphere is an open and coupling system; the electromagnetic emissions and the chemistry gases emitted by the stressed rocks in the EQ preparation zones lead to the variation in geophysical and geochemical field, thus directly impacting the ground-based instruments. The seismic signals can also couple into the atmosphere, breaking the chemical reaction balance in there, leading to the formation of local electric currents, air ionization, atmospheric gas condensation, and plasma parameter irregularity [10].
Therefore, in addition to CSES data, the CSES scientific application center also collects multi-source data. At present, the multi-source observations from infrared/hyperspectral remoting sensing satellites, ground-based electromagnetic field instruments, the vertical/oblique ionosonde, and ground-based GNSS receivers etc. can be collected for the first time to comparatively explore the seismic signals with the CSES data.
We pay considerable attention to the atmospheric parameters because the gas released from the crust can effectively change the conductivity condition of the atmosphere, thus leading to a series of coupling processes in the lithosphere-atmosphere-ionosphere [10]. Previous studies have suggested that the atmospheric parameters, such as the radon, methane, aerosol optical depth (AOD), skin temperature (SKT), and surface latent heat flux (SLHF), are sensitive to seismic activities at different stages of the EQ preparation process and can respond sooner to the seismic precursors [33,34] than the seismo-ionospheric parameters. For example, Figure 7 shows the variation in the methane concentration from 2016 to 2021 (provided by the NASA Goddard Earth Sciences Data and Information Services Center) over the Alaska EQ area. Compared with the variation trend from We pay considerable attention to the atmospheric parameters because the gas released from the crust can effectively change the conductivity condition of the atmosphere, thus leading to a series of coupling processes in the lithosphere-atmosphere-ionosphere [10]. Previous studies have suggested that the atmospheric parameters, such as the radon, methane, aerosol optical depth (AOD), skin temperature (SKT), and surface latent heat flux (SLHF), are sensitive to seismic activities at different stages of the EQ preparation process and can respond sooner to the seismic precursors [33,34] than the seismo-ionospheric parameters. For example, Figure 7 shows the variation in the methane concentration from 2016 to 2021 (provided by the NASA Goddard Earth Sciences Data and Information Services Center) over the Alaska EQ area. Compared with the variation trend from 2016 to 2020, the CH4 concentration increased significantly on 11 July 2021 (18 days before the M 8.1 mainshock), then fell to a low level where the M 8.1 shock occurred; after that, it started to rise to a new peak where the second shock M 7.0 occurred, and the CHR concentration continually decreased and fell back to the background trend level after September 2021. Figure 8 shows the spatial distribution of the 8 days averaged AOD values (provided by AQUA MODIS C6 of NASA) in 2021 relative to the background field calculated by data from 2001 to 2020 over the Alaska earthquake area. The AOD values became enhanced from 26 June to 21 August, especially in the 8 days time windows of 26 June to 3 July, 28 July to 4 August, and 13 August to 21 August; the AOD is significantly increased relative to the background. The significant enhancements mainly appeared in the southern part of the earthquake epicenter. 2016 to 2020, the CH4 concentration increased significantly on 11 July 2021 (18 days before the M 8.1 mainshock), then fell to a low level where the M 8.1 shock occurred; after that, it started to rise to a new peak where the second shock M 7.0 occurred, and the CHR concentration continually decreased and fell back to the background trend level after September 2021.  July to 4 August, and 13 August to 21 August; the AOD is significantly increased relative to the background. The significant enhancements mainly appeared in the southern part of the earthquake epicenter.      The TEC data also easily became disturbed over the forthcoming epicenter regions, showing a significant decrease or increase within one week before the mainshock [35]. Figure 9 shows the vertical TEC anomaly over the Alaska EQs, extracted from the global ionospheric map (GIM) published daily by the International GNSS Service (IGS). The green curve presents the median values of TEC data over the Alaska epicenter zone, and the grey lines are the upper and lower boundaries computed by the IQR method in Section 4.2. The red curve shows the current TEC observations. The red and black bars are the difference values between the red and green curves. It can be seen that the continuous positive anomaly appeared before the July 29 M 8.1 EQ, followed by a negative anomaly after the mainshock. However, it is difficult to identify any apparent anomaly during the August 14, M 7.0 EQ. and the grey lines are the upper and lower boundaries computed by the IQR method in Section 4.2. The red curve shows the current TEC observations. The red and black bars are the difference values between the red and green curves. It can be seen that the continuous positive anomaly appeared before the July 29 M 8.1 EQ, followed by a negative anomaly after the mainshock. However, it is difficult to identify any apparent anomaly during the August 14, M 7.0 EQ.

Discussion
In addition to the single case studies, some statistical studies were also recently carried out using CSES data. For example, Li et al. [16] used the O+ density and electron density data from CSES and the DEMETER satellite to search for seismo-ionospheric perturbations with different time resolutions. The results showed that the detection rate of the seismo-ionospheric disturbance increases as the data time resolution and the earthquake magnitude increase and as the focal depth decreases, proving that the CSES could effectively register seismo-ionospheric perturbations as the DEMETER satellite does [16]. Zhu et al. [36] statistically examined the ~2.5 years of Ne data from CSES during the M 4.8+ EQs worldwide. They found that significant positive Ne variations related to earthquakes mainly occurred ~1 to 7 days and ~13 to 15 days before the earthquakes.
Although numerous independent studies have suggested the seismo-ionospheric disturbances, the earthquake preparation process in the lithosphere and the seismic precursor transport process from the ground to space are still poorly understood due to various variables involved. The pioneering scientists established the Lithosphere-Atmosphere-Ionosphere Coupling (LAIC) models [4,6,10] mainly based on theoretical analyses and observations. Zhao et al. [37] reported a good correlation between CSES's observation and the full-wave-mode simulation on the VLF radio waves emitted by the ground-based transmitters, confirming the propagation capability of the electromagnetic waves from the lithosphere to the ionosphere. However, due to the complicated seismogenic process and the limited understanding of the seismo-ionospheric phenomena, there are many challenges in the physical process interpretation.
Huang et al. [38] put forward three critical problems in this field: (1) how to perform a reliable and appropriate observation with some precise physical quantities; (2) how to develop a robust methodology to reveal weak earthquake-related electromagnetic signals

Discussion
In addition to the single case studies, some statistical studies were also recently carried out using CSES data. For example, Li et al. [16] used the O+ density and electron density data from CSES and the DEMETER satellite to search for seismo-ionospheric perturbations with different time resolutions. The results showed that the detection rate of the seismo-ionospheric disturbance increases as the data time resolution and the earthquake magnitude increase and as the focal depth decreases, proving that the CSES could effectively register seismo-ionospheric perturbations as the DEMETER satellite does [16]. Zhu et al. [36] statistically examined the~2.5 years of Ne data from CSES during the M 4.8+ EQs worldwide. They found that significant positive Ne variations related to earthquakes mainly occurred~1 to 7 days and~13 to 15 days before the earthquakes.
Although numerous independent studies have suggested the seismo-ionospheric disturbances, the earthquake preparation process in the lithosphere and the seismic precursor transport process from the ground to space are still poorly understood due to various variables involved. The pioneering scientists established the Lithosphere-Atmosphere-Ionosphere Coupling (LAIC) models [4,6,10] mainly based on theoretical analyses and observations. Zhao et al. [37] reported a good correlation between CSES's observation and the full-wave-mode simulation on the VLF radio waves emitted by the ground-based transmitters, confirming the propagation capability of the electromagnetic waves from the lithosphere to the ionosphere. However, due to the complicated seismogenic process and the limited understanding of the seismo-ionospheric phenomena, there are many challenges in the physical process interpretation.
Huang et al. [38] put forward three critical problems in this field: (1) how to perform a reliable and appropriate observation with some precise physical quantities; (2) how to develop a robust methodology to reveal weak earthquake-related electromagnetic signals from the noisy background; and (3) how to develop plausible physical models based on theoretical analyses and/or laboratory experiments for the explanation of the earthquakerelated electromagnetic signals observed under the field conditions.
To obtain convincing seismo-ionospheric evidence, we pay considerable attention to data cleaning and data quality control before seismic information extraction. We have performed specific validations on CSES data, and some works about the data quality of the electromagnetic field [37,39,40] and the plasma data [31,41] have been published.
It is admitted that the seismic information analysis methods are still empirical. At present, we interpret data using four methods: the single-orbit analysis, the multi-orbit data analysis, the background map, and the comparative studies with multi-source data recorded from the ground to space. These methods are not sufficient to extract accurate seismic signals. Objectively and correctly extracting the feeble seismic information from the noisy electromagnetic background is still a challenging task, needing a complex multidisciplinary approach to clarify the short-term earthquake precursors by combining the observations at the different altitudes from the lithosphere to the ionosphere [6].

Conclusions
This paper briefly introduces the CSES mission, an electromagnetism satellite dedicated to earthquake monitoring using space technology. The routine data preprocessing and seismo-ionospheric information analysis methods are briefly introduced. The possible seismo-ionospheric disturbances that appeared during the strong shallow earthquakes (with a magnitude over 7 and a depth shallower than 30 km) were analyzed by using CSES and other multi-source data.
Based on previous studies and the CSES observations, we suggest that monitoring the short-term precursors from a stereoscopic monitoring system built from the lithosphereionosphere is possible. However, it is admitted that there is still a long way to go to identify the seismo-ionospheric precursors from the electromagnetism satellites correctly. Due to the complexity of the EQ preparation mechanism and the limitation of observation technology and data analysis methods, there are many challenges, and we need a multidisciplinary perspective to explore this topic. Finally, we state that the obtained results from the CSES are preliminary, requiring global efforts to search for pre-earthquake short-term precursors using advanced space technology.