Conventional methods for extracting seismomagnetic anomalies from ground records often 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 at （33.20°N, 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 Geomagnetic Network Center of China, Institute of Geophysics, China Earthquake Administration （within 700 km of the epicenter of Jiuzhaigou M_S7.0 earthquake） with a theoretical external field. The theoretical external field was constructed from two Swarm Level-2 products: 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 were taken 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 converted to station-wise standard scores using the quiet-day mean and standard deviation.

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 quiet-day mean and standard deviation of each harmonic parameter were used to define a 2σ empirical threshold. On 9 August 2017, one day after the Jiuzhaigou earthquake, 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 numbers of anomalies were also detected on several pre-seismic quiet days, particularly on 9 June and 21 June.

The 9 August anomalies showed a clear north-south asymmetry across the epicentral area. In the northern to northwestern part of the epicentral area, the first-order phase difference was positive, indicating a phase lag, whereas in the southern part it approached zero. The second-order amplitude ratio was generally below unity at northern and northwestern stations, indicating suppression of the observed Z-component amplitude relative to the model background; in contrast, some southern stations showed relatively higher |G| values, suggesting local enhancement relative to the northern low-value zone. The spatial fields were mapped by radial basis function interpolation and summarized using sector-based statistics, both confirming the asymmetry.

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.