Visual analysis method for shared bicycle demand mode

Abstract

The invention relates to a visual analysis method for a shared bicycle demand mode. The method comprises the steps of obtaining shared bicycle data and performing data cleaning; dividing the area of the target region into a plurality of target sub-grid regions; identifying specific sites in the cleaned shared bicycle data set, aggregating data categories of the cleaned shared bicycle data set according to functions and sites and grid information data, calculating a weather condition quantized value by adopting a space weighted average method, obtaining shared bicycle data feature representations Embedingear and Embedingday, then performing clustering and dimension reduction, and calculating a weather condition quantized value; and finally, visually displaying the result of each step. According to the method, a deep learning method and a visualization and visual analysis technology are combined, a visual analysis framework is provided, the framework can effectively calculate the shared bicycle demand mode through the deep learning method, and a user is helped to understand and relieve the defect that deep learning is poor in interpretability through visualization.

Description

Technical Field

[0001] This invention relates to time series feature representation methods and multi-source spatiotemporal data analysis methods, and particularly to a visual analysis method for shared bicycle demand patterns. Background Technology

[0002] As an environmentally friendly and efficient mode of transportation, shared bicycles have played a significant role in alleviating urban traffic congestion and promoting low-carbon travel. However, with the continuous growth in the number of shared bicycles deployed, numerous problems have gradually emerged in areas such as operation and management, user experience, and urban management. For example, some areas have experienced severe bicycle shortages, leading to a situation where users cannot find a bicycle and their travel needs cannot be met; in other areas, an oversupply of bicycles has resulted in indiscriminate parking. Existing research indicates that accurately identifying the evolution patterns of shared bicycle demand can help urban traffic managers and operators understand the dynamic changes in shared bicycle demand under different time and spatial conditions, providing them with effective guidance. Therefore, how to mine demand evolution patterns from massive, multi-dimensional shared bicycle data and find effective analytical methods is currently one of the focal points of attention in academia and the transportation field.

[0003] As a short-distance commuting tool, the demand for shared bicycles exhibits a typical cyclical pattern over time, and time series analysis can effectively capture these dynamic patterns. Contrastive learning techniques are widely used in time series analysis research, generating multiple sub-series views through data augmentation. This allows for the discovery of hidden patterns in time series data without the need for labeled data, resulting in high-quality feature representations. However, existing contrastive learning techniques often rely on random data augmentation, which can easily introduce inappropriate inductive biases when processing data with significant temporal characteristics, such as shared bicycle data. For example, if random pruning is applied to shared bicycle data, and more data is pruned during winter months, the model may be more inclined to capture winter demand characteristics, ignoring features from other seasons. Therefore, to address these issues, more adaptive feature representation methods are needed to obtain higher-quality feature representations.

[0004] Furthermore, the demand for shared bicycles is also influenced by various external spatiotemporal factors such as local weather and points of interest (POIs). Relying solely on time-series features extracted by models is insufficient to fully understand the evolution of demand, and the high-dimensional feature representations obtained by deep learning models lack direct interpretability. Visualization and visual analytics technologies, through interactive charts, help users intuitively explore and understand data patterns, effectively alleviating the problem of insufficient interpretability of feature representations. However, existing methods suffer from information overload and visual redundancy when visualizing multi-source heterogeneous spatiotemporal data of shared bicycles and their influencing factors, making it difficult to clearly demonstrate the dynamic relationship between demand changes and different factors. Therefore, it is necessary to construct a joint analysis framework that integrates deep learning feature representations and interactive visualization technologies to achieve a concise and effective visual analytics system. Summary of the Invention

[0005] To address the aforementioned problems in existing technologies, this invention utilizes visualization and visual analytics technologies, integrating deep learning-based feature representation methods, to propose a visual analytics framework for shared bicycle demand. This framework aims to help users uncover and analyze the evolution patterns of shared bicycle demand.

[0006] To solve the above-mentioned technical problems, the present invention adopts the following technical solution: a visual analysis method for shared bicycle demand patterns, comprising the following steps:

[0007] S1 acquires shared bicycle data and performs data cleaning.

[0008] S2 divides the area of ​​the target region into multiple target sub-grid regions.

[0009] S3 identifies specific stations in the cleaned shared bicycle dataset: The DBSCAN clustering algorithm is used to identify the precise location of each station in the dataset. First, the list of latitude and longitude coordinates for each station is traversed and compiled. Then, the DBSCAN algorithm is used to cluster the latitude and longitude coordinates of each station. The neighborhood radius is set to Q meters, and the minimum number of samples in the neighborhood is half the number of points in the list. Finally, the center coordinates of the largest cluster are used as the unique latitude and longitude coordinates of that station.

[0010] S4 aggregates the cleaned shared bicycle dataset into five categories based on their functions: life services, public services, housing and land, transportation, and leisure and tourism.

[0011] S5 Site and Grid Information Data Aggregation: The cleaned shared bike data is stored as Parquet files, partitioned by month according to the time dimension, from both the site level and the target sub-grid area level. Two new OD tables are output: one OD table records detailed order information from the origin station to the destination station for each hour, and the other OD table records detailed order information from the origin target sub-grid area layer to the destination target sub-grid area layer for each hour.

[0012] S6 uses a spatial weighted average method, taking the geographic center of each target subgrid area as the reference point, constructing an inverse distance weight based on the spatial distance between each meteorological station and the reference point, and weighting temperature and wind speed to calculate the quantitative value of weather conditions.

[0013] S7 Shared Bike Data Feature Representation: Shared bike data from different target sub-grid regions are represented using both annual and daily time granularities. In the annual analysis, the sample size N represents the effective target sub-grid regions obtained in step 2, and the time step T is 365. Four features are selected as feature dimensions F: daily demand, average riding time, average riding distance, and peak-hour demand percentage. After data standardization, the model input X∈R is obtained. N×T×F The shape is N×365×4, and the feature extraction process is represented as Embedding = Model(X). The feature representation of each data sample represents the 320-dimensional feature representation of that sample in the embedding space, i.e., the shape of the Embedding is N×320, serving as the annual scale feature representation. year .

[0014] The daily scale represents the demand evolution pattern of the same target sub-grid area on different dates, where N' represents the number of different dates (365), T' is 24 hours, and the feature dimension F' includes five features: hourly demand, average cycling time and average cycling distance, whether it is a holiday, and the weather index of the day. The resulting model input X∈R N′×T′×F′ After the same feature extraction process as on an annual scale, the daily scale feature representation (Embedding) is obtained. day .

[0015] S8 Clustering and Dimensionality Reduction: Embedding the feature representation obtained in step 7 year and Embedding day K-means clustering was used for cluster analysis, and the optimal number of clusters was determined using the elbow method to identify the macroscopic and microscopic demand patterns of different target subgrid regions. Clusters were divided into 5 classes at the annual scale and 3 classes at the daily scale. Principal component analysis was then used to analyze the feature representation obtained in step 8. year and Embeddingday Dimensionality reduction is performed to obtain its two-dimensional projection.

[0016] S9 visualizes all results from S1 to S8; at the annual scale, they are clustered into five categories: high-demand mode, medium-demand mode, low-demand mode, anomalous mode 1, and anomalous mode 2. The visualized clustering diagram shows that the high-demand mode has higher ridership and peak-hour percentages than the other two modes, while its riding distance and duration are similar to the medium-demand mode. The low-demand mode has similar ridership to the medium-demand mode, but its riding distance and duration are significantly higher than the other two modes. Anomalous modes 1 and 2 are two outlier clusters belonging to anomalous target sub-grid regions. Observing their annual demand changes and information summary diagrams reveals that the annual demand and number of stations in these anomalous target sub-grid regions are both less than the set threshold; therefore, these stations are removed.

[0017] The three daily-scale patterns were named Category 0, Category 1, and Category 2, respectively. Observing their visualized feature radar charts revealed that Category 0 had higher cycling distance and duration than the other two categories, Category 2 had higher number of rides than the other two categories, while Category 1, except for the number of rides, had similar metrics to Category 2 in all three aspects.

[0018] Furthermore, the data cleaning process in S1 involves filtering data that meets one of the following three conditions:

[0019] 1) Cycling time less than 1 minute;

[0020] 2) The cycling distance is less than 100 meters;

[0021] 3) The starting site contains abnormal sites. Sites with fewer than 50 records per year are defined as abnormal sites.

[0022] Furthermore, the calculation process in S2 of dividing the area of ​​the target region into multiple target sub-grid regions is as follows:

[0023]

