Methods for extracting seismomagnetic anomalies from ground records depend on empirically selected reference stations and the local electromagnetic environment. This work builds a reference-station-independent diurnal background from satellite external-field models, and uses it to extract coseismic geomagnetic changes of the 8 August 2017 Jiuzhaigou earthquake. The Jiuzhaigou earthquake occurred on 8 August 2017 at （33.20°N and 103.82°E）, with a magnitude of M_S7.0 and a focal depth of 20 km. The rupture had a strike of 150°, a dip of 80°, and a rake of −20°, and the source region has complex geological structures with numerous active faults such as Tazang fault, Minjiang fault, and Huya fault, and is historically prone to strong earthquakes.

We compared hourly means of the vertical component （Z） from 39 stations of the China Geomagnetic Network （within 700 km of the epicenter of Jiuzhaigou M_S7.0 earthquake） with a theoretical external field. This external field consists of two level-2 products from Swarm satellite: one is the ionospheric model MIO, derived by specialized ionospheric field inversion, which represents the solar-quiet current system at middle and low latitudes, with its variation modulated by the solar radio flux index F_10.7; the other is the magnetospheric model MMA, derived by comprehensive inversion, which represents large-scale magnetospheric current systems. Because the study area lies far from the equator, the contribution of the equatorial electrojet to the ionospheric field was neglected. The static internal field, comprising the main field up to degree 20 and the lithospheric field up to degree 185, was removed using the CHAOS-8.2 model. This process yielded a variation field dominated by external primary and induced components. The combined model results were taken as the normal diurnal background, while the polynomial-corrected observations as the local geomagnetic response.

The 39 stations comprise eight reference stations and 31 basic stations, and both baseline-corrected absolute and relative Z observations were used. For absolute data the CHAOS-8.2 main and lithospheric fields were subtracted; for relative data the recorded values were used directly. Data gaps were filled by distance-weighted K-nearest-neighbor interpolation with seven neighbors. A one-day sliding window with a fifth-order polynomial aligned the observations to the model and removed direct-current offsets left by incomplete static-field subtraction, retaining the high-frequency diurnal form of the observations. Harmonic coefficients of the observed and modeled fields were estimated by least-squares harmonic estimation （LSHE） for the first three orders, that is the 24 h, 12 h, and 8 h waves. From each set of coefficients, the amplitude and phase of every order were obtained. With the model as input and the observation as output, a complex transfer function per order yielded two anomaly indicators: the amplitude ratio |G| （the ratio of observed amplitude to modeled amplitude） and the phase difference Δφ （observed phase minus modeled phase）. Daily values of |G| and Δφ were standardized by Z-score across the 39 stations to remove differences in scale between stations and orders.

For March to October 2017, days with disturbance storm-time index |D_st|＞20 nT or planetary index K_p＞2 were removed to form a quiet-day baseline. On non-seismic quiet days the corrected observations matched the model field in form and amplitude, so the satellite-model background reproduces the diurnal field at middle and low latitudes at hourly resolution without a reference station. For each station the daily mean and standard deviation （σ） of |G| and Δφ defined a 2σ threshold. On 9 August 2017, one day after the main shock of Jiuzhaigou earthquake swarm, the first-order phase difference exceeded the threshold at 32 of 39 stations; the second-order |G| and Δφ exceeded it at 22 and 20 stations, respectively; and the third-order |G| exceeded it at 14 stations. Phase and amplitude were distorted across the first three orders on the same day, coinciding in time with the earthquake. Elevated anomaly counts also occurred on several pre-seismic quiet days, especially on 9 June and 21 June.

The 9 August anomalies were asymmetric across the epicenter. To the north of the epicenter （azimuth 255° to 30°） the first-order phase difference was positive （phase lag）, while to the south of the epicenter it approached zero. The second-order amplitude ratio was below unity （suppression） at northern and northwestern stations and above unity （enhancement） at some southern stations. The spatial fields were mapped by radial basis function interpolation and summarized by azimuth, both confirming the asymmetry. The band of |G|＜1 trended from NW-SE to nearly E-W, which coincided with the southern side of the Qilian-Haiyuan-Liupanshan-Qinling fault zone along the northeastern margin of the Tibetan Plateau, where the crust is thick and electrically conductive. Published thermal infrared and total electron content anomalies for this earthquake concentrate in the region south of the epicenter, which is linked to fluid and gas release at the basin margin.

The combined pattern is consistent with a two-zone response controlled by the regional structure: to the south of the epicenter, at the Sichuan basin margin, shallow micro-fracturing and fluid upwelling raise near-surface conductivity and strengthen the induced, thermal, and ionospheric signals; to the north of the epicenter, the thick high-conductivity crust of the plateau margin shields and suppresses the induced field, so the observed vertical component falls below the modeled background. The harmonic anomalies are frequency-selective, concentrating in the first-order phase difference and the second-order amplitude ratio, which points to induced currents at different depths. Their spatial distribution follows the regional fault system rather than radial distance from the source, which links the coseismic geomagnetic change to the deep electrical structure as well as to local stress release.