The investigated region (89°–105°E, 31°–43°N) encompasses a pivotal segment of the Tibetan Plateau’s northern tectonic transition belt, and as a leading-edge response area to the ongoing collision of the Indo-Eurasian plate, it is a typical area of intra-terrestrial tectonic deformation and strong seismic activity (Figure 1). The region is characterized by the rigid base units of the Qaidam Massif and the Bayan Har Massif, and is controlled by a number of active tectonic systems, such as the Qaidam North Rim Fracture System, the East Kunlun Fault zone, and the Yushu-Ganzi Fault zone, and is a key transition zone between uplift of the Tibetan Plateau and the extrusion of material to the east.

The North Qaidam Fault zone, acting as a progressive deformation front accommodating the plateau’s northward crustal extrusion, transitions northward into the Qilian orogenic belt along the Tibetan Plateau’s northeastern margin, adjacent to the more rigid Qaidam Massif on its southern side, and forms a tectonic convergence with the Arjin left-trending strike-slip fracture system on its western end. The surface cover of the fracture zone is dominated by the Gobi desert and salt lake deposits, and its low thermal inertia enhances the sensitivity of the tectonic stress-induced thermal-infrared anomaly response (Peng et al., 2012), and the seismic process of the 2009 Hercynian M6.4 earthquake is related to this special surface-tectonic coupling effect; As the main control fault zone at the northern boundary of the Tibetan Plateau, the East Kunlun Fault zone extends more than 1,500 km, and its western section is connected with the West Kunlun orogenic belt and its eastern section is adjacent to the Bayan Har Massif, which shows a sliding rate of 8–10 mm/yr (Zhao et al., 2022), and the 2001 M8.1 Kunlunshan earthquake generated a 426-km coseismic rupture zone, with a 96-km segment near Kunlun Pass highlighting strain partitioning dynamics of large-scale strike-slip systems (Chen et al., 2004); to the southeast, the Yushu-Ganzi Fault structurally links the East Kunlun and Xianshuihe tectonic systems while serving as a kinematic transition interface along the Bayan Har Block’s southern boundary, whose complex geometric barrier structure led to the localized accumulation of strains, which ultimately triggered the Yushu M7.1 earthquake in 2010.