[0024] Where, N col and N row These represent the number of columns and rows that can be obtained, respectively. 111.32 means that each degree of latitude and longitude is approximately 111.32 kilometers. lon and lat represent latitude and longitude, respectively.

[0025] Furthermore, the weather index calculation includes the following specific steps:

[0026] Step 6-1: Calculate the center coordinates of the target subgrid region, and then use the semi-versus formula to calculate the distance d from each weather station to the center of the target subgrid region:

[0027]

[0028]

[0029] Among them [center] lat ,center lon ] represents the center coordinates of the target subgrid region, d represents the spherical distance from the meteorological station to the center coordinates of the target subgrid region, r represents the radius of the sphere, lon and lat represent the latitude and longitude of the meteorological station, respectively. min and lat max These represent the minimum and maximum latitude values ​​of the target subgrid region, respectively. min and lon max These represent the minimum and maximum longitude values ​​of the target sub-grid region, respectively. Next, the inverse distance weights for each meteorological station are normalized. For M meteorological stations, the normalized inverse distance weight α for the i-th meteorological station... i The weighting value for weather stations is calculated as follows:

[0030]

[0031] Where p is the distance attenuation parameter;

[0032] Step 6-2: Calculate the weighted meteorological index S for each meteorological index in the target sub-grid area. agg Then perform normalization:

[0033]

[0034] Among them, S i Let S represent the meteorological index of the i-th weather station. min With S max Let represent the minimum and maximum values ​​of the weighted meteorological indicators, respectively, and V be the normalized result. Next, the Pearson correlation coefficient is used to measure the correlation between each meteorological indicator and the demand for shared bicycles; the calculation process is as follows:

[0035]

[0036] Among them, S i and Y i These represent meteorological indicators and the demand for shared bicycles, respectively. and These are their mean values. After calculation, the absolute values ​​of the correlation coefficients between the three meteorological indicators—temperature, wind speed, and precipitation—and the demand are normalized and used as the weight values ​​for these three meteorological indicators.

[0037] Step 6-3: Weather events are identified using binary numbers, with each bit representing an event: 1 indicates occurrence, and 0 indicates non-occurrence. First, the Pearson correlation coefficient is used to measure the correlation between different weather events and demand. This coefficient is then normalized and used as the weight for each weather event. Next, a score is calculated for each weather event. event :

[0038]

[0039] Among them, F ie Let β represent weather event e at the i-th weather station. e This indicates the weighting percentage of weather events.

[0040] Step 6-4: Set the weather index weighting of meteorological indicators to 0.9 and the weather index weighting of weather events to 0.1. The final weather index calculation process is as follows:

[0041]

[0042] Where K represents the number of meteorological indicators, λ j and V j Let I represent the weight percentage and normalized value of the j-th meteorological indicator obtained in step 6-2, respectively. weather The value is between 0 and 1, when I weather When the value is ≥0.7, the weather condition is considered excellent; when 0.2≤I weather When I < 0.7, the weather conditions are considered good; when I weather When the value is less than 0.2, the weather conditions are considered bad.

[0043] Furthermore, obtaining the shared bicycle data feature representation in S7 includes the following specific steps:

[0044] S7-1 constructs training data from both annual and daily scales. For the annual scale, the OD table obtained in S5 is used to calculate four features for each target sub-grid region: daily demand, average riding time, average riding distance, and peak-hour demand percentage. The feature dimension of each annual time series data sample is 365×4. All annual time series data samples are randomly divided into a 70% annual training set and a 30% annual validation set. For the daily scale, the OD table obtained in S5 is used to calculate five features for each target sub-grid region: hourly demand, average riding time, average riding distance, whether it is a holiday, and the weather index obtained in S6. The feature dimension of each daily time series data sample is 24×5. All daily time series data samples are randomly divided into a 70% daily training set and a 30% daily validation set.

[0045] The S7-2 TimeCA model is constructed in three parts: Part 1: data augmentation by randomly pruning and Gaussian perturbation to enhance time series data samples; Part 2: autoencoder to encode the augmented time series data; and Part 3: calculation of the joint loss function based on the output of the autoencoder.

[0046] Part 1 Data Augmentation: First, for each time series data sample X = {x1, x2, ..., x...} of sequence length L... L Random cropping is performed. In the annual time series data sample, X = {x1, x2, ... x} i ,...,x 365}, where x i Represents the characteristics of day i; in a daily-scale time series data sample, X = X′ = {x′1, x′2, ..., x′} j ,...,x′ 24},x′ j This represents the characteristics of the j-th hour of the day;

[0047] Set the length of the overlapping sequence of the cropped samples to l, l∈[2] t_unit [,L], where t_unit is a hyperparameter. Randomly select the starting position t2 for pruning, where t2∈[0,Ll], then t3=t2+l. Next, randomly select two time points t1∈[0,t2-1] and t4∈[t3+1,T-1] to obtain two subsequences x with repeating segments. left and x right These two subsequences generated from the same sample are defined as positive sample pairs.

[0048] The above method was used to prune each annual and daily time series data sample, resulting in two subsequences x with repeating segments for each sample. left and x right ;

[0049] Then, independent Gaussian perturbations are applied to different feature channels in the subsequence. For each feature channel k and timestamp t of the sequence, the specific perturbation process is as follows:

[0050]

[0051] in, It is a Gaussian perturbation, σ 2 This represents the variance of the noise. After the above perturbation, each sample from the yearly or daily scale will result in two enhanced subsequences x′. left and x′ rightThese enhanced subsequences serve as input to subsequent autoencoders;

[0052] Partial autoencoder: consists of an encoder and a decoder. The encoder extracts key features from augmented time-series data samples, and the decoder restores and reconstructs the original time-series data samples based on the encoder output;

[0053] The encoder consists of n stacked residual modules, each containing two dilated convolutional layers: the first layer captures dependencies over a shorter time span, and the second layer expands the receptive field to capture long-term dependencies in the time series. After each convolutional operation, the residual modules are normalized, and the output of each convolutional layer is non-linearly mapped using a non-linear activation function. Then, the input and output are summed through residual connections to accelerate the convergence process and obtain the feature representation. The i-th residual module, ResBlock... i The information transmission process is represented as follows:

[0054]

[0055] Where f represents the normalization and nonlinear activation functions in the residual module, and g represents the convolution operation. and Let these represent the weight matrices of the first and second convolutional layers in the i-th residual module, respectively. and is the dilation coefficient of the two dilated convolution layers in the i-th residual module. The encoder output r, obtained through n stacked residual modules, is used as the input to the decoder.

[0056] The decoder consists of n stacked convolutional modules. Each module contains a dilated convolutional layer and a one-dimensional convolutional layer. First, the dilated convolutional layer matches the encoder's receptive field to recover patterns over a longer time range in the time series. Then, the one-dimensional convolutional layer makes the reconstructed time series more closely resemble the original sequence data in terms of local details. The decoder output x is obtained after passing through the n convolutional modules. Recon ;

[0057] Partial three-way joint loss function: This includes soft contrast loss and reconstruction loss, which jointly constrain the model training process. The soft contrast loss is calculated separately at the instance layer and the temporal layer.

[0058] At the instance level, for each pair of original time series data samples Q i Q j After min-max normalization, the distance is calculated, and the soft-assigned weights w are defined using the sigmoid function. I (i,j):

[0059] w I (i,j)=2α·σ(-τI ·D(Q i Q j ))#

[0060] Where σ is the sigmoid function, which converts the distance into a weight between 0 and 1, and τ I α is a hyperparameter controlling the sharpness of soft assignment. D is a dynamic time programming algorithm used to measure the similarity between two time series data samples. α controls the weights of these similar but different time series data sample pairs. Next, the similarity p between different time series data sample pairs is calculated. I :

[0061]

[0062] Where, r i,t and r j,t Q represents i Q j The characteristic representation of the encoder output timestamp t.

[0063] The loss of the comparison between the i-th time series sample and other samples at time stamp t in the instance layer can be expressed as:

[0064]

[0065] Soft instance contrast loss The first term is the loss for positive sample pairs, where i+N represents the positive sample corresponding to the i-th time series sample, and the second term is the loss for other time series sample pairs;

[0066] Perform the above-described operations at the instance level on annual-scale time series data samples and daily-scale time series data samples respectively;

[0067] At the time level, soft allocation is calculated using the difference between timestamps, w T The calculation process is shown below:

[0068] w T (t,t′)=2·σ(-τ T ·|tt′|)#

