Detection of oriented fractal scaling components in anisotropic two-dimensional trajectories

We propose a novel class of mixed fluctuations with different orientations and fractal scaling features as a model for anisotropic two-dimensional (2D) trajectories hypothesized to appear in complex systems. Furthermore, we develop the oriented fractal scaling component analysis (OFSCA) to decompose such mixed fluctuations into the original orientation components. In the OFSCA, the original orientations are detected based on the principle that the original angles are orthogonal to the angles with the minimum and maximum scaling exponents of the mixed fluctuations. In our approach, the angle-dependent scaling properties are estimated using the Savitzky–Golay-filter-based detrended moving-average analysis (DMA), which has a higher detrending order than the conventional moving-average-filter-based DMA. To illustrate the OFSCA, we demonstrate that the numerically generated time-series of mixed fractional Gaussian noise (fGn) processes with non-orthogonal orientations and different scaling exponents is successfully decomposed into the original fGn components. We demonstrate the existence of oriented components in the 2D trajectories by applying OFSCA to real-world time-series, such as human postural fluctuations during standing and seismic ground acceleration during the great 2011 Tohoku-oki earthquake.

Noise and fluctuations frequently display fractal-like scaling properties associated with long-range correlations in real-world complex systems, such as biological, geophysical, and economical systems 1,2 . The appearance of such fluctuations has been considered as a key marker associated with a universal principle hidden in the complex system dynamics 3 . Hence, various analysis methodologies have been developed to provide a detailed characterisation of such complex fluctuations. These methods have demonstrated that the scaling properties emerge not only in an auto-correlation of the time series [4][5][6] , but also in higher moment correlation (e.g. multi-fractality) [7][8][9] and in a cross-correlation between multivariate time series 10,11 .
The fractal scaling behaviour in two-dimensional (2D) trajectories is also of interest in real-world applications. For instance, human postural fluctuations during quiet standing 12 and animal movement patterns 13 have been modelled using 2D fractional Brownian motion (fBm) and fractional Gaussian noise (fGn). Most of the previous studies have characterised the 2D trajectories using the projections onto the orthogonal directions and otherwise under the assumption of 2D isotropic properties 12,13 . Such approaches cannot provide any information on geometrically anisotropic structure of the system. In this study, we hypothesise that real-world 2D trajectories display anisotropy that is different from conventional, isotropic 2D fBm and fGn 14 . Moreover, we propose a novel class of mixed fluctuations with different orientations and scaling features as a model for anisotropic 2D trajectories. In our model, the mixed fluctuations are assumed to consist of two independent components with non-orthogonal orientations as shown in Fig. 1. Such anisotropy has not yet been systematically investigated. The open problem here is to detect the hidden orientation components in the observed trajectories. To solve this problem, we develop the oriented fractal scaling component analysis (OFSCA) to characterise the directional properties of such trajectories and decompose the mixed fluctuations into the original components. In the OFSCA, the original orientations are detected based on the principle that the original angles are orthogonal to the angles with the minimum and maximum scaling exponents of the mixed fluctuations.

Scientific Reports
| (2020) 10:21892 | https://doi.org/10.1038/s41598-020-78807-z www.nature.com/scientificreports/ In our approach, we introduce an extension of the detrended moving-average analysis (DMA) 15 : the directional DMA (DDMA). The purpose of DDMA is to estimate the Hurst scaling exponent associated with the fBm-and fGn-like property in the projection of the observed time-series onto an axis forming angle θ with the abscissa axis. In real-world time-series analysis, it has been shown that non-stationary trends embedded in the time series harmfully affect the scaling exponent estimation and induce a misinterpretation of the correlation properties 16 . Therefore, scaling analysis methods such as detrending procedures have been widely used instead of conventional Fourier power spectral, rescaled range 4 , and structure-function 17 analyses. For instance, the practical options are the wavelet-decomposition-based method with a wavelet having vanishing moments 5 and detrending-operation-based scaling analysis methods, such as detrended fluctuation analysis (DFA) 16,18 and DMA 19 . Recently, we have established the theoretical foundation for DMA including higher-order DMA 6,15,20 and developed a fast implementation algorithm for DMA 21 to increase the reliability and applicability of DMA. In this study, we further extend the application of DMA to 2D trajectory analysis.
The remainder of this paper is organised as follows. First, we introduce the analysis procedure of anisotropic 2D trajectories. We illustrate our approach using a numerical sample, assuming a mixed fGn model displaying anisotropic long-range correlations. We then analyse real-world time-series, such as human postural fluctuations during standing and seismic ground acceleration during the great 2011 Tohoku-oki earthquake and discuss the interpretation of the oriented scaling components in the 2D trajectories. Finally, we provide a summary and some perspectives on future research.