Based on the China Earthquake Networks Center dataset (https://news.ceic.ac.cn/history.html), this study constructed an earthquake activity catalog for the period 2009–2024 (Table 1), 37 M ≥ 5.0 seismic events were cumulatively recorded in the study area, and their magnitude-frequency distributions conformed to the modified Gutenberg-Richter relationship (b = 0.89 ± 0.05) (Ren, 2012). It should be noted that some consecutive earthquakes in the same area within a short period of time (within 5 days), such as the 5.3-magnitude and 5.5-magnitude earthquakes that occurred in Zaduo County, Qinghai Province, China, on March 5 and 7 March 2024, respectively, had an interval of only 2 days between the two events. After comprehensive consideration, it was classified as a moderate-strength earthquake event, with the first occurrence time and magnitude used as the basis for recording. The seismic activity shows significant graded characteristics: seven M6.0-7.0 medium-strong earthquake events and two M ≥ 7.0 strong earthquake events. Among them, the 2021 Maduo M7.4 earthquake, as the largest example, was revealed to have generated a co-seismic surface rupture system of about 170 km along and near the East Kunlun Fault Zone through satellite image interpretation and field verification, resulting in a maximum surface radar line-of-sight deformation of about 2.2 m, with spatial spreading characteristics that showed a high degree of coupling with the geometry of the fracture zone (Wang et al., 2022).

The Outgoing Longwave Radiation (OLR) datasets analyzed in this study were retrieved from Fengyun geostationary meteorological satellites, operated by the CMA National Satellite Meteorological Center (https://satellite.nsmc.org.cn/DataPortal/cn/home/index.html). OLR represents the broadband electromagnetic energy radiated by Earth’s surface into space, commonly referred to as thermal radiative flux, with its spectral range spanning 4 μm to longer wavelengths, and the value mainly depends on the temperature of the emitted subsurface (Wu, 2007). Operational production of OLR products commenced in June 2005 through the National Satellite Meteorological Center (NSMC) satellite data processing system, utilizing multispectral infrared (IR1, IR2) and water vapor channel radiometric data from FY-2C and subsequent geostationary satellites to generate eight daily time steps of thermal emission measurements, and at the same time integrates the data into a time series and produces OLR data files with different temporal resolutions of days, 5-day period, 10-day period, and months. OLR data files and store them in formats such as HDF for archiving (Guo et al., 2019).

Fengyun (FY) geostationary satellite-derived OLR datasets provide extensive spatial coverage with high-frequency sampling capabilities, covering the period from January 2009 to December 2024. These observations integrate multi-platform measurements from FY-2C/E/G and FY-4A/B satellites, leveraging their unique orbital configurations for continuous Earth radiation monitoring. To ensure data path consistency, satellite orbital intervals with sub-satellite point longitudes approximating 105°E are prioritized. This selection criterion stabilizes radiation reception geometry across observation cycles. Meanwhile, considering the direct influence of solar radiation and the range of time zones in China, the data from 15:00–21:00 UTC are selected daily, and the ground gas radiation can be basically balanced in this time. Given the voluminous nature of the datasets, analytical workflows were spatially constrained to the latitudinal (5°N-50°N) and longitudinal (55°E−150°E) domains for computational efficiency, and the latitude and longitude steps were all set to 0.05° to establish the original database of geostationary satellite remote sensing OLR. Our technical flowchart is shown in Figure 2, which details the logical relationships between our research objects, methods, content, and research objectives. In the research method module, we sequentially employed key algorithms including wavelet transform, power spectrum estimation, and least squares estimation.

In this study, all M ≥ 5.0 seismic events in the northern Tibetan Plateau study area (89°-105°E, 31–43°N) during the period 2009–2024 were categorized into three magnitude bins: M5.0–6.0 (n = 28), M6.0–7.0 (n = 7), and M ≥ 7.0 (n = 2), and the efficacy of OLR power spectrum anomalies in the detection of the retrospectivity of the seismic events was evaluated separately.

Based on the OLR relative power spectrum analysis method proposed by Guo et al. (2014), seven OLR feature bands (bands 1–7) are extracted, whose corresponding periods are 64 days, 32 days, 21 days, 16 days, 13 days, 11 days, and 9 days. Combining with the consensus of the previous studies on thermal infrared seismic anomalies (Liu et al., 1999; Qiang et al., 1999), we rule out the duration of less than 15 days of transient noise signals, focusing on the persistent anomalies associated with the gestation seismic process.

Anomaly determination criteria integrate regional OLR data characteristics with tectonic activity patterns (Luo et al., 2009; Yang et al., 2020; Zhang et al., 2024; Guo et al., 2024):

1. Intensity threshold: relative power spectrum amplitude exceeding six times the annual average power spectrum amplitude;

Based on the above discriminative criteria, the OLR anomalies of all feature frequency bands in the study area are extracted, and the corresponding seismic events are restricted to be M ≥ 5.0 earthquakes occurring in the anomaly area and its neighboring 200 km (spatial range: 89°-105°E, 31–43°N). Table 2 demonstrates the retrospective test accuracies of all feature bands for different magnitude files. The statistical results show that for the M5.0–6.0 and M6.0-7.0 magnitude brackets, band 6 (i.e., with a characteristic period of 11 days) has the highest detection accuracy, 82.1% (23/28) and 100% (7/7), respectively; for the M ≥ 7.0 magnitude bracket, band 3 (i.e., with a characteristic period of 21 days) has the highest detection accuracy, 100% (2/2); for all magnitude brackets combined, band 6 has the best overall accuracy (81.1%) followed by band 3 (62.2%). Based on this, band 6 is preferred for the prediction models of M5.0–6.0 and M6.0–7.0 earthquakes, while band 3 is more appropriate for M ≥ 7.0 earthquakes, and subsequent studies and analyses for different magnitudes will be carried out based on the corresponding bands.

As shown in Figure 4, the study of the spatial and temporal evolution process of the OLR relative power spectrum anomaly in band 6 and its associated seismic events indicates that the anomaly evolution of the 2023 Jishishan M6.2 earthquake presents three stages of characteristics: ① Initial stage (21 September 2023): the OLR relative power spectrum anomaly is first visible in the Sunan-Qilian fault zone northwest of the epicenter; ② Development stage (September 22-October 23): the anomaly showed an eastward migration and expansion trend along the Sunan-Qilian fault zone, and reached its first peak on October 7, with the main body of the anomaly located in the northwestern part of the epicenter, and the area covered extended to 1.0 × 10⁵ km², with the maximum ratio reaching 16 times the regional annual average; after a short period of attenuation of 3 days, the anomaly showed a secondary peak on October 23, with the main body of the anomaly located in the northeastern part of the epicenter (104.5°-105°E, 37.5°-38°N), and its spatial distribution (about 0.7 × 10⁵ km²) and the anomaly intensity (about 14 times, Figure 5) were lower than the first peak, which is presumed to be related to the phased release of the precursor energy of the pregnant earthquake system; ③ Decay phase (October 23-October 30): the anomaly intensity and spatial extent of the synchronous rapid decay, until October 30 completely disappeared, the total duration of 40 days, forming a typical “double-peak” evolutionary pattern. It is noteworthy that the epicenter (102.79°E, 35.70°N) is located at the edge of the southwest gradient zone of the OLR power spectrum anomaly area (subpeak), whose observation aligns closely with established statistical patterns in fault zone systems, where geophysical anomalies exhibit spatiotemporal migration along tectonic corridors, and seismic epicenters predominantly localize within anomaly gradient zones (Zhang et al., 2024). Spatial analyses show that the anomaly migration paths are significantly consistent with the strike of the Sunan-Qilian fault zone (NW-SE direction), which further confirms that tectonic activities drive the OLR anomaly.

Regarding Figure 5, it should be noted that there was an abnormal peak (approximately 4 times higher than normal) in this area in August. However, as can be seen from Figures 7c,d, the maximum ratio of local anomalies before the moderate-to-strong earthquake generally exceeded 8 times the normal value. Based on this, it is speculated that the August peak may not be a precursor to a moderate to strong earthquake, but is more likely related to minor seismic activity or meteorological factors. In October, the amplitude of anomalies in this region reached 14 times the normal level, causing a slight increase in the mean value and a significant deviation of twice the standard deviation. This fluctuation pattern indirectly suggests that the anomalies in October may be a potential precursor signal for a moderate to strong earthquake.

Due to differences in the spatiotemporal conditions of each earthquake, the spatiotemporal distribution of OLR anomalies preceding earthquakes exhibits diversity. Due to space constraints, this paper selects the 2023 Jishishan M6.2 earthquake (epicenter 102.79°E, 35.70°N, a recent strong earthquake of magnitude 6 or above) for a 2023 time-series curve analysis (Figure 5). It should be noted that this case study is only an individual example. The regional distribution characteristics of OLR anomalies preceding different earthquakes are similar, but also diverse and complex.

Moderate-magnitude seismicity (M5.0–6.0) predominantly clusters within the Qilian orogenic belt, Qaidam Basin, and Bayan Har block. Chronologically, the average annual frequency of this magnitude gear is 1.3 times in 2009–2019, while it enters the period of high-frequency activity in 2020–2024, and the average annual frequency rises to 2.8 times, with an increase of 115%. The OLR relative power spectral analysis shows (Table 2): For the M5.0-6.0 magnitude class, band 6 has the best retrospective detection performance with 82.1% (23/28) accuracy, followed by band 3 with 64.3% (18/28) accuracy. The statistics of spatial and temporal parameters show that 82.6% (19/23) of the earthquakes in the anomaly-generation time relationship occurred after the end of the anomaly, with the longest time from the beginning of the anomaly to the generation of earthquakes being 205 days, and the longest time from the end of the anomaly to the generation of earthquakes being 176 days, and that the anomaly window (beginning/end to the generation of earthquakes) for about 69.6% (16/23) of the seismic cases was centered within about 90 days. Spatially, the epicenters were distributed within a 200-km radius of the edge of the anomaly area, which is in line with the law of the “3-month window of seismicity”: within 3 months after the anomaly appeared, there was potential for earthquakes to be generated within the 200-km periphery of the anomaly area.

Larger events (M6.0–7.0) exhibit a bimodal spatial propagation pattern: east-west migration along the Qaidam Basin–eastern Tibetan Plateau tectonic corridor and north-south convergence toward the Qilian compressional front. In time, it gradually enters the “active period” from the “relatively quiet period”, with a low average of 0.4 earthquakes per year in 2009–2018, increasing to 0.6 earthquakes per year after 2019, as shown by the OLR relative power spectrum analysis (Table 2): For the M6.0-7.0 magnitude bracket, band 6 is able to achieve a 100% (7/7) detection rate, and band 7 is suboptimal with an accuracy of 85.7% (6/7), but with a wider spatial coverage of OLR power spectrum anomalies. In the time dimension, all earthquakes occurred after the end of the anomaly, with 85.7% (6/7) of the earthquake cases having an anomaly start-to-seismicity time ≤ about 120 days, and end-to-seismicity time ≤ about 100 days, with the only special case being the 2022 Delingha M6.0 earthquake, which had an anomaly start-to-seismicity interval of 216 days. Spatial analysis revealed all seismic events occurred within a 200-km radius of detected anomaly zones, indicating that the risk of M6.0–7.0 earthquakes is present in this range within 4 months after the appearance of the band 6 anomaly.

A total of two M ≥ 7.0 seismic events were recorded from 2009–2024: the 2010 Yushu M7.1 and the 2021 Maduo M7.4 earthquakes, both in the Bayan Har Massif. Band 3 OLR power spectral anomalies exhibited distinct precursory and post-seismic signatures during two major earthquakes. The Yushu M7.1 event displayed pre-seismic anomalies initiating 51 days before rupture, persisting for 21 days with peak intensity reaching 9 times the regional annual mean (3.0 × 10⁴ km² anomaly zone within 200 km of the epicenter). Conversely, Maduo M7.4 post-seismic anomalies emerged 11 days post-rupture, lasting 37 days with 16 times mean intensity and spatial expansion to 1.5 × 10⁵ km², where the epicenter resided within 200 km of the anomaly gradient boundary. Given the limited sample size (n = 2), these magnitude-specific patterns require validation through additional seismic events.

Based on the above laws, this study proposes the OLR power spectrum split-level prediction strategy, i.e., for M5.0–6.0 and M6.0–7.0 earthquakes: band 6 (with a characteristic period of 11 days) is prioritized for the precursor prediction of the OLR power spectra, and the seismic events mostly occur within 3 and 4 months after the anomalies, respectively, and need to pay attention to a 200-km radius around the periphery of anomalous zones; For M ≥ 7.0 earthquakes: anomalous signals in reference band 3 (with a characteristic period of 21 days), which need to be combined with an analysis of the regional tectonic stress field to improve reliability.

The relative power spectrum anomaly of the top-of-atmosphere outward longwave radiation, as a core parameter characterizing the radiative energy balance of the earth-air system, which can sensitively reflect the coupled surface-atmosphere perturbations caused by tectonic activities (Qiang et al., 1999). In order to deeply analyze the mechanism of OLR band 6’s pregnant seismic response to M5.0–6.0 and M6.0–7.0 earthquakes of two seismic grades in the northern region of the Tibetan Plateau, we performed regression modeling on the quantitative relationships between the OLR anomaly parameters (anomaly area A, the maximum ratio MR, and anomaly duration D) and the seismic energy parameters (seismicity generation energy E and seismicity generation grade M), respectively, and the results are shown in Figure 6.

For the M5.0–6.0 earthquakes, the anomalous area A (Figure 6a) and the maximu ratio MR (Figure 6c) of the OLR power spectra show an inflection point at about 17 × 10¹² J with the increase of the seismic energy E, and the overall trend is “increasing first and then decreasing”, and the correlation curves of both of them correspond to R² = 0.6734 and R² = 0.8785. In comparison, the latter fit better, its regression model was able to explain 87.85% of the variance and passed the significance test with p < 0.05 in the F-test, while the former fit was not as good as the latter, but its explanatory power also reached 67.34%, which also passed the significance test. The anomaly duration D (Figure 6e) increases with the increase of the seismic energy E, and there is an obvious decreasing trend around the seismic energy E = 15 × 10¹² J, showing the trend of “first flattening and then decreasing”, and the coefficient of determination is R² = 0.1518, which can only explain a small part of the changes, the fitting effect is poor, and in the F-test, p > 0.05, it does not pass the significance test. Therefore, the anomalous area A and the maximum ratio MR of the OLR power spectrum can be used as effective indicators for predicting the generation energy E of the M5.0-6.0 seismic event in the northern region of the Tibetan Plateau, and their empirical relationships are given by Equations 5, 6: A = − 1.38 × 10 − 2 E 3 + 0.42 E 2 − 2.26 E + 8.60 , R 2 = 0.6734 (5) M R = − 7.31 × 10 − 3 E 3 + 0.19 E 2 − 0.50 E + 11.09 , R 2 = 0.8785 (6) where A represents the anomalous area, MR represents the maximum ratio, and E represents the seismic energy, which can be obtained from I g E = 1.5 M w + 4.8 , where Mw is the seismic moment magnitude (Gutenberg and Richter, 1942).

For the M6.0–7.0 earthquakes, the anomalous area A (Figure 6b) and the maximum ratio MR (Figure 6d) of the OLR relative power spectra show an inflection point around the seismic energy E = 100 × 10¹² J and E = 360 × 10¹² J with the increase of the seismic energy E. The overall trend is “increase, then decrease, then increase”. The coefficient of determination of the correlation curve for the former is R² = 0.6593, which can explain 65.93% of the variance, while the coefficient of determination of the latter is R² = 0.1941, which is not a good fit, but neither of them passes the significance test of F. The anomalous duration D of the relative power spectrum of OLR (Figure 6f) shows an inflection point with the increasing of the energy of seismicity E, which occurs at around E = 320 × 10¹² J, and the overall trend is “first increasing and then decreasing”, and the coefficient of determination of the fitted curve is R² = 0.5043, which does not pass the significance test. Meanwhile, for the correlation fitting curve between anomaly area A and seismic energy E (Figure 6b), there is a long period of negative value of anomaly area A at seismic energy E = 250 × 10¹² J to 450 × 10¹² J, which is not physically meaningful. Therefore, there is no stable statistical correlation between OLR anomaly parameters and seismic energy for the M6.0-7.0 seismic events in this region, which may be limited by the sample size (n = 7) or the complexity of the generating mechanism, and it is not recommended to be used for energy prediction for the time being.

Therefore, it can be obtained that for the M5.0-6.0 magnitude class in the northern Tibetan Plateau, there is a good correlation between the anomalous parameters of the OLR power spectrum and the seismic energy E, and the seismic energy E can be predicted by the anomalous area A and the maximum ratio MR of the OLR power spectrum, which is of some significance for the prediction of the M5.0-6.0 magnitude class of the earthquakes in the northern Tibetan Plateau. Current analytical frameworks fail to establish statistically stable relationships between M6.0-7.0 seismic events and energy parameter E. This magnitude-specific correlation deficit persists despite rigorous empirical validation efforts. For this reason, this study further explores the statistical relationship between the OLR anomaly parameters and the corresponding seismic magnitude M, in order to reveal the differential response mechanism of the precursor signals of earthquakes of different magnitude classes.

For the M5.0–6.0 earthquakes, the anomalous area A of the OLR power spectrum (Figure 7a) shows a monotonically increasing trend with the increase of the seismicity M, and the anomalous duration D (Figure 7e) shows an overall trend of “increasing and then decreasing”, and the coefficients of determination of the related fitted curves are R² = 0.6281 (p < 0.05, passed the F test of significance), respectively, R² = 0.2646 (p > 0.05, did not pass F significance test). There is some statistical correlation between anomaly area A and earthquake magnitude M, but its goodness of fit to earthquake magnitude M is slightly lower than the goodness of fit of the regression model with seismic energy E, and its explanatory power decreases by about 5% compared with that of the seismic energy model (based on R² comparison). The goodness of fit of the anomaly duration D to the earthquake magnitude M is very low, indicating that it is not statistically significantly correlated with the earthquake magnitude M. The anomaly peak MR (Figure 7c) shows a three-stage nonlinear change, i.e., “decreasing, then increasing, then decreasing”, with a coefficient of determination of R² = 0.2602 (p > 0.05), which fails the F significance test, and the model can only explain 26.02% of the variance, which is insufficient to support its use as a reliable prediction indicator indicator. In conclusion, for the M5.0–6.0 earthquakes, the OLR power spectrum anomaly parameter has limited retrospective detection of the earthquake magnitude M. It is recommended to follow the strong correlation model with the seismic energy E for energy-constrained detection.

For the M6.0–7.0 earthquakes, the OLR relative power spectrum anomaly area A (Figure 7b) shows an overall monotonically increasing trend with increasing magnitude M, especially after magnitude M = 6.5. The coefficient of determination of this fitted curve is R² = 0.6576 (p > 0.05), which does not pass the F significance test. The maximum ratio MR (Figure 7d) and anomalous duration D (Figure 7f) show a three-phase trend of “first increase, then decrease, then increase” with the increase of seismic magnitude M, and both of them form a double inflection point at M = 6.2 and M = 6.7, and the coefficient of determination of the former is 0.9495 (p < 0.05), and the explanatory power of the polynomial model reaches nearly 95%, which is a good fit and passes the significance test, indicating that there is a strong nonlinear coupling between the maximum ratio MR and the magnitude M, while the latter has a coefficient of determination of R² = 0.6930 (p > 0.05), which does not pass the F significance test. Therefore, for the M6.0–7.0 earthquakes in the northern region of the Tibetan Plateau, the earthquake magnitude M can be predicted by the OLR relative power spectrum maximum ratio MR, which has an empirical relationship expressed by Equation 7: M R = 58.55 M 3 − 1.14 × 10 3 M 2 + 7.33 × 10 3 M − 1.57 × 10 4 , R 2 = 0.9495 (7) where MR is the maximum ratio and M is the earthquake magnitude.

Combining the above results, the preferred scheme of remote sensing precursor indicators for sub-seismic levels is proposed. For M5.0–6.0 earthquakes, the core indicators are anomaly area A (R² = 0.6734), the maximum ratio MR (R² = 0.8785), and the prediction target is the seismic energy E (the critical thresholds are both 17 × 10¹² J), which can be combined with the moment magnitude-energy conversion relationship to realize the energy-constrained magnitude interval estimation; For M6.0–7.0 earthquakes, the core indicator is the anomalous maximum ratio MR (R² = 0.9495), and the prediction target is the earthquake magnitude M. Segmented model and polynomial model joint constraints can be used to improve the magnitude prediction accuracy.

A systematic statistical analysis of 37 moderate-to-strong seismic events (2009–2024) was conducted across four principal tectonic units in the northern Tibetan Plateau: the Qilian Mountains orogenic belt, Bayan Har block, Qaidam Basin, and eastern Tibetan Plateau (Table 3). These geotectonic domains exhibit marked regional disparities in seismicity patterns and corresponding OLR anomaly responses, revealing unit-specific crust-atmosphere coupling mechanisms.

The most seismically active tectonic units in the region are the Qilian Mountains and the Bayan Har Mountains, both of which recorded 12 moderately strong seismic events (M ≥ 5.0), but there are two undetected seismic events in each region with an OLR power spectrum accuracy of 83.3%, and this good detection efficacy may stem from the fact that pre-seismic thermal radiation anomalies are more easily captured by remote sensing in arid zones. Among them, the Qilian Mountains area, as an important tectonic boundary on the northeastern margin of the Qinghai-Tibetan Plateau, is in the intersection area of the Qilian Mountains-Haiyuan Fault Zone and the Alashan Massif, and is subject to the tectonic coupling of the continuous eastward thrusting of the plateau and the blocking of the Ordos Massif, which has resulted in the formation of a belt of strong seismic activity dominated by the reverse type. Three M6.0–7.0 earthquakes were recorded during the study period, accounting for 42.9% of the total number of earthquakes of this magnitude file. The Bayan Har region shows anomalous OLR power spectra in band 3 of the OLR power spectrum for the two M ≥ 7.0 earthquakes, which we hypothesize may be related to the long-period radiative changes caused by the migration of deep fluids prior to the rupture of the large rupture zones (Ma et al., 2024).

The Qaidam Basin shows unique seismic detection characteristics, with 100% detection accuracy of 10 M5.0-7.0 seismic events (including the 2009 Haixi M6.4 earthquake) in the 6-band of the OLR power spectrum. This high seismicity and high OLR power spectrum detection accuracy may be related to its “basin effect”, i.e., the anomalies associated with the basin may be due to the oil and gas resources in the basin, and when the stress accumulates to a certain extent, the active tectonic zones and some microfractures at the periphery of the basin become the upward influx of greenhouse gases, which leads to a significant enhancement of the radiation effect. Enhancement effect is obvious (Zhang et al., 2018). Seismicity in the eastern Tibetan Plateau region exhibits obvious grading characteristics. The OLR power spectrum achieves complete detection for two M6.0-7.0 magnitude grade seismic events, but there is a missed detection for the M5.0–6.0 earthquake (Xiahe M5.7, 2019). This magnitude detection threshold may originate from the double effect of topography and cloud cover.

Comprehensive analysis reveals distinct seismogenic patterns across northern Tibetan Plateau’s tectonic units. The Qilian Mountains exhibit high-frequency moderate-to-strong seismicity (M5.0–7.0), where Band 6 OLR spectral ratios demonstrated 83.3% retrospective detection accuracy. In contrast, the Bayan Har block requires dual-band monitoring—Band 6 anomalies for general seismicity and Band 3 signatures for M ≥ 7.0 events. The Qaidam Basin showed predominant M5.0–6.0 activity (single M6.0–7.0 event since 2009) with 100% Band 6 detection efficacy, likely enhanced by basin amplification effects. Conversely, the eastern Tibetan Plateau presents detection challenges for M5.0–6.0 events via OLR spectral analysis, necessitating prioritized monitoring of M6.0–7.0 seismicity through integrated geophysical approaches.

Analyzed from a fault zone perspective (Table 4), the Sunan-Qilian fault zone and the northern Qaidam fault zone show the strongest seismicity. Among them, 11 seismic events occurred in the Sunan-Qilian fault zone during the observation period, including three M6.0–7.0 earthquakes, and their high-frequency activities may be related to the composite kinematic characteristics of this rift zone, which combines both thrust tectonics and strike-slip components. It is noteworthy that the OLR power spectrum anomaly showed spatial migration characteristics before the 2022 Menyuan M6.9 earthquake and formed a thermal radiation anomaly zone about 150 km long along the fault zone. Five seismic events, including the 2009 Haixi M6.4 earthquake, were recorded in the northern Qaidam fault zone. The eastern section of this fault zone (near Delingha) exhibits more intensive medium- and strong-seismic activity, which may be related to the non-uniform distribution of tectonic stresses by the rigid base of the Qaidam Basin (Peng et al., 2012).

The Yushu-Ganzi Fault Zone and the Maduo-Gande Fault Zone show similar differences in detection effectiveness. Both of them detected one M5.0–6.0 earthquake and one M ≥ 7.0 earthquake during the study period. On the other hand, the Northwestern Qinling Fault Zone showed different detection efficacy, i.e., it accurately tested for two M6.0–7.0 earthquakes, but the test results were unsatisfactory for the M5.0–6.0 earthquakes. We consider that this difference may reflect the influence of different fracture zones’ sliding rates on the thermal radiation precursor signals - the sliding rates of the Yushu-Ganzi fault zone and the Maduo-Gande fault zone are 5–6 mm/yr and 4–6 mm/yr, in that order, which are much lower than the sliding rate of 12 mm/yr for the fracture zone at the northern margin of the West Qinling Mountains (Qiao et al., 2022; Huang et al., 2011), and the fracture zones with high slip rates are more prone to long-term thermodynamic perturbations due to continuous friction and creep-slip activities. Therefore, the prediction of medium-sized earthquakes (Xiahe M5.7 in 2019) is poor, whereas low-rate faults are usually accompanied by high-stress confinement with energy accumulating in the form of elastic strain, and thermal anomalies in such regions may be concentrated in the pre-slip or microfracture extension phase in the short term (weeks or months, e.g., anomalies of the Yushu M7.1 earthquake in 2010 started 51 days prior to the earthquake) prior to strong earthquakes (Song et al., 1944).

The year was divided into 4 seasons based on natural months: spring from March-May, summer from June-August, fall from September-November, and winter from December-February of the following year. Classify the season in which the OLR power spectrum anomaly is located based on when it begins and ends (if the anomaly duration spans 2 seasons, it is classified according to the season in which it is located for the greater amount of time) (Zhang et al., 2024).

Table 5 presents seasonal seismicity distributions in the study area during 2009–2024. Table 6 analyzes seasonal occurrence frequencies of OLR power spectral anomalies associated with individual seismic events. According to the statistical analysis in Table 5, analysis reveals quasi-uniform seasonal seismicity distribution: (1) spring (10 events), summer (9), autumn (8), and winter (9); (2)Moderate-to-strong earthquakes (M6.0–7.0) demonstrate winter prevalence (42.9% frequency), while major events (M ≥ 7.0) exclusively cluster in spring; (3) The OLR relative power spectrum anomaly detection has full capture capability for springtime seismic events, but partial misses for other seasons, with the most significant misses in the fall, where the miss rate reaches 25% of the total.

Seasonal analysis of OLR power spectrum anomalies (Table 6) reveals pronounced spatiotemporal heterogeneity in their distribution patterns: (1) in the seasonal dimension, summer anomalies dominated the observed events, constituting 34.4% of total occurrences, while winter displayed the lowest incidence (6.3%), which may stem from the high temperature and relative atmospheric stability in summer, which is more favorable for thermal infrared anomaly identification; (2) In the magnitude dimension, M5.0-6.0 seismic anomalies occur mostly in spring (34.8%), M6.0-7.0 seismic anomalies are concentrated in summer and fall (42.9% each), and M ≥ 7.0 seismic anomalies are seen only in spring and summer (50% each).

Integrated with tectonic subdivision analysis (Table 7), OLR anomalies demonstrate distinct seasonal preferences across regions: (1) Qaidam Basin exhibited year-round anomalies with summer dominance (40% frequency), followed by spring, winter, and autumn; (2) Qilian Mountains showed concentrated summer occurrences (50%), complemented by autumn events (30%), with no winter anomalies detected; (3) Bayan Har Mountains displayed primary spring anomalies (50%), followed by autumn and summer, potentially associated with spring temperature increases and associated surface thermal radiation changes during snowmelt periods; (4) Eastern Qinghai-Tibet Plateau exclusively presented two autumn anomalies, likely attributable to enhanced thermal infrared detection efficiency under reduced cloud cover conditions during this season (Liu et al., 2025).

Taken together, the OLR power spectrum anomalies occur most frequently in summer and least frequently in winter. This may be attributed to: (1) the generally higher ambient temperature in summer, which enhances the thermal radiation of the object; (2) the better atmospheric transparency in summer and the relative stability of the atmosphere, which is less prone to complex meteorological conditions such as inversions at low temperatures, which reduces the scattering and absorption of the signals, and is conducive to the propagation of the thermal infrared signals; (3) Vegetation is capable of emitting thermal infrared radiation, which can be enhanced to some extent by the luxuriant growth of vegetation in summer. While the second highest frequency of anomalies in fall may be attributed to the reduced cloud cover in fall and the lithosphere-atmosphere coupling effect induced by the seasonal adjustment of atmospheric circulation, the scarcity of anomalies in winter may be related to the attenuation or even shielding effect of snow cover on the OLR radiation as well as the attenuation of the surface radiation by the inversion layer (Wang et al., 2018). However, there are some limitations in this part of the study. On the one hand, the small sample of winter earthquake cases (n = 2) affects the statistical significance of this study, and on the other hand, the modulation effect of topographic elevation on the thermal radiation transmission is not taken into account, which needs to be further validated by combining the InSAR surface deformation data with the observation of geothermal temperature field in a multi-parameter verification.

Based on the statistics of the source mechanism solutions for 37 earthquakes in the northern Tibetan Plateau region from 2009 to 2024 (Table 8), this study reveals the correspondence between different fault types and OLR power spectrum anomalies. The statistical results show that (1) strike-silp earthquakes account for the highest proportion (45.9%, 17/37), and their OLR anomaly detection accuracy is 82.4% (14/17), and the missed cases are concentrated in the M5.0-6.0 magnitude, with a total of three undetected cases; (2) the reverse-fault type earthquakes were the next most common (27%, 10/37), with a detection accuracy of 80% (8/10); (3) the detection accuracy of the anomalies of the normal-fault type (n = 7) and thrust-fault type (n = 3) earthquakes was 100%, with the normal-fault type containing a single M6.0–7.0 earthquake and the thrust-fault type being distributed across the magnitude brackets. In terms of detection effectiveness, the OLR power spectra predicted 43.8%, 9.4%, 21.9%, and 25% of the strike-slip-type, thrust-type, normal-faulting, and reverse-faulting earthquakes, respectively.

The study area is dominated by strike-slip type earthquakes (46%), which is closely related to the tectonic context of the northern Tibetan Plateau subjected to eastward extrusion shortening and lateral slip of the block (Zheng et al., 2016). Moreover, its high anomaly detection rate (82.4%) may be significantly related to the frictional heat generation effect during rapid sliding of the strike-slip faults, and such heat radiation signals are easily captured by satellite remote sensing. Although there are fewer samples of normal and thrust-strike earthquakes, the accuracy of OLR detection reaches 100%, which may reflect that the shallow thermal releases associated with tensional/compressive ruptures are more likely to penetrate the surface cover. The missed detection rate of 20% for the reverse earthquakes is presumably related to the dispersed spatial distribution of thermal anomalies due to the slower dip of the reverse faults, and such signals may be masked by cloud cover or topographic noise (Qin et al., 2010).

The seismic magnitude prediction model proposed in this study reveals differences in the response of OLR anomalies in different frequency bands to earthquakes of different magnitudes. The high detection rate (85.7%, 30/35) of M5.0–7.0 earthquakes in frequency band 6 (characteristic period of 11 days) may be related to its periodic characteristics and the energy release pattern of the earthquake incubation process. According to rock fracture experiments, the prenecessary stage of moderate to strong earthquakes (M5.0–7.0) is typically accompanied by accelerated microfracturing and expansion of the fracture network. Such processes may trigger short-period (7–15 days) thermal radiation fluctuations through surface-atmosphere coupling effects (Pulinets and Ouzounov, 2011). Abnormal signals from earthquakes with magnitudes greater than 7.0 tend to occur more frequently in band 3 (with a characteristic period of 21 days), which may reflect long-term disturbances caused by pre-slip or deep fluid migration in the locked segment of large fault zones during the late stages of stress accumulation. That is, due to changes in strain within the pores of the Earth’s lithosphere, large fault zones may experience pre-slip or restricted fluid stretching, leading to interstitial fluid flow and resulting in the redistribution of deep fluids (Fleischer, 1981). Dobrovolsky et al. (1979) provided an alternative description of the long-range stresses associated with earthquakes, which reasonably explains this phenomenon. For example, the anomaly that occurred 51 days before the M7.1 earthquake in Yushu was related to plastic deformation in the low-velocity layer of the lower crust of the Qinghai-Tibet Plateau. The thermal anomaly may have been caused by the combined effects of friction heat from the fault zone and fluid upwelling.

It is worth noting the magnitude-dependent features of the anomaly parameters. The anomalous area A of the OLR power spectrum and the anomalous maximum ratio MR of the M5.0–6.0 earthquakes show a nonlinear correlation with the energy E (R2 of 0.6734 and 0.8785, in that order). Among these, energy E is calculated based on the Gutenberg-Richter formula using the moment magnitude Mw (Gutenberg and Richter, 1942). This model is similar to the dislocation model proposed by Fleischer (1981), which derives the relationship between strain distance and magnitude using an empirical strain-energy equation to aid in predicting the time of earthquake occurrence (Fleischer, 1981). However, the difference lies in the fact that the dislocation model calculates energy E using the surface wave magnitude Ms according to Stacey (1969) formula E ergs = 10 12.24 + 1.44 M . In contrast, this paper uses the moment magnitude Mw, which can accurately reflect the energy released by an earthquake, thereby improving the reliability and accuracy of the magnitude prediction model. Whereas the anomalous maximum ratio MR of the M6.0–7.0 earthquakes exhibits a curvilinear response to the seismic magnitude M (R2 of 0.9495). This feature is similar to the empirical model of precursor anomalies and seismic magnitude proposed by Ramola et al. (1988) based on the radon transport mechanism. By analyzing 17 major seismic data sets worldwide, the model found that in the empirical calculation results of seismic magnitude based on changes in the shape of radon anomalies, eight seismic events deviated from the actual seismic magnitude by more than 0.3 magnitudes, of which only two were M6.0–7.0 earthquakes (accounting for 25%). This suggests that, for the M6.0-7.0 magnitude range, empirical relationships between magnitude and seismic precursors based on abnormal seismic precursors may be a potential avenue for improving the effectiveness of earthquake prediction. The differences between these prediction models may reflect differences in the scale of earthquake precursors for earthquakes of different magnitudes: thermal anomalies associated with small and medium-sized earthquakes are mainly caused by local stress adjustments, while precursory signals for large earthquakes may involve the coordinated action of tectonic systems over a larger area.

Virk and Singh (1994) studied the 1991 M7.0 earthquake in Uttarkashi, India, and found that the abnormal radon concentration prior to the earthquake was caused by local fluid migration triggered by distant stress. This study elucidated the mechanism of radon as a precursor to earthquakes from a structural physics perspective. In the 1966 M5.4 earthquake in the Garm region of the Soviet Union, the P/S wave velocity ratio (Vp/Vs) decreased significantly before the earthquake, and the duration of the anomaly was related to the magnitude of the earthquake. This phenomenon was attributed to changes in the microfractures of the crustal rocks and fluid permeability caused by stress (Rikitake, 1976). The above study reveals that the mechanism of earthquake precursors is related to changes in medium properties (such as microfractures and fluid migration) caused by crustal stress. Therefore, when using OLR power spectra to predict earthquakes, in addition to the model based on magnitude and energy described in this paper, analyzing the coupling mechanism between OLR anomalies and crustal stress-strain will be a key focus of future research.

The spatial variability of OLR power spectrum anomalies in the northern Tibetan Plateau region is significantly controlled by the active tectonic pattern. The high detection rates (100%, 83.3%, and 83.3%, respectively) in the Qaidam Basin, Qilian Mountains, and Bayan Har Mountains can be attributed to the following mechanisms: (1) the stress-shielding effect of the rigid basement of Qaidam on the peripheral ruptures, which concentrates the energy release at the basin margins and creates a significant gradient of thermal radiation (Zhang et al., 2018); (2) The arid areas of the Qilian Mountains and the Bayan Har Mountains have low vegetation cover, and the surface thermal radiation signals are highly penetrating and more easily captured by satellite remote sensing. On the contrary, the missed detection phenomenon in the eastern Tibetan Plateau may be related to the higher surface thermal inertia of alpine meadows, and such environments will have a certain damping effect on thermal anomalies (Liu et al., 1999).

For the fracture zones, the modulation effect of their sliding rates on the OLR anomalies should not be neglected. High-rate slip fractures (e.g., the Northwestern Qinling Fault, with a sliding rate of 12 mm/yr) are more prone to long-term dynamical perturbations due to the release of some strain energy from continuous creep-slip, while low-rate fracture zones (e.g., the Ganzi-Yushu Fault Zone, with a sliding rate of 5–6 mm/yr) are dominated by short-term strong anomalies. In addition, the 1.5 × 10⁵ km² anomaly area observed after the M7.4 earthquake in Maduo in 2021 may be the result of large-scale CO₂ emissions within the tensile section of the East Kunlun Fault Zone. This phenomenon implies that the combination of gas geochemical observations can enhance the interpretation reliability of OLR anomalies.

The seasonal distribution of OLR anomalies (34.4% in summer) is characterized by a clear connotation of the dynamics of the earth-air system. Increased solar radiation in summer, more energy received at the surface, and significant surface warming enhance the release of longwave radiation from the surface to the atmosphere (Wei et al., 2024). At the same time, the atmosphere is relatively stable in summer, with fewer complex meteorological conditions such as inversions at low temperatures, which reduces the scattering and absorption of signals and facilitates the propagation of thermal infrared signals (Chen, 2015). In addition, summer is full of sunlight and the atmosphere is more transparent, which is conducive to the acquisition of clear remote sensing images by thermal infrared sensors.

The cause of the scarcity of winter anomalies (6.3%) may involve multi-scale coupled processes: (1) snow cover (thickness > 10 cm) can reflect up to 80% of surface longwave radiation, significantly attenuating the signal of true thermal anomalies (Wang et al., 2018); and (2) frequent occurrence of an inversion layer leads to attenuation of upward transmission of surface radiation (Qin et al., 2010); (3) The crustal loading rate slows down in winter (GPS observations show that the strain rate decreases by about 30%), and the energy accumulation of the seismogenic system enters into a relatively silent period (Wang et al., 2018). However, the winter anomaly case of the 2022 Menyuan M6.9 earthquake suggests that large-scale strike-slip earthquakes may break through the seasonal limitation, and its thermal anomaly generation mechanism needs to be further investigated.

Differences in the thermal radiation response of different earthquake source types reveal the role of fault kinematics in modulating precursor signals. The high detection rate (82.4%) of strike-slip type earthquakes may be due to two aspects: (1) the steeper dip angle of strike-slip faults (usually > 70°), which makes it easier for heat flux generated by frictional heat generation to be vertically conducted to the surface; (2) the shear heating effect accompanying strike-slip motion can increase the temperature of the fault zone, significantly changing the local thermal radiation background field (Wu et al., 2004). In contrast, the high miss detection rate (20%) of reverse-type earthquakes may be due to the horizontal diffusion of thermal anomalies due to low-angle retroflex (dip < 30°), and satellite data with insufficient spatial resolution (2 km spatial resolution for Feng Yun geostationary meteorological satellites) have difficulty in recognizing such diffuse signals. Notably, the anomalous signal intensity (maximum ratio exceeding 11 times the annual average) is significantly higher for the normal fault type of earthquake than for the other types, which may reflect that the rapid opening of the fissure accompanying tensional rupture promotes efficient upwelling of fluids from the depths.

Regarding seismic thermal radiation anomalies, experts and scholars have proposed various hypotheses regarding their causes, such as the gas heat theory and the stress-induced heat theory. Qiang et al. (1999) “Earth gas release theory” posits that fractures in fault zones allow gases accumulated underground to diffuse to the surface, where they emit invisible infrared light under the influence of low-altitude electric fields, absorb solar and ground radiation, and trigger a greenhouse effect that causes regional warming. Zhang et al. (2010) proposed the “Earth greenhouse effect-atmospheric coupling” mechanism for the cause of seismic thermal infrared anomalies, pointing out that these anomalies are caused by the combined effects of basic and indirect causes. The basic cause (Earth greenhouse effect) is related to rapid accumulation of regional stress, active tectonics, and geological environment; the indirect cause (atmospheric coupling) is related to regional atmospheric water vapor content. Cai et al. (1987) used a fault creep model to infer that the abnormal ground temperature during the Tangshan earthquake may have been caused by accelerated fault creep prior to the earthquake. Xu et al. (1991) verified the “electromagnetic warming effect theory” and found that rock conductivity is highly sensitive to temperature. Increased crustal activity can cause the release of large amounts of greenhouse gases, forming an electromagnetic anomaly field. Under the influence of an electric field, greenhouse gases can produce a warming effect. Wu et al. (2004) experimentally confirmed that rock fractures are preceded by thermal infrared precursors, including two forms: abnormal average infrared radiation temperature-time curves and thermal infrared radiation anomalies. Furthermore, thermal infrared radiation anomalies exhibit spatial migration characteristics prior to fracture, developing in the sequence of “appearance-intensification-calm-fracture.”

In summary, currently, regarding the physical mechanisms underlying seismic thermal radiation anomalies, scholars have proposed a wide range of hypotheses and explanations from multiple perspectives, including Earth outgassing, greenhouse effect-atmospheric coupling, fault creep, electromagnetic heating, and rock fracturing. These research findings provide diverse perspectives for our deeper understanding of the phenomenon of seismic thermal radiation anomalies. However, due to the complexity and uncertainty of the seismic system, a unified and definitive theoretical framework has not yet been established, and further exploration and verification are still needed in the future.

Previous studies have thoroughly confirmed that pre-earthquake thermal radiation anomalies are evident in a large number of observations. From a spatial perspective, thermal radiation anomalies generally exhibit a dynamic process of “onset-intensification-peak-attenuation-disappearance.” In terms of anomaly intensity, when reaching its peak, the area exceeding twice the standard deviation of the long-term average can span tens of thousands to hundreds of thousands of square kilometers, with the maximum relative power spectrum amplitude typically exceeding six times the average. From a temporal perspective, such anomalies typically emerge within 3 months prior to an earthquake, with durations ranging from 15 to 120 days (Guo et al., 2010; Wei et al., 2013; Zhang et al., 2024).

However, most studies have focused on only one or a few aspects of pre-earthquake precursors, and during the analysis process, they often use only a single parameter of OLR anomalies (such as the area of anomaly) and simply compare it with information from a single earthquake case. While this research approach can reveal some characteristics of pre-earthquake thermal radiation anomalies to a certain extent, it lacks comprehensiveness and systematicness. As a result, it is difficult to form a complete understanding of pre-earthquake thermal radiation anomaly phenomena, and there are certain limitations.

Therefore, this paper focuses its research perspective on the entire northern region of the Qinghai-Tibet Plateau. By conducting a comprehensive review and summary of all OLR anomalies preceding major earthquakes in this region over the past 15 years, the study not only statistically analyzed multiple key parameters of OLR anomalies but also introduced the Gutenberg-Richter relationship formula to precisely calculate the energy released prior to earthquakes using the moment magnitude (Mw) scale. The calculated results were then deeply integrated with OLR anomaly parameters. Building on this foundation, a magnitude-classified prediction model specifically tailored to the northern Tibetan Plateau has been developed. This research represents a significant breakthrough in terms of scope, parameter application, and the construction of prediction paradigms, providing new insights and methods for understanding the seismic precursors of the northern Tibetan Plateau and enhancing earthquake prediction capabilities.

Kang et al. (2011) used OLR grid data from NOAA polar-orbiting meteorological satellites to extract seismic anomaly information through vorticity quantification calculations. Their anomaly determination criteria relied solely on a significant increase in vorticity values 1–2 days before the earthquake. This method can be used for strong earthquake monitoring in many regions around the world; Fu et al. (2020) used daily OLR data obtained from NOAA satellites and calculated the difference between daily vorticity and the average vorticity of the previous 16 days to characterize OLR anomalies. This method took into account the influence of the hypocenter depth and focused on monitoring abnormal signals that appeared 2–15 days before the earthquake. Unlike previous studies, this paper utilizes OLR data from the FY series of geostationary meteorological satellites, dividing the signal into seven characteristic frequency bands. Based on anomalies in different frequency bands, the magnitude of the earthquake can be preliminarily determined. Further analysis combines the constructed fitting model and the coefficient of determination R² to determine the earthquake magnitude and energy. The coefficient of determination R² can reach up to 0.9495, indicating high fitting accuracy. Additionally, this method integrates various influencing factors such as regional tectonic units, seasonal changes, and seismic source mechanism solutions, enabling effective monitoring of OLR anomalies 1–4 months prior to moderate to strong earthquakes. Compared to the aforementioned two methods, more influencing factors are considered, resulting in higher fitting accuracy and providing more reliable precursory evidence.

Since 2007, the FY satellite has entered a stable operational phase and achieved routine OLR data collection. However, the limited number of strong earthquakes (M ≥ 7.0) in the study area (n = 2) casts doubt on the statistical significance of Band 3, and the generalizability of the constructed model lacks sufficient support. More data is needed for validation, and future research should utilize global datasets for extensive validation to enhance the model’s accuracy and practicality. Additionally, the terrain of the Qinghai-Tibet Plateau is complex, with the Bayan Kara Mountains exceeding 5,000 m in elevation, and the attenuation effect of OLR radiation transmission has not yet been quantified. Future studies should incorporate terrain correction models such as the Minard model to accurately assess terrain impacts. The current study relies on a single data source, using only OLR anomaly parameters. Future research should integrate multi-source remote sensing data such as surface temperature and aerosol optical depth to construct multi-parameter coupled prediction indices, thereby enhancing the scientific rigor of predictions.