[0069] Where |tt′| represents the difference between two timestamps t and t′, τ T It is a hyperparameter that controls the soft-assignment sharpness. According to w T (t,t′), the contrast loss at the time level can be expressed as:

[0070]

[0071] Where, r i,t and r i,t′p is a feature representation of the same time series at different timestamps t and t′. T (i,(t,t′)) represents the similarity of the same time series at different timestamps. Soft temporal contrast loss. The first term represents the loss of the positive sample pair, T represents the time dimension length of a single time series sample, t+T represents the timestamp of the positive sample corresponding to time t, and the second term is the loss of other sample pairs;

[0072] The reconstruction loss is calculated by measuring the Euclidean distance between the reconstructed sequence output by the decoder and the corresponding original sequence:

[0073] l Recon =MSE(Q,x Recon )#

[0074] Where MSE is the mean squared error function, Q is the original time series sample, and x Recon This is the corresponding reconstructed sequence output by the decoder module. Finally, the joint loss function is expressed as:

[0075]

[0076] Where λ is the balance coefficient of the reconstruction loss function, 2N represents the sample size after data augmentation, 2T represents the total number of candidate timestamps for a sample to participate in the time layer comparison after data augmentation, and the loss calculated during the training phase is called the training loss.

[0077] Perform the above-described time-level operations on annual-scale and daily-scale time-series data samples respectively;

[0078] S7-3 TimeCA Model Training: A balance coefficient λ is set, and the Adam optimizer is used to complete the optimization process. When the reconstruction loss function value no longer changes, a conditioned TimeCA model is obtained. After training, data augmentation is performed on the validation set samples in the same manner as during training. The encoder extracts features and performs joint loss calculation to obtain the validation loss. The relative error between the training loss and the validation loss should not exceed a preset threshold δ.

[0079] S7-4 uses a pre-trained TimeCA model for feature extraction. After training, the encoder module is used for feature extraction of the original sequence, outputting a feature representation vector. The original annual-scale time series data outputs a 320-dimensional feature representation, i.e., the embedding. year The shape is N×320; the original diurnal time series data output diurnal feature representation embedding. day The shape is N×20.

[0080] Furthermore, the S9 visualization display includes five parts: an overview analysis module, a feature analysis module, a pattern analysis module, a comparative analysis module, and a console.

[0081] Overview Analysis Module: Includes an overview map of regional POI information that displays the five major categories aggregated in S4, and a spatiotemporal overview map of demand information that displays the clustering results at the S8-year scale.

[0082] Feature Analysis Module: Used to display line graphs showing significant fluctuations in demand across target subgrid areas at different time granularities, and region-site feature analysis graphs showing demand changes between different subgrid areas.

[0083] Pattern Analysis Module: Includes a regional macro-pattern clustering map for dimensionality reduction display of clustering results at the S8-year scale, a pattern feature radar chart for multi-dimensional display of statistical characteristics of different categories, and a regional pattern evolution map for displaying different demand patterns of the same regional grid on different dates.

[0084] Comparative Analysis Module: Contains regional demand comparison charts to display regional demand patterns for different regions, the same regional grid on different dates, and different regional grids.

[0085] Compared with the prior art, the present invention has at least the following advantages:

[0086] 1. TimeCA employs a hybrid data augmentation strategy to mitigate inductive bias caused by random augmentation. Traditional contrastive learning models rely on a single random pruning data augmentation strategy, which is prone to introducing inappropriate inductive bias due to the randomness of pruning, thus disrupting the global regularity of the original time series data. TimeCA uses a hybrid data augmentation strategy of "random pruning + Gaussian perturbation," which can effectively mitigate the negative impact of inductive bias on feature quality and obtain higher-quality positive and negative samples.

[0087] 2. TimeCA employs a joint loss function to constrain the model training process. Standard contrastive learning techniques often overlook the inherent correlation between similar time series instances and adjacent timestamp values. Therefore, TimeCA uses a soft contrastive learning strategy to calculate the contrastive learning losses at the instance and time layers. The instance layer uses dynamic time programming to calculate sample distances and assigns soft weights, while the time layer assigns soft weights based on timestamp differences, flexibly handling the similarity and differences between samples. Furthermore, it uses the Euclidean distance similarity between the reconstructed sequence and the original sequence as the reconstruction loss, ensuring that the reconstructed sequence is as similar as possible to the unenhanced time series, resulting in higher quality feature representations of the time series data learned by the model.

[0088] 3. A quantitative method for weather index was proposed, which calculates quantitative values ​​by weighting various meteorological indicator monitoring data from multiple meteorological stations, allowing users to intuitively judge weather conditions. Furthermore, it integrates shared bicycle data and POI data to better understand the relationship between weather conditions and POI demand for shared bicycles.

[0089] 4. This invention designs a complete visual analytics workflow, combining deep learning methods with visualization and visual analytics technologies. The proposed visual analytics framework mainly includes four modules: data preparation, data processing, feature extraction and analysis, and visualization and visual analytics, forming a general framework for shared bicycle demand analysis. This framework can not only effectively calculate shared bicycle demand patterns through deep learning methods, but also help users understand through visualization, mitigating the drawback of poor interpretability in deep learning.

[0090] 5. This invention designs and implements a visual analysis system for shared bicycle demand analysis. This system can help users efficiently analyze multi-source heterogeneous spatiotemporal data composed of shared bicycles and related weather and POI data. Through a concise and novel view and interaction design, it alleviates problems such as information overload, visual redundancy, and the lack of direct interpretability of high-dimensional feature representations obtained by deep learning models. It can improve the efficiency of shared bicycle management and also play a role in its decision-making and scheduling. Attached Figure Description

[0091] Figure 1 This is a system display diagram of the present invention.

[0092] Figure 2 This is a diagram of the TimeCA model architecture proposed in this invention.

[0093] Figure 3 This is a flowchart of the present invention.

[0094] Figure 4 This describes the evolution of patterns in regions with different demand patterns.

[0095] Figure 5 This represents the model evaluation results.

[0096] Figure 6 The results are from the model ablation experiment.

[0097] Figure 7 This is the result of user research for this invention.

[0098] Figure 8 This is a summary diagram of the annual demand changes and information for areas with abnormal patterns.

[0099] Figure 9 This is an information graph of the characteristic indicators for evolutionary pattern categories 0, 1, and 2. Detailed Implementation

[0100] The present invention will now be described in further detail.

[0101] A visual analysis method for shared bicycle demand patterns includes the following steps:

[0102] S1 Data Cleaning: Shared bike data is typically collected from devices such as smart locks and user mobile phones, which may lead to problems such as weak GPS signals and high network latency, resulting in data loss or anomalies. This involves acquiring shared bike data and performing data cleaning.

[0103] S2 divides the area of ​​the target region into multiple target sub-grid regions. Region grid division: In order to capture the spatial distribution of shared bicycle demand in more detail, the present invention divides the target region into target sub-grid regions with an area of ​​1km×1km.

[0104] S3 identifies specific stations in the cleaned shared bicycle dataset: In the shared bicycle dataset, the latitude and longitude coordinates of the same station vary across different order data. Therefore, the DBSCAN clustering algorithm is used to identify the accurate location of each station in the shared bicycle dataset. First, the list of latitude and longitude coordinates for each station is traversed and counted. Then, the DBSCAN algorithm is used to cluster the latitude and longitude coordinates of each station. The neighborhood radius is set to Q meters [Q = 200 meters], and the minimum number of samples in the neighborhood is half the number of points in the list. Finally, the center coordinates of the largest cluster are used as the unique latitude and longitude coordinates of that station.

[0105] S4 aggregates the cleaned shared bicycle dataset into five categories based on their functions: life services, public services, housing and land, transportation, and leisure and tourism.

[0106] POI Data Category Aggregation: In order to more efficiently identify the dominant functions of grid areas, the 14 categories of POI data provided by OSM are aggregated into 5 major categories according to their functional categories.

[0107] S5 Site and Grid Information Data Aggregation: The cleaned shared bicycle data is stored in Parquet files, partitioned by month according to the time dimension, from both the site level and the target sub-grid area level. This file format uses columnar storage, which can effectively support the management and analysis of large-scale shared bicycle data. Two new OD (Origin-Destination) tables are output. One OD table records detailed order information from the origin site to the destination site for each hour, and the other OD table records detailed order information from the origin target sub-grid area layer to the destination target sub-grid area layer for each hour.

[0108] S6 Weather Index Calculation: The meteorological dataset contains observation data from multiple different stations. The spatial weighted average method is used, with the geographic center of each target subgrid area as the reference point. Based on the spatial distance between each meteorological station and the reference point, an inverse distance weight is constructed, and temperature and wind speed are weighted to calculate the quantitative value of the weather condition.