Analysis method for anisotropic 2D trajectories
Let us consider a 2D time-series (trajectory) represented as the ordered series of points in the Cartesian coordinate plane, (x (1) [i], x (2) [i]) ( i = 1, 2, . . . N , where N is the length of the time series). If the time series is a sample path of 2D fGn given by two independent fGn sample paths, x (1) [i] and x (2) [i] , which are two series of fBm increments with the same Hurst exponent H, the auto-correlation properties are isotropic, independent of the orientation. Thus, the scaling property of each angular component (projection onto a rotated axis) is identical and not affected by any rotational transform. We introduce methods for detailed anisotropic characterisation assuming an anisotropic 2D time-series different from such 2D fGn. In our approach, we first evaluate the angledependent scaling properties using higher-order DMA 21 (henceforth, called DDMA) and then decompose the observed 2D time series into two components with different orientations and scaling properties.
Angle-dependent scaling analysis. The observed time series is projected onto an axis forming angle θ with the positive direction of the abscissa axis to detect the anisotropic scaling behaviour. Projected time series x (θ ) [i] is given as where θ is varied across the range of 0 ≤ θ < π.
Time series x (θ ) [i] for each θ is analysed using the Savitzky-Golay-filter-based DMA 6,22 . In the DDMA, we first calculate the integrated series of Anisotropic properties can be evaluated using the θ-dependent heterogeneity in F (θ ) (s) . In addition, the long-range correlation observed in each orientation, θ , is evaluated using the power-law increase of F (θ ) (s) as and quantified using scaling exponent α(θ) estimated as the slope of the double-logarithmic plot of F (θ ) (s) against s.
Note that F (θ ) (s) can be directly linked with auto-correlation function C (θ ) (k) and power spectrum 11,23 . That is, we have The analytical forms of kernels L(k, s) and G s (f ) 2 were demonstrated in 11,23 . Using these relations and assuming Eq. (4), we can show the scaling relations, In particular, when C (θ ) (k) decays exponentially to zero and S (θ ) (f ) shows the low-frequency plateau indicating short-term correlation, the scaling exponent results asymptotically in α = 0.5. In a previous study 11 , it was analytically shown that time-scale distortion between the scale in the time domain of DMA and the frequency in the Fourier spectral domain is induced by the higher-order DMA. Thus, we use the corrected time scale, s , instead of s. Although scale s in the zeroth-order DMA corresponds well to frequency f in the Fourier spectral domain, i.e., s = s/1.00 , s are given by s = s/1.93 and s = s/2.74 in the second-and fourthorder DMAs, respectively. Note that a straightforward implementation of the aforementioned procedure has a high computational complexity. Thus, a fast algorithm of DMA should be employed in the practical analysis 21 .
Decomposition of mixed long-range correlated fluctuations. We introduce a new class of 2D longrange correlated processes displaying special orientations and an orientation decomposition method for such processes. We illustrate our approach by assuming a mixed fGn model displaying anisotropic long-range correlations ( Fig. 1). Our model consists of two independent fGn processes {ǫ 1 [i]} and {ǫ 2 [i]} with two different Hurst exponents H 1 and H 2 ( H 1 > H 2 ), respectively ( Fig. 1a,b). It mixes these processes oriented at angles θ 1 and θ 2 to the x-axis (Fig. 1c,d) as Now, we consider the problem of decomposing the observation of Notably, the independent component analysis (ICA) is not properly applicable to our model, because both x (1) [i] and x (2) [i] in our model follow Gaussian distributions. Two linearly mixed independent Gaussian processes cannot be separated uniquely into two independent components (ICA for Gaussian processes is equivalent to principal component analysis).
Our approach is described as follows (Fig. 2). In this model [Eq. (6)], the projected time series is given by Fig. 2c).
That is, we obtain which is proportional to the original ǫ 2 [i] with H 2 . For the same reason, In contrast, when θ = θ 1 ± π/2 and θ = θ 2 ± π/2 , F (θ ) (s) shows a crossover of scaling exponents (see Fig. 2b). The forced linear fit to the broken lines in the log-log plot of F (θ ) (s) vs. s yields a slope in the range (H 1 , H 2 ) . Therefore, seeking two main orientations, θ min and θ max , respectively, with the minimum and maximum values of α(θ) , we can estimate the original orientations of ǫ 1 [i] and ǫ 2 [i] as θ 1 =θ min ± π/2 and

