为了通过1200–2000重建普华永道的可变性，我们采用了多种方法的多种方法。现代全球降水Δ18O与PWC高度相关，这是各种良好机制和远程连接的结果22。我们利用了使用全球分布的水 - 同位素代理记录网络的关系，来自五种不同的代理档案类型：冰川冰，木材，湖泊沉积物，珊瑚和Speleothem。使用不同的存档类型降低了档案特异性偏见的风险，例如，偏见“温暖”或“湿”季节值，同时还允许包括最多可能的记录。这很重要，因为与单个站点相比，网站网络对非平稳远程连接更强大38,54。我们使用了八种统计方法（“重建方法”部分）来隔离PWC信号，从而考虑了特定方法的偏见。我们还使用了几个目标数据集来说明观察不确定性（“观察数据源”部分）的影响，并包括对年代不确定性的强大治疗方法（“时间顺序不确定性的结合”部分）。
重建靶标的是反太平洋赤道ΔSLP，定义为中部太平洋中部均值SLP之间的差异（160°–80°W，5°S – 5°N）和西太平洋/西太平洋/东部印度大洋洲（80°–160°–160°–160°– 160°– 5°，5°S-5°S-5°n – 1 – 5°）。ΔSLP与PWC8,9,11,22,47,55的强度密切相关，并且与更复杂的基于PWC的强度高度相关，PWC的强度仅适用于1950年，可用于现场56。
我们使用两种网格的观察产物（HADSLP和ICOADS）和一个大气再分析产物（ERA-20C）计算了ΔSLP。HADSLP在跨越1850年的5°分辨率下可用，源自质量控制的海洋和陆地SLP观测结果25。在1900 - 2004年，我们使用了“ HADSLP2”产品；从2005年开始，我们使用了“ HADSLP2R”产品。来自ICOAD的SLP数据可在2°分辨率下可用，跨越1800到现在，并源自Surance Marine观察数据26。ERA-20C将表面压力和海洋风异常吸收到大气中的通用循环模型27中，并以大约1°的分辨率获得，跨越1900年至2010年。
大多数代理数据来自ISO2K数据库（n = 50），这是Water23稳定同位素组成的代理记录的多座汇编23。经过广泛的文献搜索，我们采购了九个记录，其中作者描述了代理记录与普华永道或ENSO之间的牢固关系。该类别中的speleothemΔ18O数据来自Sisal数据库57。对于ISO2K记录，我们仅保留每个record23的指定的“主要”时间序列，然后仅考虑每年或次生解析的代理记录，数据至少扩展到2000。（参考文献58,59,60,61,62,63,64,65,66,66,67,68,69,70,71,71,72,73,73,74,75,76,77,77,78,79,79,80,81,81,82，83,84,85,86,87,89,90,91,92,92,93,94,95,96,97,97,98,99,99,100,101,102,102,103,104,105,106）。
我们使用了CESM1 LME（参考文献40）以及五个PMIP3（参考文献42）和三个PMIP4（参考文献43）贡献的数据，以与我们的PWC重建进行比较。我们使用了以下CMIP5/PMIP3模型的千年模拟：BCC-CSM1.1（参考文献107），CCSM4（参考文献108），Fgoals-S2（参考文献109）（参考文献109），GISS-E2-R（参考文献110）和MRI-CGCM3（参考111）。我们使用了以下CMIP6/PMIP4模型的千年模拟：INM-CM48（参考文献112），MiroC-ES2L（参考文献113）和MRI-ESM2.0（参考文献114）。
我们仅使用CESM1 LME成员使用所有人为和天然的外部强迫因子，即完全强制的集合成员（n = 13）。PMIP3数据包括一个额外的单模集合（GISS-E2-R，n = 8）。
我们选择了重建的1200-2000间隔，因为这在代理数据的可用性和长期强迫和内部变异性的采样之间取得了最佳平衡。
我们使用了两个校准窗口。对于主要文本中提出的重建，我们使用了1900-2000，以最大程度地减少非平稳远程连接的影响18。1900年是ERA-20C最早覆盖的年度，2000年的结束年份提供了最大化校准窗口长度和随附的代理记录的数量的最佳平衡。
我们使用较短的校准窗口（1951-2000）重新计算了完整的ΔSLP重建集合，从1900 - 1950年间执行的独立验证测试（“评估重建技能”部分）提供了重建技能的最低估计值。
我们重建了年度ΔSLP，允许表征长期的年际PWC变异性和低频变异性。大多数重建方法都需要在共同的时间步骤上进行数据，因此每年将次年解决的记录归为日历年（1月至12月）。装订后，我们将记录保留在校准窗口（1900-2000）中的三分之二的垃圾箱中。我们使用数据插值经验正交函数（DINEOF）方法估算了校准窗口中的任何丢失年度，该方法以维持潜在共同点的方式插值丢失值115。
我们的三种重建方法要求贡献记录与目标索引相关。在这种情况下，我们仅保留了大量记录（P< 0.1) correlated with ΔSLP over the calibration window.
Because the number of available records that extend to 2000 decreases with increasing record length (Extended Data Figs. 1 and 2), we performed all reconstructions over five temporal subsets: 1860–2000, 1800–2000, 1600–2000, 1400–2000, and 1200–2000, following the ‘nested’ approach of previous studies24,116,117.
For each subset, we only included records with data in greater than or equal to two-thirds of the years spanning the entire interval. All methods except one require continuous data, so we interpolated missing data using the DINEOF method (Extended Data Fig. 1). To avoid spurious jumps when appending segments, we aligned each older segment (for example, 1400–2000) with the adjacent newer segment (for example, 1600–2000) by matching the mean of the first 20 years of the newer segment with the mean of the corresponding interval in the adjacent older segment (for example, 1600–1620). This nested approach allowed us to incorporate proxy records that do not span the full reconstruction interval.
All reconstructions except pairwise comparison (MATLAB) were performed in R (ref. 118).
For each proxy record included in the reconstructions, we used the ‘simulateBam’ function from the geoChronR package119 to calculate a 100-member banded age–depth model28, assuming a 1% counting error. We explicitly incorporated this chronological uncertainty and its influence on the variance structure of the reconstruction by calculating the 800-year ΔSLP reconstruction 200 times, at each iteration randomly sampling one realization from the age–depth model ensemble for each record. This was done separately for each combination of reconstruction method and gridded product used to calculate the observational ΔSLP target index (hereafter ‘target index’; eight reconstruction methods, three target indices, 200 age–depth model iterations = 4,800 ensemble members). This incorporates the probability distribution of the age–depth model ensembles, providing a robust treatment of age uncertainty.
To incorporate uncertainty arising from the possibility that some records have an outsized influence on the reconstruction, before each iteration, we randomly removed up to 15% of all possible contributing records.
To quantify uncertainty arising from the ΔSLP reconstruction method, we used eight different methods. These have various requirements for the input data—some require proxy records correlated with the target index (‘Data preparation’ section), whereas others use all available records. Records significantly correlated with the target indices in the calibration interval are denoted with a black outline in Extended Data Fig. 2.
PCA: Reconstructions based on PCA assume that the underlying gradient common to a group of time series significantly correlated with ΔSLP in the calibration window should be equivalent to ΔSLP (refs. 117,120,121). For PCA-based reconstructions, we therefore only used records that are correlated with ΔSLP in the calibration window.
For opPCA reconstructions, we performed PCA on the calibration window (that is, in which we know that the proxy records are correlated with the PWC) and then multiplied the loading of each proxy record on PC1 by the complete time series of the proxy records. The contribution of each record to the PCA was weighted according to the strength of its correlation with ΔSLP.
The direction of a PC axis is arbitrary. To align the temporal subsets, we flipped (if necessary) PC1 of the 1860–2000 subset to make it positively correlated with the target index and then aligned PC1 of subsequent temporal subsets to be positively correlated with their predecessors.
The fiPCA reconstructions were performed identically to the opPCA, except that PC1 was calculated over the full length of the proxy time series for each temporal subset.
CPS: The CPS method has been used in many multi-proxy palaeoclimate reconstructions24,122,123,124. In our implementation, all proxy records were scaled to unit variance and zero mean and then weighted according to their correlation with ΔSLP in the calibration window. The scaled and weighted records were composited and the composite was scaled to match the mean and variance of the ΔSLP in the calibration window. CPS reconstructions were performed using all available proxy records.
For CPSns, we repeated all steps for the CPS method but first filtering to only include records preserving an annually integrated signal (33 records; Supplementary Table 1). For Iso2k records, this determination was made on the basis of the ‘isotopeInterpretation1_seasonality’ metadata field23. For all other records, this was inferred from the primary publications. We only excluded records with a known (reported) seasonal bias.
For CPScoa, we repeated all steps for the CPS method but first filtering to only include records in a tropical Pacific ‘centre of action’, that is, only records between 40° S and 40° N, and 50° E and 50° W (42 records). This removes records with a higher potential for non-stationary teleconnections.
PCR: PCR is a multivariate regression method that has been used for palaeoclimate reconstructions of the past millennium44,125,126. PCR targets ΔSLP by performing PCA, calculating a linear regression of ΔSLP on the PCs and then retaining the minimum number of PCs required to maximize the correlation with ΔSLP. The number of retained PCs was determined using root mean squared error (RMSE) of prediction, estimated from cross-validation127. We chose the model with the fewest PCs that was still less than one sigma from the overall best model128. We performed PCR reconstructions using all proxy records (PCRall) and also on a subset that only included records significantly (P < 0.1) correlated with ΔSLP in the calibration window (PCRcor). We performed PCR reconstructions using the ‘pls’ R package128. Models were fitted to data in the calibration window and then values predicted for the full length of each temporal subset.
PaiCo: The non-linear PaiCo method was developed for use with multi-proxy palaeoclimate datasets129. The underpinning assumption of PaiCo is that an increase in a proxy record indicates an increase in the target index (ΔSLP) and the strength of agreement among proxy records on the change between two time points relates to the magnitude of reconstructed change in the target129. PaiCo reconstructions were performed in MATLAB, using records significantly (P < 0.1) correlated with ΔSLP in the calibration window.
The mean and variance of all reconstructed temporal subsets was adjusted so that:
For ease of comparison, we adjusted all reconstruction time series to match the mean and variance of ΔSLP calculated from HadSLP, although the results are not sensitive to this choice. When adjusting the variance of each reconstruction time series, we applied a single variance-scaling factor to the entire time series. That is, temporal variability in variance was maintained, potentially allowing for similar changes as seen in reconstructions of tropical Pacific SST130,131.
We repeated all reconstruction steps but with all correlations calculated on detrended datasets. This did not make any meaningful difference to the ΔSLP reconstruction ensemble, reconstruction skill or post-reconstruction analyses.
We calculated the following skill metrics for the reconstruction ensemble presented in the main text:
We performed skill tests on all 4,800 ensemble members, which are reported by reconstruction method and ΔSLP index (Extended Data Fig. 4a,b).
For ease of comparison with existing reconstructions of tropical Pacific variability, we calculated all skill metrics for the reconstruction median (Extended Data Table 1), as well as r for the median reconstructions for each reconstruction method and target index (Extended Data Fig. 3b). To estimate changes in reconstruction skill back through time, we calculated the same validation statistics for each temporal subset (Extended Data Fig. 5c).
To provide a minimum independent estimate of reconstruction skill, we calculated the same validation statistics across the 1900–1950 interval, using an otherwise exactly equivalent reconstruction ensemble calculated using a shorter calibration window (1951–2000) (Extended Data Fig. 5a,b and Extended Data Table 1). We also calculated the coefficient of efficiency for the reconstruction medians (Extended Data Table 1).
To assess internal consistency among reconstruction ensemble members, we considered all possible combinations of reconstruction method and ΔSLP training data and calculated the 30-year running correlation among each pair of ΔSLP time series (Extended Data Fig. 4c). When agreement is high among reconstruction ensemble members, this probably reflects a strong ΔSLP signal in the proxy datasets regardless of reconstruction method and target index choice.
To estimate the overall contribution of individual palaeoclimate records, we calculated the correlation of each component record (on its published chronology) with the ΔSLP reconstruction ensemble median across the interval to which that record contributed (Fig. 1b). Correlations were deemed significant if P < 0.05, and were calculated from the start of the earliest temporal segment to which each record contributed.
We calculated the full distribution of values in the 4,800-member ΔSLP reconstruction ensemble as well as for the preindustrial (1200–1849) and industrial-era (1850–2000) sections of the reconstruction (Fig. 2d). We performed two-sample Kolmogorov–Smirnov tests on the preindustrial versus industrial-era segments of all 4,800 individual ensemble members. We adjusted the P values to account for false discovery rate133. For 81% of ensemble members, the difference between the two time periods was not significant (P ≥ 0.05; Fig. 2d).
We calculated the temporal power spectrum for each ensemble member and determined frequencies at which each ensemble member has significant (P < 0.05) power. Our spectral analysis was based on the geoChronR (ref. 119) implementation of multitaper spectral analysis, by means of the ‘mtmPL’ function from the R package 'astrochron' (ref. 134). Significance of spectral peaks was established through a power-law null135. We report the proportion of ensemble members with a significant peak at each period below 75 years. Beyond 75 years, a maximum of 3% of ensemble members have significant power at lower frequencies (maximum n = 122 ensemble members, at period length 148 years). For comparison, we performed the same analysis on instrumental ΔSLP (Extended Data Fig. 7b).
To determine whether the industrial-era power spectrum is different from that of the preindustrial, we assessed the distribution of spectral power in only the most recent 150 years of the reconstruction (1850–2000) (Fig. 2c). To ensure a fair comparison with spectral densities in the preindustrial, we compared this with the distribution of spectral power in all possible 150-year periods before 1850 (Fig. 2b), still showing the proportion of ensemble members with power in each period.
To assess whether the power spectra are influenced by the ‘nesting’ reconstruction approach, we repeated the above analysis across the 1600–2000 interval, using a reconstruction ensemble derived only from proxy records with full coverage across that interval (otherwise identically constructed). In this way, we test (1) whether our nesting approach dampens low-frequency (decadal to multidecadal) variability and (2) whether differences between the power spectra of the preindustrial and industrial era are because of changing contributions from different proxy records (Extended Data Fig. 7).
To assess whether the 1992–2011 PWC strengthening13 is anomalous, we calculated the distribution of 20-year trends in the 4,800-member ΔSLP reconstruction ensemble, for comparison with the observed trend from 1992–2011. We provide the full distribution, as well as individual distributions for reconstructions trained on each gridded SLP product. The observed 1992–2011 trend is shown as a red bar on each distribution in Fig. 3b–e. ERA-20C data only go to 2010, so for the ERA-20C-only distribution, we show the distribution of 19-year trends.
To isolate potential influence of the 1991 Mount Pinatubo eruption on the 1992–2011 strengthening, we also calculated the distribution of 20-year trends that start in the year following volcanic eruptions equal to or greater in magnitude than the Mount Pinatubo eruption. We similarly compared the recent observed trends with these post-eruption distributions (Extended Data Fig. 8). We identified volcanic eruption years using global mean SAOD, a dimensionless metric for the scattering of solar radiation by aerosol particles, calculated in ref. 16 from the ‘eVolv2k’ ice-core reconstruction of volcanic sulfate aerosol loading35. Eruption years are defined as the maximum of each SAOD peak. Following findings from several recent studies36,136, we reassigned the year of the major Kuwae eruption to 1452 (as opposed to 1458 as per eVolv2k). The 1991 eruption of Mount Pinatubo had an estimated maximum SAOD of around 0.1.
For each ensemble member, we calculated the linear trend (regression coefficient) across two time intervals: 1200–1849 and 1850–2000. We show the distribution of trends in Extended Data Fig. 12; panel a shows the full distributions and panels b–d split the results according to the ΔSLP target index. We did not differentiate between significant and non-significant trends.
To assess the PWC response to volcanic forcing, we composited the ΔSLP response to all large volcanic eruptions intersecting the reconstruction interval. This technique, known as SEA, treats volcanic eruptions as replicate cases of the same process. This allows assessment of whether the PWC responds in a consistent manner to volcanic forcing. SEA is commonly used to assess the ENSO response to volcanic eruptions16,34.
For each time series (that is, each ΔSLP reconstruction ensemble member), we isolated 10-year segments spanning each eruption—3 years before and 6 years following each eruption. This resulted in n ten-year segments, in which n is the number of volcanic eruptions included in the SEA. We centred each 10-year segment according to its 3-year pre-eruption mean and then took the mean of all n segments. This provided a single 10-year composite time series, in which any consistent response in a particular year relative to the eruptions is concentrated and intrinsic variability should cancel out to an anomaly around zero. This replicates the SEA parameters of ref. 16, although our results are insensitive to the addition of several years either side.
We identified eruption years using SAOD as described in the ‘Calculating distribution of 20-year trends’ section. We restricted eruptions to those with SAOD ≥ 0.05. We performed SEA on all 4,800 ΔSLP reconstruction ensemble members and determined the significance of the results using the ‘double-bootstrap’ method of ref. 37. Specifically, we used the ‘random-bootstrapping’ approach, with confidence intervals generated from 1,000 pseudo-composite matrices. These pseudo-composites are also centred on the pre-eruption mean, resulting in relatively narrow confidence intervals before the eruption year. In Fig. 4 and Extended Data Figs. 9 and 11, we report the proportion of ensemble members with a significant (P < 0.05) positive or negative ΔSLP response to volcanic eruptions in each year of the analysis.
Twenty-five volcanic eruptions between 1203 and 1993 exceeded our 0.05 SAOD cutoff. We performed SEA using all 25 eruptions, as well as two sequences:
We also repeated sequence 2 but first removing the three most recent eruptions.
To directly compare the reconstructed and model-simulated tropical Pacific response to volcanic forcing, we replicated the analysis described in the previous section ('Assessing the PWC response to volcanic eruptions'), using data from all fully forced CESM1 LME members and eight PMIP3/4 models. ΔSLP calculated from climate models was scaled to match the variance of ΔSLP calculated from HadSLP. We performed SEA on three subsets of eruptions. In the first subset, we retained the same number of eruptions as input to the SEA performed on the ΔSLP reconstruction, that is, the 25 strongest eruptions during 1200–2000. We assessed two further subsets, the 12 strongest eruptions and the four strongest eruptions, allowing for comparison with the similar analysis performed in ref. 16, that is, Fig. 4B,C in that reference (although note that this reconstruction covered a different time interval). Eruption magnitudes were determined using the volcanic forcing reconstruction used to drive the model. For the CESM1 LME and PMIP3 models, this is ref. 53. For PMIP4 models, this is ref. 35.
We performed SEA on ΔSLP calculated from the PSL field of the atmospheric models, as well as relative SST (RSST) in the Niño 3.4 region (5° S–5° N, 170° W–120° W). RSST is the residual signal after removing mean tropical (20° N–20° S) SST anomalies from raw SST anomalies. We used RSST rather than raw SST anomalies because of the expectation that volcanic aerosols will cause cooling globally and mask the tropical Pacific response52,137. This allowed us to compare our findings with previous work investigating the effect of explosive volcanism on ENSO (in terms of SST anomalies), as well as comparing the oceanic and atmospheric responses over the tropical Pacific. We acknowledge that SEA is a suboptimal method for assessing the climatic response to explosive volcanism in climate models, which have full spatial and temporal data coverage and hence allow more nuanced analyses. However, performing the same analysis on model-derived and proxy-derived ΔSLP allows us to directly compare results.
We evaluated the relationship of the PWC with GMST by comparing our ΔSLP reconstruction ensemble with the PAGES 2k multi-proxy, multi-method ensemble (n = 7,000) reconstruction of GMST throughout the Common Era122. To assess temporal variability in the relationship between ΔSLP and GMST, we calculated correlations between the ΔSLP and GMST ensemble medians in many different time periods, starting between 1200 and 1990, spanning 10 to 800 years in duration (Fig. 6a).
We assessed uncertainty in the long-term relationship between ΔSLP and GMST by computing correlations between 4,000 unique combinations of individual members from both ensembles, over the full 1200–2000 interval (Fig. 6b).
We compared our ΔSLP reconstruction with published annually resolved reconstructions of tropical Pacific variability extending back to at least 1600 (Extended Data Fig. 6). Reconstructed climate modes include ENSO117,125,138,139,140,141,142 and the IPO31. ENSO reconstructions have different targets, for example, Niño 3, Niño 3.4 or ENSO indices incorporating several regions. If a study provided reconstructions of SST in several regions, we used the Niño 3.4 reconstruction. For the Last Millennium Reanalysis139, we used the Niño 3.4 reconstruction median. We clipped reconstructions to their common time period 1600–1978. Note that reconstructions have different reconstruction target seasons. We calculated 30-year running correlations between each ENSO reconstruction and the ΔSLP reconstruction median (Extended Data Fig. 6c), as well as correlations between all reconstructions across the 1600–1978 interval (Extended Data Fig. 6d). To compare ΔSLP with the IPO, we applied a 13-year Gaussian kernel low-pass filter to all ΔSLP ensemble members (following ref. 31) and then calculated the correlation of each smoothed ensemble member with the IPO reconstruction (1) over 1200–2000 and (2) only 1900–2000. For comparison, we correlated mean smoothed observational ΔSLP (from ERA-20C, ICOADS and HadSLP) with observed IPO variability over the 1900–2000 period (Extended Data Fig. 6b). In Extended Data Fig. 6b, we only show significant (P < 0.05) correlations.
The ensemble approach to this reconstruction allows estimation of reconstruction skill at several levels of detail. The simplest possible tests compare ensemble median reconstructed ΔSLP with mean ΔSLP from the three observational products (Extended Data Table 1). In this test, the reconstruction is highly correlated with observations. There is only a small difference in skill scores for tests on reconstructions using the entire calibration window (r = 0.81, P < 0.05) versus independent calibration-validation tests (r = 0.77, P < 0.05), whereby validation is performed on a 1900–1950 window, using reconstructions trained only on observational data from 1951–2000. The RMSE is low in both cases (0.27 for the full calibration window and 0.26 on the independent validation window). RE can range from negative infinity to one; reconstructions are generally considered skilful if RE >0。RE在所有情况下都是积极的。
比较亚气组中位数与目标指数（n = 3）和重建方法的独特组合（n = 8）揭示了与相关目标指数的相关性差异和亚集合中位数之间的一致性（扩展数据图3B）。对于所有三个目标指数，PCRALL子集合中位数与观测值高度高度相关。FIPCA亚汇编中位数始终与观测值相关。PAICO子集结中值通常显示与其他重建中值的最低相关性。
对于使用不同的培训指数计算的集合成员的技能得分分布之间的差异很小（扩展数据图4A），但重建方法之间的差异较大（扩展数据图4B）。如子安装中值（扩展数据图3B）所示，使用PAICO计算的集合成员倾向于表现最差，而使用基于PCR的方法计算的集合成员倾向于表现最佳。尽管基于PCA的方法具有最大的总体分布（偏斜到低分），但所有其他方法都有相似的中位数和四分位间范围。
当在独立窗口上计算技能得分（校准1951-2000，验证1900-1950）时，基于PCR的方法仍然表现最佳（扩展数据图5B），但是两种基于PCA的方法的表现最差，尤其是长的低分尾巴（扩展数据图5B）。
通过计算有助于重建的单个时间子集的技能得分（“所有方法共同的重建步骤”部分），我们可以估计随着时间的推移重建技能的变化（扩展数据图5C）。随着年龄的增长，技能降低了，这并不奇怪，鉴于代理数据的可用性从1600年左右迅速下降（扩展数据图1和2）。
通过比较CPS重建方法变体的技能得分，我们估算了（1）远离热带太平洋的代理档案的影响（因此很大程度上依赖于远程连接）和（2）具有已知偏见的代理人对特定季节的影响。当在整个1900–2000间隔中计算时，子集合中位数的技能得分（即使用CPS，CPSCOA和CPSN计算的重建集合成员中位数）非常相似（扩展数据表1，第一列）。在独立的1900-1950验证窗口（扩展数据表1，第四列）上计算时，技能得分之间存在较大差异。与其他两个变体中的任何一个相比，独立的CPSCOA重建具有更高的R和RMSE。这表明，远离热带太平洋的记录可能会对重建技能产生负面影响。但是，排除具有已知季节性偏见的记录并不能提高重建技能。
值得注意的是，CPSN的重建通常与其他方法的重建最少相似，通常具有更大的可变性，有时会显示出与其他方法重建相反符号的变化（图1A）。这可能是由于：（1）创纪录的季节性对重建的重大影响；（2）从特定档案馆（树纤维素）中丢失了许多记录；或（3）造成重建的记录数量减少。
扩展数据图4C证明了ΔSLP重建集合成员之间的一致性随着时间的流逝而变化，并且在大约1600年的集合一致性变化，与降低代理数据的可用性相吻合（扩展数据图1）。我们可以使用该协议来估计ΔSLP可以从该代理数据组合中恢复的程度。在1600年至2000年之间，合理的一致性表明，在此间隔内，ΔSLP信号强烈支持代理数据。在1600年之前，合奏成员之间的一致性较小，可能表明ΔSLP与某些代理记录之间关系中存在时间非平稳性。