[0109] S7 Shared Bike Data Feature Representation: To comprehensively reveal the demand patterns of shared bikes and meet the decision-making needs of long-term urban macro-planning and daily operation management, shared bike data from different target sub-grid areas are represented using features at two time granularities: annual and daily. In the annual-scale analysis, the sample size N is the effective target sub-grid area obtained in step 2, and the time step T is 365. To balance the long-term spatial needs of urban transportation planning with the short-term response of operation management to time-slot scheduling, four features are selected as feature dimensions F: daily demand, average riding time, average riding distance, and the proportion of demand during peak hours (8-10 am and 5-7 pm). After data standardization, the model input X∈R is obtained. N×T×F The shape of the embedding is (N×365×4), and the feature extraction process is represented as Embedding = Model(X). The feature representation of each data sample represents the 320-dimensional feature representation of that sample in the embedding space, i.e., the shape of the embedding is (N×320), serving as the annual-scale feature representation. year ;

[0110] The daily scale represents the demand evolution pattern of the same target sub-grid area on different dates, where N' represents the number of different dates (365), T' is 24 hours, and the feature dimension F' includes five features: hourly demand, average cycling time and average cycling distance, whether it is a holiday, and the weather index of the day. The resulting model input X∈R N′×T′×F′ After the same feature extraction process as on an annual scale, the daily scale feature representation (Embedding) is obtained. day ;

[0111] S8 Clustering and Dimensionality Reduction: Embedding the feature representation obtained in step 7 year and Embedding day K-means clustering was used for cluster analysis, and the optimal number of clusters was determined using the elbow method to identify the macroscopic and microscopic demand patterns of different target subgrid regions. Clusters were divided into 5 classes at the annual scale and 3 classes at the daily scale. Principal component analysis was then used to analyze the feature representation obtained in step 8. year and Embedding day Dimension reduction is performed to obtain its two-dimensional projection;

[0112] S9 visualizes all the results from S1 to S8;

[0113] At the annual scale, the demand is clustered into five categories: high-demand, medium-demand, low-demand, anomalous pattern 1, and anomalous pattern 2. The visualized clustering diagram shows that the high-demand pattern has a higher number of rides and a higher proportion of peak hours compared to the other two patterns, while its riding distance and duration are similar to the medium-demand pattern. The low-demand pattern has a similar number of rides to the medium-demand pattern, but its riding distance and duration are significantly higher than the other two patterns. Anomalous patterns 1 and 2 are two outlier clusters, belonging to the anomalous target sub-grid region. Their annual demand changes and information summary diagrams are observed. Figure 8 If it is found that the annual demand and number of sites in these abnormal target sub-grid areas are less than the set threshold, then such sites will be removed in order to save on operating and management costs.

[0114] The clustering in S8 classifies the annual scale into 5 patterns and the daily scale into 3 patterns. This is achieved by observing the regional macroscopic pattern clustering diagrams in the visualization system. Figure 1 -D1) reveals that the three main clusters at the annual scale cover most of the samples, which, combined with the pattern feature radar map ( Figure 1 -D2) can be used to define these three main clusters as high-demand, medium-demand, and low-demand modes: the high-demand mode is significantly higher than the other two modes in terms of the number of rides and the proportion of peak hours, while the riding distance and riding time are similar to the medium-demand mode. The low-demand mode is not much different from the medium-demand mode in terms of the number of rides, but its riding distance and riding time are significantly higher than the other two modes. The other two outlier clusters are defined as anomalous mode 1 and anomalous mode 2, respectively. Observe the annual demand changes and information summary diagrams of the anomalous grids ( Figure 8 It was found that the annual demand and number of sites under these regional grids were very low. In order to save on operating and management costs, these sites with very low demand could be cancelled.

[0115] The three diurnal patterns are designated as Category 0, Category 1, and Category 2. Their visualized feature radar plots are then observed. Figure 9 The study found that Category 0 had higher riding distance and duration than the other two categories, while Category 2 had significantly higher number of rides. Category 1, except for the number of rides, was similar to Category 2 in all three metrics. Therefore, Category 0 resembles a low-demand pattern on an annual scale, Category 1 resembles a medium-demand pattern, and Category 2 resembles a high-demand pattern. Combining this with other charts in the visualization system can provide users with a more detailed analysis of the evolution of shared bicycle demand patterns and its influencing factors.

[0116] Specifically, the data cleaning process in S1 involves filtering data that meets one of the following three conditions:

[0117] 1) Cycling time less than 1 minute;

[0118] 2) The cycling distance is less than 100 meters;

[0119] 3) The starting site contains abnormal sites. Sites with fewer than 50 records per year are defined as abnormal sites.

[0120] Furthermore, the calculation process in S2 of dividing the area of ​​the target region into multiple target sub-grid regions is as follows:

[0121]

[0122] Where, N col and N row 111.32 represents the number of columns and rows that can be obtained, where 111.32 indicates that each degree of latitude and longitude is approximately 111.32 kilometers. lon and lat represent latitude and longitude, respectively. Calculations show that, based on the target sub-grid area of ​​1km × 1km, the area covered by the Citibank shared bicycle dataset can be divided into 29 rows and 16 columns. After removing grids that do not meet the required number, the effective number of grids is 249.

[0123] Specifically, the calculation of the S6 weather index includes the following steps:

[0124] Step 6-1: Calculate the center coordinates of the target subgrid region, and then use the semi-versus formula to calculate the distance d from each weather station to the center of the target subgrid region:

[0125]

[0126] Among them [center] lat ,center lon ] represents the center coordinates of the target subgrid region, d represents the spherical distance from the meteorological station to the center coordinates of the target subgrid region, r represents the radius of the sphere, lon and lat represent the latitude and longitude of the meteorological station, respectively. min and lat max These represent the minimum and maximum latitude values ​​of the target subgrid region, respectively. min and lon max These represent the minimum and maximum longitude values ​​of the target sub-grid region, respectively. Next, the inverse distance weights for each meteorological station are normalized. For M meteorological stations, the normalized inverse distance weight α for the i-th meteorological station... i The weighting value for weather stations is calculated as follows:

[0127]

[0128] Where p is the distance attenuation parameter (usually p = 2);

[0129] Step 6-2: Calculate the weighted meteorological index S for each meteorological indicator (including average temperature, average wind speed, and precipitation) in the target sub-grid area. agg Then perform normalization:

[0130]

[0131] Among them, S i Let S represent the meteorological index of the i-th weather station. min With S max Let represent the minimum and maximum values ​​of the weighted meteorological indicators, respectively, and V be the normalized result. Next, the Pearson correlation coefficient is used to measure the correlation between each meteorological indicator and the demand for shared bicycles; the calculation process is as follows:

[0132]

[0133] Among them, S i and Y i These represent meteorological indicators and the demand for shared bicycles, respectively. and These are their mean values. After calculation, the absolute values ​​of the correlation coefficients between the three meteorological indicators—temperature, wind speed, and precipitation—and the demand are normalized and used as the weight values ​​for these three meteorological indicators.

[0134] Step 6-3: Weather events are identified using binary numbers, with each bit representing an event: 1 indicates occurrence, and 0 indicates non-occurrence. First, the Pearson correlation coefficient is used to measure the correlation between different weather events and demand. This coefficient is then normalized and used as the weight for each weather event. Next, a score is calculated for each weather event. event :

[0135]

[0136] Among them, F ie Let β represent weather event e at the i-th weather station. e This indicates the weighting percentage of weather events.

[0137] Step 6-4: Considering that the raw data of weather events uses binary representation and lacks quantitative statistics, the weather index weighting for meteorological indicators is set to 0.9, and the weather index weighting for weather events is set to 0.1. Finally, the weather index calculation process is as follows:

[0138]

[0139] Where K represents the number of meteorological indicators, λ j and V jLet I represent the weight percentage and normalized value of the j-th meteorological indicator obtained in step 6-2, respectively. weather The value is between 0 and 1, when I weather When the value is ≥0.7, the weather condition is considered excellent; when 0.2≤I weather When I < 0.7, the weather conditions are considered good; when I weather When the value is less than 0.2, the weather conditions are considered bad.

[0140] Specifically, obtaining the shared bicycle data feature representation in S7 includes the following specific steps:

[0141] S7-1 constructs training data from both annual and daily scales. For the annual scale, the OD table obtained in S5 is used to calculate four features for each target sub-grid region: daily demand, average riding time, average riding distance, and peak-hour demand percentage. The feature dimension of each annual time series data sample is 365×4. All annual time series data samples are randomly divided into a 70% annual training set and a 30% annual validation set. For the daily scale, the OD table obtained in S5 is used to calculate five features for each target sub-grid region: hourly demand, average riding time, average riding distance, whether it is a holiday, and the weather index obtained in S6. The feature dimension of each daily time series data sample is 24×5. All daily time series data samples are randomly divided into a 70% daily training set and a 30% daily validation set.

[0142] S7-2 TimeCA model construction. Figure 2 The overall architecture of TimeCA is presented, consisting of three parts: Part 1: data augmentation that enhances time series data samples through random pruning and Gaussian perturbation; Part 2: an autoencoder that encodes the augmented time series data; and Part 3: calculation of the joint loss function based on the output of the autoencoder.

[0143] Part 1 Data Augmentation: First, for each time series data sample X = {x1, x2, ..., x...} of sequence length L... L Random cropping is performed. In the annual time series data sample, X = {x1, x2, ... x} i ,...,x 365}, where x i Represents the characteristics of day i; in a daily-scale time series data sample, X = X′ = {x′1, x′2, ..., x′} j ,...,x′ 24},x′ j This represents the characteristics of the j-th hour of the day;

[0144] Set the length of the overlapping sequence of the cropped samples to l, l∈[2]t_unit [,L], where t_unit is a hyperparameter. Randomly select the starting position t2 for pruning, where t2∈[0,Ll], then t3=t2+l. Next, randomly select two time points t1∈[0,t2-1] and t4∈[t3+1,T-1] to obtain two subsequences x with repeating segments. left and x right These two subsequences generated from the same sample are defined as positive sample pairs.

[0145] The above method was used to prune each annual and daily time series data sample, resulting in two subsequences x with repeating segments for each sample. left and x right [Here, two subsequences are generated for each sample, effectively doubling the total sample size. Therefore, 2N appears in the joint loss function formula later. Annual and daily data are processed separately, but the processing methods are the same. The subsequent autoencoder and joint loss function are also calculated separately for annual and daily data. This can be understood as the model building process itself, performing this data processing procedure on both the input annual and daily data.]

[0146] Then, independent Gaussian perturbations are applied to different feature channels in the subsequence. For each feature channel k and timestamp t of the sequence, the specific perturbation process is as follows:

[0147]

[0148] in, It is a Gaussian perturbation, σ 2 This represents the variance of the noise. After the above perturbation, each sample from the yearly or daily scale will result in two enhanced subsequences x′. left and x′ right These enhanced subsequences serve as input to subsequent autoencoders;

[0149] Partial autoencoder: consists of an encoder and a decoder. The encoder extracts key features from augmented time-series data samples, and the decoder restores and reconstructs the original time-series data samples based on the encoder output;

[0150] The encoder consists of n stacked residual modules. In this invention, n=5. Each residual module contains two dilated convolutional layers: the first layer captures dependencies over a shorter time span, and the second layer expands the receptive field to capture long-term dependencies in the time series. After each convolutional operation, the residual modules are normalized, and the output of each convolutional layer is non-linearly mapped using a non-linear activation function. Then, the input and output are summed through residual connections to accelerate the convergence process and obtain the feature representation. The i-th residual module, ResBlock... i The information transmission process is represented as follows:

[0151]

[0152] Where f represents the normalization and nonlinear activation functions in the residual module, and g represents the convolution operation. and Let these represent the weight matrices of the first and second convolutional layers in the i-th residual module, respectively. and is the dilation coefficient of the two dilated convolution layers in the i-th residual module. The encoder output r, obtained through n stacked residual modules, is used as the input to the decoder.

[0153] The decoder consists of n stacked convolutional modules. In this invention, n=5. Each convolutional module contains a dilated convolutional layer and a one-dimensional convolutional layer. First, the dilated convolutional layer matches the receptive field of the encoder to recover patterns over a longer time range in the time series. Then, the one-dimensional convolutional layer makes the reconstructed time series closer to the original sequence data in terms of local details. The decoder output x is obtained after passing through the n convolutional modules. Recon ;

[0154] Partial three-way joint loss function: This includes soft contrast loss and reconstruction loss, which jointly constrain the model training process. The soft contrast loss is calculated separately at the instance layer and the temporal layer.

[0155] At the instance level, for each pair of original time series data samples Q i Q j After min-max normalization, the distance is calculated, and the soft-assigned weights w are defined using the sigmoid function. I (i,j):

[0156] w I (i,j)=2α·σ(-τ I ·D(Q i Q j ))#

[0157] Where σ is the sigmoid function, which converts the distance into a weight between 0 and 1, and τ Iα is a hyperparameter controlling the sharpness of soft assignment. D is a dynamic time programming algorithm used to measure the similarity between two time series data samples. α controls the weights of these similar but different time series data sample pairs. Next, the similarity p between different time series data sample pairs is calculated. I :

[0158]

[0159] Where, r i,t and r j,t Q represents i Q j The characteristic representation of the encoder output timestamp t.

[0160] The loss of the comparison between the i-th time series sample and other samples at time stamp t in the instance layer can be expressed as:

[0161]

[0162] Soft instance contrast loss The first term is the loss for positive sample pairs, where i+N represents the positive sample corresponding to the i-th time series sample, and the second term is the loss for other time series sample pairs;

[0163] Perform the above-described operations at the instance level on annual-scale time series data samples and daily-scale time series data samples respectively;

[0164] At the time level, a soft loss strategy is used to compare the characteristics of different timestamps within a time series. Adjacent timestamps in time series data are typically similar; the soft assignment, w, is calculated based on the differences between timestamps. T The calculation process is shown below:

[0165] w T (t,t′)=2·σ(-τ T ·|tt′|)#

[0166] Where |tt′| represents the difference between two timestamps t and t′, τ T It is a hyperparameter that controls the soft-assignment sharpness. According to w T (t,t′), the contrast loss at the time level can be expressed as:

[0167]

[0168] Where, r i,t and r i,t′ p is a feature representation of the same time series at different timestamps t and t′. T (i,(t,t′)) represents the similarity of the same time series at different timestamps. Soft temporal contrast loss. The first term represents the loss of the positive sample pair, T represents the length of the time dimension of a single time series sample (number of timestamps), t+T represents the timestamp of the positive sample corresponding to time t, and the second term is the loss of other sample pairs;

[0169] The reconstruction loss is calculated by measuring the Euclidean distance between the reconstructed sequence output by the decoder and the corresponding original sequence:

[0170] l Recon =MSE(Q,x Recon )#

[0171] Where MSE is the mean squared error function, Q is the original time series sample, and x Recon This is the corresponding reconstructed sequence output by the decoder module. Finally, the joint loss function is expressed as:

[0172]

[0173] Where λ is the balance coefficient of the reconstruction loss function, 2N represents the sample size after data augmentation, 2T represents the total number of candidate timestamps for a sample to participate in the time layer comparison after data augmentation, and the loss calculated during the training phase is called the training loss.

[0174] Perform the above-described time-level operations on annual-scale and daily-scale time-series data samples respectively;

[0175] Through continuous iterative training, the quality of the time series feature representations learned by the autoencoder can be improved. In subsequent analysis, the encoder structure of the autoencoder module in the trained model is directly used as a feature extractor to obtain the feature representations of the shared bicycle time series data.

[0176] Training the S7-3 TimeCA model. During training, the balance coefficient λ of the reconstruction loss function was set to 1, the Adam optimizer was used to complete the optimization process, the learning rate was set to 0.001, the batch size was set to 8, and the maximum number of optimization iterations was set to 200.

[0177] A balance coefficient λ is set, and the Adam optimizer is used to complete the optimization process. When the reconstruction loss function value no longer changes, the acclimatized TimeCA model is obtained. After training, data augmentation is performed on the validation set samples in the same manner as during the training phase. The encoder extracts features and performs joint loss calculation to obtain the validation loss. The relative error between the validation training loss and the validation loss should not exceed a preset threshold δ (δ = 15% in this invention).

[0178] S7-4 uses a pre-trained TimeCA model for feature extraction. After training, the encoder module is used for feature extraction of the original sequence, outputting a feature representation vector. The original annual-scale time series data outputs a 320-dimensional feature representation, i.e., the embedding. year The shape is (N×320); the original diurnal time series data output diurnal feature representation embedding. day The shape is (N×20).