Real-world anisotropic scaling
We demonstrate the existence of anisotropic scaling by analysing real-world time-series, such as human postural fluctuations during standing (Figs. 4 and 5) and seismic ground acceleration during the great 2011 Tohoku-oki earthquake (Fig. 6).
Python source code for computing OFSCA and obtaining results of examples presented in the paper are public available for download in 24 .
Centre-of-pressure fluctuations of postural sway. We apply the proposed method to the human centre-of-pressure (CoP) trajectories in a standing posture 25 . The CoP is the projection of the centroid of the vertical force distribution on the ground plane. The CoP characteristics have been studied extensively to evaluate balance dysfunction for patients with neurological and motor disorders and to understand the postural control mechanism 26,27 .
To extract the oriented scaling components from the CoP characteristics, we analysed CoP trajectories measured in two standing posture conditions: quiet standing with eyes open (Fig. 4a) and a voluntary sway movement from the steady position to right-backward position with eyes open (Fig. 5a); the recording duration was 30 s with a sampling frequency of 100 Hz (for more details of the measurement, see 28,29 ). Recordings of CoP trajectories were obtained from a project approved by the Institutional Review Board of the National University of Physical Education and Sport of Ukraine. The procedures complied with the Declaration of Helsinki regarding human experimentation. Written informed consent was obtained from all the participants prior to the test. The angle dependence of F (θ ) (s) was evaluated over the range of 0 ≤ θ < π in increments of π/64 rad (Figs. 4b,c, 5b,c). We set the scaling range 1.2 < log 10s < 2.5 (from 0.16 to 3.2 s) and estimated the slopes of linear regressions (Fig. 4d) to find two representative orientations.
In the quiet standing condition, the minimum and maximum slopes were observed at θ min = 6 • and θ max = 124 • , respectively (Fig. 4d). Thus, the estimated orientations were θ 1 = 96 • and θ max = 34 • . As shown in Fig. 4d, the estimated minimum and maximum slopes (1.92 and 2.07, respectively) were almost equal. Nevertheless, the scaling behaviours of the decomposed components ǫ 1 [i] and ǫ 2 [i] (Fig. 4e) showed a nontrivial difference. As shown in Fig. 4f, a crossover point around log 10s = 2.2 was observed only in the component ǫ 2 with an orientation of 34 • . The observed scaling exponent α ≈ 2.0 indicates an fBm-like behaviour with longrange-correlated increments, which suggests a rather anomalously expanding and unstable CoP trajectory. In contrast, the smaller scaling exponent ( < 1.0 ) observed at the larger scales in the ǫ 2 orientation suggests a bounded and stable CoP trajectory. Therefore, the ǫ 2 orientation would correspond to the most stable direction of the posture control. This result is reasonable from an anatomical point of view, as a standing human body has more stability against left-right oscillations than against forward-backward oscillations while attempting to stand still. The ǫ 1 orientation almost coincided with the face-front (anteroposterior) direction, whereas the ǫ 2 orientation was considerably deviated from the shoulder-width (mediolateral) direction. This deviation originates from the motor asymmetry of the lower limbs. The mean CoP position of the subject was shifted slightly to the right from the coordinate centre of the platform, i.e. the subject body weight was distributed more on his dominant right leg than on his left one. Hence, we conclude that the subject controlled the postural balance mainly with his right foot 30 .
In the condition of voluntary sway movement, the subject was asked to perform swaying movements in the backward-right direction (Fig. 5a). In this condition, the minimum and maximum slopes were observed at θ min = 37 • and θ max = 132 • , respectively (Fig. 5d). Thus, the estimated orientations were θ 1 = 127 • and θ max = 42 • . The orientation with θ 1 almost coincided with the body sway direction, and the orientation with θ 2 was almost perpendicular to it. The scaling exponent of ǫ 1 along the body sway direction was 2.5 ( Fig. 5f), indicating a smoothly spreading fluctuation. In contrast, the scaling exponent of ǫ 2 at larger scales ( > 1 s) was 1.5 (Fig. 5f), indicating a Brownian-motion-like fluctuation.
The decomposed components ǫ 1 and ǫ 2 (Figs. 4e and 5e) appeared to be moving along each axis in an independent manner. Moreover, the cross-correlations between ǫ 1 and ǫ 2 were almost zero in both conditions. The Pearson cross-correlation coefficients were −0.07 and −0.12 , respectively, in the quiet standing and leaning motion conditions. Hence, our approach could detect and reconstruct the reasonable main orientations in the posture control dynamics. However, the kinetic and mathematical mechanism generating the observed fluctuations is still unclear. The source of the earthquake was not a single large earthquake but several interdependent earthquakes linked to each other. We analyse the recordings 32 of a seismic sensor located at Shin-Machi, Wakuya-cho, Toda-gun Miyagi Prefecture (Latitude: 38.5401 • N, Longitude: 141.1272 • E) to characterise the complicated earthquake process. From this location, the epicentre was located approximately 111 km east and the earthquake sources were distributed over the azimuth angle range of approximately 102 • to 140 • , where the azimuth angle φ is defined as the angle from the north direction counted clockwise, (i.e. φ ≈ π/2 − θ ). The analysed 2D time series is the acceleration on the ground surface with the north-south (NS) and east-west (EW) axes at the sampling frequency of 100 Hz. The ground acceleration has been considered the most important factor in determining the stress induced to structures during earthquakes. The top two panels in Fig. 6 show 3-min time series. This time series could be separated into four representative phases: P-, first S-, second S-, and aftershock waves (Fig. 6a,b). The previous studies 31 using inversion analysis have identified three main sources generating the representative waveforms. The first and second S-waves were generated, respectively, by the first and second main sources in the east-southeast direction from the measurement point. In addition, even in the aftershock phase, the earthquake wave induced www.nature.com/scientificreports/ by the third main source in the south-southeast direction was overlapped. However, the corresponding S-wave was not dominant because of the distance decay effect. In each phase, the angle dependence of F (θ ) (s) was estimated over the range of 0 ≤ θ < π in increments of π/64 rad using fourth-order DDMA. F (θ ) (s) and the local slopes are shown in Fig. 6c,d, together with the estimated representative orientations (red and blue dashed lines). We set the scaling range 1.3 < log 10s < 2.2 (from 0.20 to 1.6 s) and estimate the slopes in this range using linear regressions (Fig. 6e) to find two representative orientations.
A remarkable transition was observed in the orientations and the correlation properties before and after the arrival of the first S-wave. In the P-wave phase, the estimated orientations of ǫ 1 and ǫ 2 were θ 1 ≈ 85 • and θ 2 ≈ 170 • , respectively (the first column in Fig. 6 c,d). The corresponding slopes in the range of 1.3 < log 10 s < 2.2 were close to 0.5 and less than 1.0 (Fig. 6e). The fluctuation functions showed an evident difference between ǫ 1 and ǫ 2 (Fig. 6f). A change in the orientations was observed immediately after the arrival of the first S-wave (the second column in Fig. 6c,d). In the first and second S-waves, the orientations of ǫ 1 and ǫ 2 were θ 1 ≈ 135 • and θ 2 ≈ 45 • , respectively. In the aftershock wave, the orientations were rotated approximately 20 • in the counter clockwise (the fourth column in Fig. 6c,d). In addition, gradual increases of the slopes in the range of 1.3 < log 10s < 2.2 were observed through the earthquake process (Fig. 6f). The slopes approached 1.5 through and after the second S-wave phase.
All the plots of the fluctuation functions (Fig. 6f) showed a (non-increasing) plateau in the range of log 10s > 2.1 , indicating oscillating behaviour. To illustrate such behaviour, we also analysed the numerical time-series of a purely sinusoidal wave (Fig. 7a) and a second-order auto-regressive (AR(2)) model (Fig. 7b) using DMA. Here, the periods in the sinusoidal wave and the AR(2) model are set to 10 2.1 points. In Fig. 7d, the plot of log 10 F (θ ) (s) vs. log 10s shows a power-law increase in scales shorter than the oscillating period ( log 10s < 2.1 ) and a plateau in scales longer than the oscillating period. Notably, the power-law (scaling) exponents of the purely sinusoidal wave (Fig. 7d) cannot be linked with the power-law correlation properties and the exponent depends on the order of the DMA. In contrast, the AR(2) model (Fig. 7b) shows a short-term power-law correlation characterised by α = 2.5 (dashed lines in Fig. 7e) because of the power-law of the power spectrum as S(f ) ∼ f −4 for the high-frequency range. However, as shown in Fig. 7d, zeroth-order DMA cannot detect a scaling exponent larger than 2 6 . That is, zeroth-order DMA (commonly called DMA so far) cannot distinguish between a purely sinusoidal wave (Fig. 7a) and an AR(2) model. Therefore, the consistency of the scaling behaviour needs to be tested using higher-order DMAs to detect meaningful scaling behaviour. In our earthquake analysis (Fig. 7c,f), the scaling behaviour in the range of log 10s > 1.3 in higher-order DMAs showed consistent results.
The observed transition before and after the arrival of the first S-wave would reflect the mechanism difference. The P-wave is a compressional wave with no shear components and it travels faster than other seismic waves. In contrast, the S-wave is a shear wave and propagates faster than the P-wave. Although a detailed interpretation of our observation is not possible at present, our approach would provide a new method for earthquake analysis and advance the existing earthquake research activities, such as Refs. [33][34][35] .  (Fig. 6). The slopes ( ≈ 1.2 ) in the range smaller than the plateau regime characterise the correlation properties of the noisy fluctuations.