[0179] Specifically, the S9 visualization display includes five parts: an overview analysis module, a feature analysis module, a pattern analysis module, a comparative analysis module, and a console.

[0180] Overview Analysis Module: Includes an overview map of regional POI information displayed by aggregating S4 data into 5 major categories. Figure 1 -B1) and a spatiotemporal overview diagram of the requirements for displaying clustering results at the S8-year scale ( Figure 1 -B2).

[0181] POI information reflects the distribution of various important functional areas in a city, and different categories of POIs directly affect the demand for shared bicycles. The regional POI information overview map uses a tabular heatmap (each table corresponds one-to-one with the target sub-grid area) to display the specific categories of POIs in each area (5 major categories aggregated by S4), and uses color coding for POI categories. Users can click on a table to view a bar chart showing all POI categories and their descending order of quantity for the corresponding target sub-grid area.

[0182] The spatiotemporal overview map of demand information also adopts a tabular heatmap, overlaid on the target area map. Users can click to select a target sub-grid area, and each sub-grid area displays demand quantity, station quantity information, and macro-pattern category through color coding, helping users quickly identify high-demand areas, densely populated station areas, or specific macro-patterns (clustering results at an S8-year scale). Furthermore, to reveal the dynamic changes in shared bicycle demand over time, this invention utilizes the multi-dimensional information integration advantages of a time-axis bar chart, correlating daily demand changes (length coding) with daily weather indices (color coding), and forming a combined view with the tabular heatmap to display demand under different weather conditions.

[0183] Feature Analysis Module: Used to display line graphs showing demand changes in target sub-grid regions where demand fluctuates significantly at different time granularities. Figure 1 -C1) and a region-site feature analysis plot used to display demand variations between different subgrid areas ( Figure 1 -C2).

[0184] The demand for shared bicycles exhibits significant temporal fluctuations. The demand for a target sub-grid area varies considerably across different time granularities. For example, during weekdays, shared bicycle demand shows a tidal distribution, with higher demand during morning and evening peak hours. The regional demand change line chart allows users to select a date, change the time granularity to a daily scale, and drill down to view the 24-hour demand changes for the selected sub-grid area. Analyzing demand changes at different time granularities helps users understand more granular patterns of demand evolution.

[0185] Demand levels between different subgrid areas reflect travel connections and demand intensity between different areas of the city, helping users identify areas with travel interactions, locate major cycling routes and hotspots, and provide a basis for vehicle scheduling and urban planning. The area-site feature analysis map provides demand, speed, distance, user type, and vehicle usage type for the selected area's associated external areas and internal sites. The statistical distribution values ​​of speed and distance for commuting between areas reflect cycling efficiency and help users evaluate their travel experience.

[0186] Pattern Analysis Module: Contains a regional macroscopic pattern clustering diagram for dimensionality reduction and visualization of clustering results at the S8-year scale. Figure 1 -D1), a pattern feature radar chart used to display statistical characteristics of different categories in multiple dimensions. Figure 1 -D2) and a regional pattern evolution diagram used to display different demand patterns of the same regional grid on different dates ( Figure 1 -D3).

[0187] To analyze the annual demand patterns of each sub-grid region, a scatter plot is used to visualize the results of S8 clustering and dimensionality reduction, providing an intuitive overview of the clustering results. Each scatter point in the view represents a data sample (i.e., a target sub-grid region), and different colors represent different clustering categories. The interactive design supports hovering over each point to display its corresponding region ID and pattern category, and also supports mouse wheel zooming. Users can use the lasso tool to select data samples to generate a regional pattern evolution diagram for the selected region, enabling evolution pattern analysis.

[0188] The pattern feature radar chart displays the statistical features of different categories in macro-scale patterns (annual scale) and evolution patterns (daily scale) in multiple dimensions, including the number of rides, riding distance, peak period proportion, and riding duration. These four features can effectively cover multiple aspects such as spatial scale and travel behavior, helping users identify the differences in statistical features of different patterns.

[0189] The same grid area may exhibit different demand patterns on different dates. This temporal evolution characteristic reflects the complex evolution of user demand for shared bicycles over time and with external factors. The pattern evolution diagram uses different colors to encode different evolution patterns and employs a horizontal axis layout to display the daily pattern evolution information of sub-grid areas at a daily scale, helping users analyze the evolution patterns of shared bicycle demand in each sub-grid area.

[0190] Comparative Analysis Module: Includes regional demand comparison charts that display demand patterns for different regions, the same region's grid on different dates, and different regions' grids. Figure 1 -E).

[0191] The demand patterns for shared bicycles vary across different regions, both temporally and spatially. The same regional grid may exhibit different micro-demand patterns on different days; different regional grids may also exhibit different macro-demand patterns. In-depth analysis of these differences helps users understand the reasons behind these patterns and provides data support for developing operational strategies and optimizing resource allocation. In the regional demand comparison chart, each region is represented by an information summary card, using a three-column vertical layout. The top column displays the regional grid number and time range; the middle column displays multiple characteristic indicators, including average daily ridership, number of stations, and average riding time; the bottom column displays the region's POI aggregation categories (S4 aggregation results) and its top five detailed POI categories. This view integrates these key characteristic indicators from each sub-grid region, presenting a core information summary related to shared bicycles for each region, helping users conduct horizontal comparisons and correlation analyses of demand across different regions.

[0192] Console: Contains a console view ( Figure 1 -A). This can be used to adjust the time range of all other modules, and can add information summary cards for specified regions and time ranges (displayed in the regional demand comparison chart in the comparative analysis module).

[0193] The TimeCA model proposed in this invention can effectively capture the differences between various demand patterns. Furthermore, by combining weather and POI information, it can help users better understand the underlying factors influencing shared bicycle demand patterns. For example, ... Figure 4 As shown, the weather index was poor at the beginning and end of the year, making it difficult for high-demand areas to maintain high-frequency shared bicycle demand, resulting in frequent changes in demand patterns. In contrast, the POI categories in medium-demand areas were mostly transportation areas, so their evolution patterns were relatively stable.

[0194] To verify the effectiveness of the proposed TimeCA model, a time series classification task was performed using the UEA dataset. The selected baseline methods included the traditional time series feature representation method PAA and other self-supervised methods TimesURL and TS2Vec. The classification accuracy of each model is shown below. Figure 5 As shown, the results indicate that TimeCA achieves high, or even the highest, classification task accuracy on most datasets.

[0195] Meanwhile, to verify the effectiveness of the three modules—data augmentation, autoencoder, and soft contrastive learning loss—ablation experiments were conducted. Four datasets were selected from the UEA dataset as the benchmark datasets for the ablation experiments, including ERing (ER), FingerMovements (FM), Heartbeat (HB), and SelfRegulationSCP1 (SCP1). These four datasets differ significantly in feature dimensions, covering low-dimensional to high-dimensional cases. They also cover sequence lengths from long to short and include both binary and multi-class classification tasks, making them representative to a certain extent.

[0196] [Ablation Experiment Methodology] The ablation experiments conducted in this invention involved four variations of TimeCA: removing data augmentation from TimeCA, removing residual connections and dilated convolutions from the encoder, removing the decoder, and removing the soft contrastive learning loss strategy. After removing data augmentation, two identical complete time series were directly used as inputs to TimeCA. Removing residual connections and dilated convolutions from the encoder involved removing residual connections and replacing dilated convolutions with one-dimensional convolutions. The decoder removal experiment did not utilize reconstruction loss to constrain model training. Removing the soft contrastive learning loss strategy involved adjusting the hyperparameter τ... I and τ T All values ​​are set to 0, employing a hard contrastive learning loss strategy. Experimental results are as follows: Figure 6 As shown, each module can improve model performance on different types of datasets.

[0197] To further verify the effectiveness of this invention, a Likert scale was used to construct a questionnaire to evaluate user experience and collect feedback. The survey recruited 13 volunteers, including 4 undergraduate students with introductory experience in visualization and visual analytics, 6 researchers in visualization and visual analytics, and 3 researchers in the transportation field. Based on the visual analytics task, the following eight questions were developed:

[0198] Q1: In the overview analysis, I grasped the spatiotemporal distribution characteristics of shared bicycle demand.

[0199] Q2: In pattern analysis, I understood the meaning and causes of different categories of macro patterns.

[0200] Q3: In the pattern analysis, I learned about the daily pattern evolution patterns and causes in different macro-pattern regions.

[0201] Q4: I discovered points of interest and gained valuable data insights in the overview or pattern analysis.

[0202] Q5: During the analysis, I learned about the impact of external information such as weather and POI on demand patterns.