Discussion
We have proposed a sum of mixed fluctuations with different orientations and fractal scaling features as a model for anisotropic 2D trajectories. A decomposition method called OFSCA was developed to decompose a 2D trajectory into the original oriented components. Through the analysis of a sample anisotropic 2D fGn model, we have demonstrated that the OFSCA estimates the model parameters (the orientation angle of each component and its scaling parameter) successfully. Moreover, the analysis of real-world anisotropic 2D trajectories, such as human CoP and the great 2011 Tohoku-oki earthquake, has demonstrated the existence of their oriented fractal scaling components and provided detailed characterisation of their properties. These findings are merely exploratory, and further studies of the oriented scaling components will provide new insights on the complex fluctuations. Our observations on earthquake data analysis based on OFSCA are consistent with the findings of Ref. 31 , where the motion characteristics of the tectonic movement during the Tohoku earthquake were studied using inverted teleseismic P-wave data. The authors showed the relationship between slip distribution (tectonic movements) and strong ground motion while considering a long slip duration, large stress drop, and extensional (normal faulting) aftershocks. We consider that there could be a relationship between the Hurst exponent and strong ground motion, since as discussed in the Ref. 36 the power law in the earthquake statistics and non-linear nature of the earthquake propagation comes from the fractal nature of rough crack surfaces of crust and tectonic plate dynamics. We assume that OFSCA approach might provide basic information about the rupture process for real-time monitoring (real-time seismology). Further application of OFSCA to seismic data and the interpretation of the fluctuation characteristics of tectonic movements are required to understand the earthquakes more deeply and to predict their occurrence and propagation in the future.
Further extension of OFSCA is feasible for postural sway analysis, which has been used for the medical diagnosis and prognosis of patients with neurological and motor disorders. When a subject has a motor/sensory disorder (such as hemiplegia or incomplete spinal cord injury), which has an asymmetrical effect on the standing stability 37 , the evaluation of the directional dependency of the CoP trajectory could provide useful information on the stabilising ability and dysfunction.
The proposed OFSCA can be used to characterise other real-world 2D trajectories. For instance, in Ref. 38 , an automated continuous monitoring system for tracking the behaviour of pigs was used to monitor their movements with a depth camera. Then, using 2D trajectories of pig movement, the directional behaviour and characteristics of the locomotion dynamics could be quantified with OFSCA and mapped onto the animal health issues to detect warning signs. This approach can be used to characterise the trajectories of other animals, such as fish, cows, and unicellular organisms. It can also be applied in several other behavioural, ecological, and veterinarian studies, and to amend the existing techniques of 2D and 3D trajectory analyses 39,40 .
Another promising application is the analysis of the 2D trajectories of individuals encountered in group games, such as soccer, in both humans 41 and robots 42 . The group dynamics and interconnections of the player's movement dynamics may provide interesting insights for sports studies, and the orientation of the fluctuation characteristics can be a valuable feature of single, group, and event-specific trajectories 43 . Furthermore, the OFSCA approach can be extended to higher-dimensional data analysis, e.g. studying complex evolutionary dynamics 44 , or understanding the structure and patterns in high-dimensional RNA sequencing datasets 45 . Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.