[0203] Q6: The visual analytics system is simple to operate, user-friendly, and easy to use.

[0204] Q7: The view is clear, the color coding is reasonable, and the readability is high.

[0205] Q8: The visual analytics workflow is logical and clear, and can guide me to complete data exploration and analysis.

[0206] Before the formal experiment began, each volunteer was given 20 minutes to be introduced to the spatiotemporal scope of the research dataset, the meaning of the feature indicators, and the basic usage and analysis process of the visual analytics system. Afterwards, volunteers had 30 minutes for free exploration. Upon completion, they were required to rate the questions raised, with ratings from highest to lowest as follows: 5 points (strongly agree), 4 points (agree), 3 points (uncertain), 2 points (disagree), and 1 point (strongly disagree). The results are as follows: Figure 7 As shown, most respondents gave the system a rating of 4 points or higher (agree), indicating that volunteers generally recognized the system's visual and process design. In the interviews, they agreed that the visual analytics system effectively helped them discover the spatiotemporal distribution characteristics of shared bicycle demand, the macro-level patterns of regional demand, and the evolutionary patterns of daily demand, and enabled detailed analysis of points of interest.

[0207] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit it. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can be made to the technical solutions of the present invention without departing from the spirit and scope of the technical solutions of the present invention, and all such modifications or substitutions should be covered within the scope of the claims of the present invention.

Claims

1. A visual analysis method for shared bicycle demand patterns, characterized in that: Includes the following steps: S1 acquires shared bicycle data and performs data cleaning; S2 divides the area of ​​the target region into multiple target sub-grid regions; S3 identifies specific stations in the cleaned shared bicycle dataset: The DBSCAN clustering algorithm is used to identify the accurate location of each station in the shared bicycle dataset. First, the list of latitude and longitude coordinates of each station is traversed and counted. Then, the DBSCAN algorithm is used to cluster the latitude and longitude coordinates of each station. The neighborhood radius is set to Q meters, and the minimum number of samples in the neighborhood is half the number of points in the list. Finally, the center coordinates of the largest cluster are used as the unique latitude and longitude coordinates of the station. S4 aggregates the cleaned shared bicycle dataset into five categories based on their functions: life services, public services, housing and land, transportation, and leisure and tourism. S5 site and grid information data aggregation: The cleaned shared bicycle data is stored as Parquet files by month according to the time dimension from two perspectives: the site layer and the target sub-grid area layer. Output two new OD tables. One OD table records the order information from the starting station to the destination station for each hour in detail. The other OD table records the order information from the starting target subgrid area layer to the destination target subgrid area layer for each hour in detail. S6 uses a spatial weighted average method, taking the geographic center of each target subgrid area as the reference point, constructing an inverse distance weight based on the spatial distance between each meteorological station and the reference point, and weighting temperature and wind speed to calculate the quantitative value of weather conditions. S7 Shared Bike Data Feature Representation: Shared bike data from different target sub-grid areas are represented using two time granularities: annual and daily. In the annual analysis, the sample size N is the effective target sub-grid area obtained in step 2, and the time step T is 365. The four features selected as feature dimensions F are daily demand, average riding time, average riding distance, and peak hour demand percentage. After standardizing the data, we obtain the model input X∈R. N×T×F The shape is N×365×4, and the feature extraction process is represented as Embedding = Model(X). The feature representation of each data sample represents the 320-dimensional feature representation of that sample in the embedding space, that is, the shape of the Embedding is N×320, which serves as the annual scale feature representation. year ; The daily scale represents the demand evolution pattern of the same target sub-grid area on different dates, where N' represents the number of different dates (365), T' is 24 hours, and the feature dimension F' includes five features: hourly demand, average cycling time and average cycling distance, whether it is a holiday, and the weather index of the day. The resulting model input X∈R N′×T′×F′ After the same feature extraction process as on an annual scale, the daily scale feature representation (Embedding) is obtained. day ; S8 Clustering and Dimensionality Reduction: This involves embedding the feature representation obtained in step S7. year and Embedding day K-means clustering algorithm was used for cluster analysis, and the elbow method was used to determine the optimal number of clusters to identify the macroscopic and microscopic demand patterns of different target subgrid regions. Clustering was performed at the annual scale (5 clusters) and at the daily scale (3 clusters). Principal component analysis was then used to analyze the feature representation obtained in step 8. year and Embedding day Dimension reduction is performed to obtain its two-dimensional projection; S9 visualizes all the results from S1 to S8; At the annual scale, the clusters are divided into five categories: high-demand mode, medium-demand mode, low-demand mode, abnormal mode 1, and abnormal mode 2. The visualized clustering diagram shows that the high-demand mode has a higher number of rides and a higher proportion of peak hours compared to the other two modes, while its riding distance and duration are similar to the medium-demand mode. The low-demand mode has a similar number of rides to the medium-demand mode, but its riding distance and duration are significantly higher than the other two modes. Abnormal mode 1 and abnormal mode 2 are two outlier clusters belonging to abnormal target sub-grid regions. Observing their annual demand changes and information summary diagrams reveals that the annual demand and number of stations in these abnormal target sub-grid regions are both less than the set threshold; therefore, these stations are removed. The three daily-scale patterns were named Category 0, Category 1, and Category 2, respectively. By observing their visualized feature radar charts, it was found that Category 0 had higher cycling distance and cycling time than the other two categories, Category 2 had higher number of rides than the other two categories, while Category 1 had similar values ​​to Category 2 in all three indicators except for the number of rides.

2. The visual analysis method for shared bicycle demand patterns as described in claim 1, characterized in that: The data cleaning process in S1 involves filtering data that meets one of the following three conditions: 1) Cycling time less than 1 minute; 2) The cycling distance is less than 100 meters; 3) The starting site contains abnormal sites. Sites with fewer than 50 records per year are defined as abnormal sites.

3. The visual analysis method for shared bicycle demand patterns as described in claim 2, characterized in that: The calculation process in S2, which divides the area of ​​the target region into multiple target sub-grid regions, is as follows: Where, N col and N row These represent the number of columns and rows that can be obtained, respectively. 111.32 means that each degree of latitude and longitude is approximately 111.32 kilometers. lon and lat represent latitude and longitude, respectively.

4. The visual analysis method for shared bicycle demand patterns as described in claim 3, characterized in that: The calculation of the S6 weather index includes the following specific steps: Step 6-1: Calculate the center coordinates of the target subgrid region, and then use the semi-versus formula to calculate the distance d from each weather station to the center of the target subgrid region: Among them [center] lat ,center lon ] represents the center coordinates of the target subgrid region, d represents the spherical distance from the meteorological station to the center coordinates of the target subgrid region, r represents the radius of the sphere, lon and lat represent the latitude and longitude of the meteorological station, respectively. min and lat max These represent the minimum and maximum latitude values ​​of the target subgrid region, respectively. min and lon max Let represent the minimum and maximum longitude values ​​of the target subgrid region, respectively; then, the inverse distance weights for each meteorological station are normalized; for M meteorological stations, the normalized inverse distance weight α of the i-th meteorological station is... i The weighting value for weather stations is calculated as follows: Where p is the distance attenuation parameter; Step 6-2: Calculate the weighted meteorological index S for each meteorological index in the target sub-grid area. agg Then perform normalization: Among them, S i Let S represent the meteorological index of the i-th weather station. min With S max Let represent the minimum and maximum values ​​of the weighted meteorological indicators, respectively, and V be the normalized result. Next, the Pearson correlation coefficient is used to measure the correlation between each meteorological indicator and the demand for shared bicycles; the calculation process is as follows: Among them, S i and Y i These represent meteorological indicators and the demand for shared bicycles, respectively. and These are their mean values; after calculation, the absolute values ​​of the correlation coefficients between the three meteorological indicators—temperature, wind speed, and precipitation—and the demand are normalized and used as the weight values ​​for these three meteorological indicators. Step 6-3: Weather events are identified using binary numbers, with each bit representing an event (1 indicating occurrence, 0 indicating non-occurrence). First, the Pearson correlation coefficient is used to measure the correlation between different weather events and demand. This coefficient is then normalized and used as the weight of each weather event. Next, the score for each weather event is calculated. event : Among them, F ie Let β represent weather event e at the i-th weather station. e This indicates the weighting percentage of weather events. Step 6-4: Set the weather index weighting of meteorological indicators to 0.9 and the weather index weighting of weather events to 0.1; finally, the weather index calculation process is as follows: Where K represents the number of meteorological indicators, λ j and V j These represent the weight percentage and normalized value of the j-th meteorological indicator obtained in step 6-2, respectively; I weather The value is between 0 and 1, when I weather When the value is ≥0.7, the weather condition is considered excellent; when 0.2≤I weather When I < 0.7, the weather conditions are considered good; when I weather When the value is less than 0.2, the weather conditions are considered bad.

5. A visual analysis method for shared bicycle demand patterns as described in claim 4, characterized in that: The specific steps for obtaining the shared bicycle data feature representation in S7 are as follows: S7-1 constructs training data from annual and daily scales respectively. For the annual scale, the OD table obtained from S5 is used to calculate four features for each target subgrid area: daily demand, average riding time, average riding distance, and peak hour demand ratio. The feature dimension of each annual scale time series data sample is 365×4. All annual scale time series data samples are randomly divided into 70% annual scale training set and 30% annual scale validation set. For the daily scale, the OD table obtained from S5 is used to calculate five features for each target subgrid area: hourly demand, average riding time and average riding distance, whether it is a holiday, and the weather index obtained from S6. The feature dimension of each daily scale time series data sample is 24×5. All daily scale time series data samples are randomly divided into 70% daily scale training set and 30% daily scale validation set. The S7-2 TimeCA model is constructed in three parts: Part 1: data augmentation by randomly pruning and Gaussian perturbation to enhance time series data samples; Part 2: autoencoder to encode the augmented time series data; Part 3: calculation of joint loss function based on the output of the autoencoder. Part 1 Data Augmentation: First, for each time series data sample X = {x1, x2, ..., x...} of sequence length L... L } Perform random cropping; in the annual time series data sample, X = {x1, x2, ... x} i ,...,x 365 }, where x i Represents the characteristics of day i; in In a daily-scale time series data sample, X = X′ = {x′1, x′2, ..., x′} j ,...,x′ 24 },x′ j This represents the characteristics of the j-th hour of the day; Set the length of the overlapping sequence of the cropped samples to l, l∈[2] t_unit ,L], t_unit is a hyperparameter; randomly select the starting position t2 for pruning, and t2∈[0,Ll], then t3=t2+l; then randomly select two time points t1∈[0,t2-1] and t4∈[t3+1,T-1], to obtain two subsequences x with repeated segments. left and x right These two subsequences generated from the same sample are defined as positive sample pairs; The above method was used to prune each annual and daily time series data sample, resulting in two subsequences x with repeating segments for each sample. left and x right ; Then, independent Gaussian perturbations are applied to different feature channels in the subsequence. For each feature channel k and timestamp t of the sequence, the specific perturbation process is as follows: in, It is a Gaussian perturbation, σ 2 This represents the variance of the noise; after the above perturbation, each sample from the yearly or daily scale will result in two enhanced subsequences x′. left and x′ right These enhanced subsequences serve as input to subsequent autoencoders; Partial autoencoders consist of an encoder and a decoder; the encoder extracts key features from augmented time-series data samples, and the decoder restores and reconstructs the original time-series data samples based on the encoder output. The encoder consists of n stacked residual modules, each containing two dilated convolutional layers: the first layer captures dependencies over a shorter time span, and the second layer expands the receptive field to capture long-term dependencies in the time series. After each convolutional operation, the residual modules are normalized, and the output of each convolutional layer is non-linearly mapped using a non-linear activation function. Then, the input and output are summed through residual connections to accelerate the convergence process and obtain the feature representation. The i-th residual module, ResBlock... i The information transmission process is represented as follows: Where f represents the normalization and nonlinear activation functions in the residual module, and g represents the convolution operation. and Let these represent the weight matrices of the first and second convolutional layers in the i-th residual module, respectively. and is the dilation coefficient of the two dilated convolutions in the i-th residual module; the encoder output r is obtained through n stacked residual modules and used as the input of the decoder; The decoder consists of n stacked convolutional modules. Each module contains a dilated convolutional layer and a one-dimensional convolutional layer. First, the dilated convolutional layer matches the encoder's receptive field to recover patterns over a longer time range in the time series. Then, the one-dimensional convolutional layer makes the reconstructed time series more closely resemble the original sequence data in terms of local details. The decoder output x, obtained after passing through the n convolutional modules, is then obtained. Recon ; Partial three-way joint loss function: including soft contrast loss and reconstruction loss, which jointly constrain the model training process; soft contrast loss is calculated separately from the instance layer and the temporal layer. At the instance level, for each pair of original time series data samples Q i Q j After min-max normalization, the distance is calculated, and the soft-assigned weights w are defined using the sigmoid function. I (i,j): w I (i,j)=2α·σ(-τ I ·D(Q i ,Q j ))# Where σ is the sigmoid function, which converts the distance into a weight between 0 and 1, and τ I α is a hyperparameter controlling the soft-assignment sharpness; D is dynamic time programming used to measure the similarity between two time series data samples; α controls the weights of these similar but different time series data sample pairs; then, the similarity p of different time series data sample pairs is calculated. I : Where, r i,t and r j,t Q represents i Q j The characteristic representation of the encoder output timestamp t; The loss of the comparison between the i-th time series sample and other samples at time stamp t in the instance layer can be expressed as: Soft instance contrast loss The first term is the loss for positive sample pairs, where i+N represents the positive sample corresponding to the i-th time series sample, and the second term is the loss for other time series sample pairs; Perform the above-described operations at the instance level on annual-scale time series data samples and daily-scale time series data samples respectively; At the time level, soft allocation is calculated using the difference between timestamps, w T The calculation process is shown below: w T (t,t′)=2·σ(-τ T ·|t-t′|)# Where |tt′| represents the difference between two timestamps t and t′, τ T It is a hyperparameter that controls the soft-sharpening sharpness; according to w T (t,t′), the contrast loss at the time level can be expressed as: Where, r i,t and r i,t′ p is a feature representation of the same time series at different timestamps t and t′. T (i,(t,t′)) represents the similarity of the same time series at different timestamps; soft temporal contrast loss The first term represents the loss of the positive sample pair, T represents the time dimension length of a single time series sample, t+T represents the timestamp of the positive sample corresponding to time t, and the second term is the loss of other sample pairs; The reconstruction loss is calculated by measuring the Euclidean distance between the reconstructed sequence output by the decoder and the corresponding original sequence: l Recon =MSE(Q,x Recon )# Where MSE is the mean squared error function, Q is the original time series sample, and x Recon This is the corresponding reconstructed sequence output by the decoder module; finally, the joint loss function is expressed as: Where λ is the balance coefficient of the reconstruction loss function, 2N represents the sample size after data augmentation, 2T represents the total number of timestamp candidates for a sample to participate in the time layer comparison after data augmentation, and the loss calculated during the training phase is called the training loss. Perform the above-described time-level operations on annual-scale and daily-scale time-series data samples respectively; S7-3 TimeCA Model Training: Set the balance coefficient λ, use the Adam optimizer to complete the optimization process, and obtain the acclimatized TimeCA model when the reconstruction loss function value no longer changes; after training, perform data augmentation on the validation set samples in the same way as in the training phase, extract features from the encoder and perform joint loss calculation to obtain the validation loss; verify that the relative error between the training loss and the validation loss is not higher than the preset threshold δ. S7-4 uses a pre-trained TimeCA model for feature extraction. After training, the encoder module is used for feature extraction of the original sequence, outputting a feature representation vector. The original annual-scale time series data outputs a 320-dimensional feature representation, i.e., the embedding. year The shape is N×320; the original diurnal time series data output diurnal feature representation embedding. day The shape is N×20.

6. The visual analysis method for shared bicycle demand patterns as described in claim 5, characterized in that: The S9 visualization display includes five parts: overview analysis module, feature analysis module, pattern analysis module, comparative analysis module, and console. Overview Analysis Module: Includes an overview map of regional POI information that displays the five major categories aggregated in S4, and a spatiotemporal overview map of demand information that displays the clustering results at the S8-year scale. Feature Analysis Module: Used to display line graphs showing significant fluctuations in demand across target subgrid areas at different time granularities, and region-site feature analysis graphs showing demand changes between different subgrid areas; Pattern Analysis Module: Includes a regional macro-pattern clustering map for dimensionality reduction display of clustering results at the S8-year scale, a pattern feature radar chart for multi-dimensional display of statistical characteristics of different categories, and a regional pattern evolution map for displaying different demand patterns of the same regional grid on different dates. Comparative Analysis Module: Contains regional demand comparison charts to display regional demand patterns for different regions, the same regional grid on different dates, and different regional grids.

Images