extfieldsofstudy:
Computer Science
Medicine

text:
Automated Multistep Parameter Identification of SPMSMs in Large-Scale Applications Using Cloud Computing Resources

Parameter identification of permanent magnet synchronous machines (PMSMs) represents a well-established research area. However, parameter estimation of multiple running machines in large-scale applications has not yet been investigated. In this context, a flexible and automated approach is required to minimize complexity, costs, and human interventions without requiring machine information. This paper proposes a novel identification strategy for surface PMSMs (SPMSMs), highly suitable for large-scale systems. A novel multistep approach using measurement data at different operating conditions of the SPMSM is proposed to perform the parameter identification without requiring signal injection, extra sensors, machine information, and human interventions. Thus, the proposed method overcomes numerous issues of the existing parameter identification schemes. An IoT/cloud architecture is designed to implement the proposed multistep procedure and massively perform SPMSM parameter identifications. Finally, hardware-in-the-loop results show the effectiveness of the proposed approach.

To overcome such an issue, many approaches have been proposed.
A simple method to obtain a full-rank identification problem is based on fixing some parameters to their nominal values, as in [11,12]. However, such methodology is compromised by mismatches occurring between nominal and actual values [7,8]. Moreover, those methods are not practical for large-scale systems where some nominal parameters may be unknown. Another approach is based on the use of extra sensors, such as thermal sensors and power or torque meters [13,14]. This approach is not suitable in large-scale applications due to the increased complexity and costs [7].

•
We solve the rank-deficient problem by employing a multistep procedure based on three AdNNs with mutual updating without using signal injection as in [8][9][10][15][16][17][18][19]; • Our method does not require exact knowledge of nominal values, overcoming the issues of [11,12]; • Extra sensors are not required as in [13,14]; • Unlike [20,21], the proposed method also handles zero d-axis current control schemes; • An automated solution based on the R-statistic algorithm allows one to properly identify the SPMSM steady-state operating conditions in which the multistep procedure operates; • Low computational requirements on the edge-side make the proposed method highly suitable for IoT integrations, without any change on the motor control unit; • Off-the-shelf cloud technologies implement a hardware-in-the-loop (HIL) setup while simulation results show the effectiveness of the proposed approach.
Compared to the existing solutions, we point out how the main advantage of the proposed method is its applicability to large-scale industrial processes due to an automated identification scheme and the limitation of human intervention. This paper is organized as follows. In Section 2 the rank-deficiency problem is discussed. The proposed algorithm is detailed in Section 3. Section 4 describes the HIL setup used for validation, while Section 5 presents the simulation results. Finally, concluding remarks are reported in Section 6.

Rank-Deficiency Problem
We consider a field-oriented control (FOC) of an SPMSM drive. The dynamical model of an SPMSM in the dq rotating synchronous reference frame is described by the following equations [3]: where v d , v q are the dq-axes voltages, R s and L s are the stator resistance and dq-axes inductance, respectively, i d , i q are the dq-axes currents, ω r is the electrical rotor speed, and ψ is the rotor flux linkage. Note that if a zero d-axis current control is performed, at the steady state, (1) and (2) can be rewritten in the following form: v q = R s i q + ω r ψ.
From (3), L s can be estimated independently from R s and ψ through the d-axis voltage, speed, and q-axis current measurements. Instead, R s and ψ cannot be estimated simultaneously from (4), since only one equation is available. This is commonly known as rank deficiency problem for the parameter identification of SPMSMs controlled with i d = 0. In this paper, a multistep procedure is proposed to solve such a problem.

Novel Multistep Parameter Identification Algorithm Based on Adaline NNs
The proposed identification scheme is shown in Figure 1. Starting from the measurement data, the RStatIn block recognizes the SPMSM steady-state conditions, the MovAvIn filters the measurement data, the decision-making algorithm (DMA) rules the operations of the three AdNNs, one per parameter: rotor flux linkage (AdNN1), stator resistance (AdNN2), and inductance (AdNN3). RStatOut detects the steady state of the parameter estimations which are filtered by the MovAvOut. These blocks are detailed in the following subsections.

Measurement Data
The measurement data contain the samples of i q , ω r , v d , and v q , which are gathered in signal m in Figure 1. It is assumed that the d-axis current is zero to perform a maximum torque per ampere control and that no additional signals are injected for the purpose of the parameter identification.

RStatIn: Identificator of Steady-State Conditions of the SPMSM
R-statistic is a statistical method developed to automatically distinguish transient states from steady states in noisy processes [27]. The RStatIn block employs the R-statistic method for the automated identification of the steady-state conditions of the SPMSM. This operation is required since the AdNN estimators operate only during the steady-state conditions of the motor, as will be detailed below. The SPMSM operates in steady-state conditions when its speed and electromagnetic torque are constant. Therefore, to identify the steady-state conditions of the SPMSM, the R-statistic is applied on the electrical rotor speed and q-axis current.
conditions when its speed and electromagnetic torque are constant. Therefore, to identify the steady-state conditions of the SPMSM, the R-statistic is applied on the electrical rotor speed and q-axis current. The R-statistic procedure calculates the following indices [27]: where and are the indices for and , respectively, is the k-th sample processed, and is the sample's window length. and are gathered in signal _ in Figure 1.
where , i q are noise signals based on the Box-Muller method, and i q are the standard deviations, while and are independent samples chosen from the uniform distribution in the interval [0, 1]. The noise signals are introduced to avoid numerical issues [27]. and tend to be around 1 as the process tends to the steady state. Instead, during transients, the indices are expected to be greater than 1. Hence, the process is in the steady state if ≤ , where is a critical threshold. In [27], the author suggests manually tuning according to the actual responses. The R-statistic procedure calculates the following indices [27]: where R ω r and R i q are the indices for ω r and i q , respectively, k is the k-th sample processed, and N is the sample's window length. R ω r and R i q are gathered in signal R_in in Figure 1. ω r,n and i q,n are computed as follows: where w ω r , w i q are noise signals based on the Box-Muller method, σ ω r and σ i q are the standard deviations, while r 1 and r 2 are independent samples chosen from the uniform distribution in the interval [0, 1]. The noise signals are introduced to avoid numerical issues [27]. R ω r and R i q tend to be around 1 as the process tends to the steady state. Instead, during transients, the indices are expected to be greater than 1. Hence, the process is in the steady state if R ≤ R crt , where R crt is a critical threshold. In [27], the author suggests manually tuning R crt according to the actual responses.
R crt , N, σ ω r , and σ i q are the tuning parameters of RStatIn found with trial-and-error approach. The tuning of these parameters does not need to be repeated for different electrical drives since they do not depend upon physical parameters of the machine. Note that the R-statistic algorithm acts as a filter on the sample window of length N, hence it is affected by delay in identifying the beginning or the end of a steady state. This issue has been considered and properly compensated through to the MovAvOut block, as will be shown below.

MovAvIn: Moving Average of Measurement Data
The MovAvIn block is used to filter the input signals ω r and i q as follows: where ω r,ave and i q,ave are the moving averages (gathered in signal ma_in in Figure 1). Note that the windows size, N, is the same as in RStatIn. The outputs of MovAvIn are considered as the current steady-state operating conditions of the SPMSM by the DMA.

Decision-Making Algorithm
The DMA enables the AdNNs through the Enable signals to execute the multistep parameter identification. AdNN1 and AdNN2 are enabled in different operating conditions while AdNN3 works simultaneously with AdNN1 or AdNN2. The DMA employs the information provided by the RStatIn block, R_in = [R ω r , R i q ], the current operating conditions from MovAvIn, ma_in = [ω r,ave , i q,ave ], and a convergence condition to decide which of the two estimators (AdNN1 or AdNN2) should be enabled. Moreover, the DMA employs information provided by the RStatOut block, R_out = Rψ * ave , RR * s,ave , RL * s,ave , to accept only the filtered parameter estimations, ma_out = ψ * ave ,R * s,ave ,L * s,ave , which reach the steady state. Finally, the DMA performs a stop criterion which automatically concludes the multistep procedure without human intervention. The DMA outputs are the Enable signals for the three AdNNs and est_star = [ψ * ,R * s ,L * s ]. Moreover, the DMA uses three internal states, i.e.,ψ stop (x),R stop (y), andL stop (z), called partial estimations, to computê ψ * ,R * s ,L * s when the AdNNs are not active; also, they are used in the stop criterion. Three indices, x, y, z are updated asynchronously as explained later. Finally, the DMA employs other internal states, called check variables, which are ω rR (y), ω rψ (x), i qR (y), and i qψ (x), representing respectively the speed and q-axis current at which the latest rotor flux linkage and stator resistance estimations have been performed by the AdNNs. All these variables are initialized as follows: the three Enable signals are set to OFF;ψ * ,R * s ,L * s ,ψ stop ,R stop , L stop , x, y, and z are set to zero; ω rR and i qR are set to M and ε, respectively, where M is a large number and ε is a small number; ω rψ and i qψ are set to zero.
The operation of the DMA is described in the flow chart in Figure 2 At each iteration, the DMA checks if the SPMSM is in the steady-state condition. If the SPMSM is in the steady state, i.e., R ω r (k) ≤ R crt and R i q (k) ≤ R crt , the DMA implements Algorithm A1, providing the first parameter estimation by activating AdNN1 and AdNN3 if x = 0. The first estimation initializes the multistep procedure, obtaining the first partial estimationŝ ψ stop andL stop . Once the first estimationψ stop has been achieved, the SPMSM stationary operating conditions are inspected. As proven in Appendix B, to make sure that stator resistance and rotor flux linkage estimation errors asymptotically converge to zero, the AdNN1 and AdNN2 must be activated separately to satisfy the following convergence condition (A13): where i qψ , ω rψ , i qR , and ω rR are the check variables defined above, which express the operating conditions of the SPMSM in which the AdNN1 and AdNN2 are enabled by the DMA. In order to satisfy (11), the DMA activates the AdNN1 if (ω rR /(ω r,ave (k))· i q,ave (k)/i qR < 0.95 and activates the AdNN2 if (ω r,ave (k)/ω rψ )· i qψ /i q,ave (k) < 0.95. The AdNN3 is

Appendix C
In this appendix, we prove how the learning rate of the AdNNs expressed by (14), (18), and (22) guarantee convergence.
According to (A2), to ensure the convergence of the flux estimations provided by AdNN1, the learning rate should be chosen to make 1 − 2 ( ) < 1. Therefore, the following equation must be satisfied: The learning rate which satisfies (A14) is: If R ω r (k) > R crt or R i q (k) > R crt , the DMA performs the operations reported in Algorithm A2. Firstly, the Enable signals are set to OFF, since the AdNNs should work only during steady-state operations of the SPMSM. Then, if the AdNNs were active and their filtered estimations were at the steady state in the previous step, i.e., k − 1, the DMA updatesψ stop ,R stop , andL stop with the values in ma_out(k − 1) and ω rR , ω rψ , i qR , and i qψ with the values in ma_in(k − 1). After the updating of the partial estimations, the DMA checks the stop criterion, which is described in Algorithm 1. This algorithm operates on the sets of the last N stop samples ofψ stop ,R stop , andL stop where ε stop is an arbitrary small value.
If the AdNNs were not active or the provided estimations were not in the steady state in the previous step,ψ stop ,R stop ,L stop , ω rR , ω rψ , i qR , and i qψ are not updated. This allows one to accept only the stationary estimations provided by the AdNNs, avoiding the introduction of estimation errors due to transients of the AdNNs. Note that, if the AdNNs estimations are accepted, the partial estimations assume the values of the filtered estimations provided by the MovAvOut block. The filtering reduces the errors introduced by the estimations performed during transients. In fact, as stated in Section 3.2, the R-statistic algorithm detects with a delay the loss of the steady state of the SPMSM. During this delay, the AdNNs continue to produce estimations which are affected by unwanted perturbations that must be filtered. Finally, the DMA updates the output variables, est_star, according to Algorithm A3. Note that if the AdNNs are activated, the est_star variables are set to the corresponding AdNNs estimations; otherwise, they are set to the partial estimations.
We remark how the partial estimations are updated asynchronously and only if the corresponding ma_out variables have reached their steady states, as described in Algorithm A2. Moreover, the partial estimations represent the results of the proposed algorithm, since they are the only estimations not affected by errors due to the transient states.

Adaline NNs
This section describes the neural network estimators implemented in this work and depicted in Figure 3. These estimators are driven by the Enable signals provided by the DMA. The AdNN1 receives as inputs the measurements of the q-axis current, i q , the electrical rotor speed, ω r , and the q-axis voltage, v q , and the estimationsψ * andR * s provided by the DMA. The weight of the electrical rotor speed is the estimated rotor flux linkage (ψ), updated according to the following: with η ψ as learning rate.v q is the estimated q-axis voltage, expressed as follows: According to Appendix B, to ensure the convergence of the flux estimation, the learning rate is computed with the following formula: where k ψ is a real constant in the interval [−1, 1]. Note that, as in Algorithm A3, if AdNN1 was active in the previous step, thenψ * (k − 1) coincides withψ(k − 1) in (11) providing the following [9],ψ otherwise,ψ * (k − 1) coincides withψ stop (x), allowing the AdNN1 to restart with the best rotor flux linkage estimation available. This reduces the AdNNs transients since they are where is a real constant in the interval [−1, 1]. Note that, as in Algorithm A3, if AdNN1 was active in the previous step, then * ( − 1) coincides with ( − 1) in (11) providing the following [9], otherwise, * ( − 1) coincides with ( ), allowing the AdNN1 to restart with the best rotor flux linkage estimation available. This reduces the AdNNs transients since they are not reset to their initial values. The estimated rotor flux linkage coincides with the output signal of AdNN1 in Figure 1. The AdNN2 shares the same inputs of the AdNN1. The weight of the q-axis current is the estimated stator resistance ( ), updated according to the following: with as learning rate. is the estimated q-axis voltage, expressed as follows: The AdNN2 shares the same inputs of the AdNN1. The weight of the q-axis current is the estimated stator resistance (R s ), updated according to the following: with η R s as learning rate.v q is the estimated q-axis voltage, expressed as follows: According to Appendix B, to ensure the convergence of the resistance estimation, the learning rate is computed with the following formula: where k R s is a real constant in the interval [−1, 1]. As for the AdNN1,R * s (k − 1) is replaced withR s (k − 1) in (14) if AdNN2 was active in the previous step, as follows [9] The estimated stator resistance coincides with the output signal of AdNN2 in Figure 1. The AdNN3 is used for the stator inductance's estimation. It receives as input the measurements i q , ω r , and v d , and the signalL * s from est_star. The weight of the product of the electrical rotor speed and q-axis current is the estimated stator inductance (L s ), updated according to the following: with η L s as learning rate andv d as the estimated d-axis voltage, expressed as follows: According to Appendix B, to ensure the convergence of the resistance estimation, the learning rate is computed with the following formula: where k L s is a real constant in the interval [−1, 1]. If AdNN3 at the previous step was active, thenL The estimated stator inductanceL s coincides with the output signal of AdNN3 in Figure 1.
The values of k ψ , k R s , and k L s affect the convergence speed of the AdNNs. We recommend choosing these values in the interval [0.8, 1] to mitigate the perturbations of the estimations which occur when the SPMSM changes its operating condition.

MovAvOut: Moving Average of Measurement Data
The MovAvOut block shown in Figure 1 performs a moving average of the est_star variables. The function calculates, using (24)-(26), the following values: whereψ * ave ,R * s,ave andL * s,ave (gathered in signal ma_out in Figure 1) are the averages of the est_star variables. Note that, as shown in Algorithm A3, the est_star variables coincide with the output of the AdNNs when these ones are enabled by the DMA. Therefore, the MovAvOut block performs a filtering of the estimations produced by the AdNNs. This operation is essential since it allows one to mitigate the effect of the delay of RStatIn and to properly perform the R-statistic algorithm on the parameter estimations.

RStatOut: Identificator of Steady-State Conditions for Estimated Parameters
The RStatOut block in Figure 1 operates in the same way as RStatIn but on different datasets. Its inputs are the ma_out variables and the outputs are computed using (27)-(29): In these formulas Rψ * ave , RR * s,ave , andL * s,ave (gathered in signal R_out in Figure 1) are the R-statistic indices andψ * n ,R * s,n, andL * s,n are the noisy estimations. Furthermore, in this case, a properly tuned Box-Muller noise is introduced as explained in Section 3.2. The parameters N and R crt and the noise standard deviations σψ * ave , σR * s,ave , and σL * s,ave are equal to those defined for the RStatIn block. This function block is used to reveal the stationarity of the DMA outputs, which coincide with the outputs of the AdNNs estimations during their operations, as shown in Algorithm A3. As explained above, this operation allows one to accept only the AdNNs stationary estimations, avoiding the introduction of estimation errors due to transients.

Hardware-in-the-Loop (HIL) Setup
To evaluate the performance of the proposed solution, the experimental setup in Figure 4 has been arranged. This experimental setup has been inspired by the one shown in [26]. It consists of three levels:
Internet of things (IoT) devices which elaborate and buffer the collected data in order to optimize the bandwidth towards the cloud; 3.
A cloud application where the data is stored and the multistep parameter identification algorithm is executed.
The Simulink PC acts as HIL simulator of the FOC SPMSM. An SPMSM driven by an FOC controller and fed by a pulse-width-modulation (PWM) three-phase inverter has been modeled. PWM is a widely adopted technique for the regulation of power electronic devices such as AC voltage regulators, inverters, rectifiers, etc. [28][29][30]. The FOC controller performs a speed control with zero d-axis current and is made by the cascade of speed and current PI regulators. To bring the model as close to reality as possible, measurement uncertainties are also added to the phase currents and on the DC-link voltage measurements. These uncertainties have been introduced considering datasheets of commercial transducers, i.e., LEM LA 55-P for the currents transducer and LEM LV 25-1000 for the voltage transducer. Moreover, a resolver has been modeled to provide rotor speed and angular position measurements. The SPMSM models with wyewound stator are configured as Bonfiglioli BMD 400V 65 and 170 surface magnet motors. The main motors and inverter parameters are constant during the simulation and they are summarized in Table 1. The model is used to simulate, at fixed sampling time step, two working cycles of six seconds with variable speed and load torque as shown in the next paragraph.    The signals v dq , i dq, and ω r are logged, buffered, and sent to the IoT device via Modbus TCP protocol. The IoT architecture has been implemented through the AWS IoT Greengrass technology and performs gathering, pre-cleansing, storing, and analysis operations on the data provided by the HIL simulator [31][32][33]. These operations are performed by two local IoT lambda functions [32], developed using the Python language. The first lambda function fetches and cleans the field data. Then, ready-to-send data are stored in a local database, while the second lambda function queries for the newest data and wraps it in an MQTT [33] message payload which is finally published to the cloud using an MQTT client instance. In this way, the IoT device communicates asynchronously and securely with the cloud over the MQTT protocol using a key pair provisioning mechanism [33].
The third level of the setup consists of the cloud application, where messages coming from the IoT device are collected by the AWS IoT Core [34,35]. IoT Core allows one to filter and transform the input messages. An IoT Rule, triggered by an input MQTT message, is implemented to identify, transform, and forward the message payload to Amazon S3. S3 is a hierarchical object storage service in which a source bucket is arranged to receive all the data from the IoT device [36]. The source bucket contains a folder per IoT device, and each folder collects the objects, i.e., the time series generated by the HIL simulator, represented as a red circle in Figure 4. S3 Event Notification is used to trigger the multistep parameter identification every time a new object is saved in the source bucket. The identification algorithm is implemented with a lambda function whose results are recorded in a target bucket, as shown in Figure 4. The tuning parameters of the proposed identification algorithm are reported in Table 2. Table 2. Parameters of the proposed algorithm.

Parameters
Value 5a shows a photo of the experimental setup. The IoT device is connected via ethernet cable to the Simulink PC to perform the HIL simulations. The recorded data are sent to the cloud via IEEE 802.11 connection. Note that, even though we employed a simulation environment for the motor data generation, an industrial device has been used as an IoT edge device. Figure 5b shows the IoT device implemented in a real application within an industrial environment interfaced with a data logger and power meters.

Simulation Results
Two different working cycles have been simulated for the two motors. The speed, torque, and d-axis current profiles of the simulated working cycles are reported in Figure 6, where , * , , , , and * are the measured speed, the reference speed, the electromagnetic torque, the load torque, the measured d-axis current, and the reference d-axis

Simulation Results
Two different working cycles have been simulated for the two motors. The speed, torque, and d-axis current profiles of the simulated working cycles are reported in Figure 6, where ω r , ω r * , T e , T L , i d , and i d * are the measured speed, the reference speed, the electromagnetic torque, the load torque, the measured d-axis current, and the reference d-axis current, respectively. Note that the d-axis current is different from zero only during transients of the SPMSMs. Therefore, it does not affect the performances of the AdNN estimators since they operate only during the steady states of the SPMSMs. We consider two data packets made by 108 s of measurement data recorded from 18 working cycles of the two motors.    Figure 7 shows the results of the R-statistic analysis, identifying the rotor speed steady states of the two motors. In this figure, the measured electrical rotor speed ω r , the noisy speed ω r,n obtained using the Box-Muller noise, the R-statistic values R ω r , and the critical  Figure 7a shows how the R-statistic algorithm detects with delay the end of the steady states. Similarly, in Figure 8, the measured q-axis current i q , the noisy current i q,n , the R-statistic values R i q , and the critical value R crt are shown. The R-statistic identifies three q-axis current steady-state time intervals ( Figure 9 shows the AdNNs activation signals 1, 2, and 3 in the first 12 s for the two motors. It is shown in which of the motors' steady states the rotor flux linkage, stator resistance, and stator inductance's AdNN estimators are working. In the first working cycle of the Bonfiglioli 65, the AdNN1 is activated for the first time at half of the rated speed and current and for the second time at the rated speed and at 10% of the rated current; instead, the AdNN2 is activated only at the 10% of rated speed and at the rated current. Instead, in the first working cycle of the Bonfiglioli 170, only the AdNN1 is activated: the first time at 10% of the rated speed and at the rated current, the second time at the rated speed and current, and the third time at the rated speed and half of the rated current. In both cases, the first activation of AdNN1 corresponds to the first estimation to initialize the multistep algorithm while the other two activations satisfy the other conditions expressed in Algorithm A1. The trend of the enable signals is periodical from 6s to the end of the parameter identification, i.e., the time instant in which the stop criterion is satisfied. In particular, for the Bonfiglioli 170, after 6s, the AdNN1 is activated twice during each working cycle: the first time at the rated speed and current and the second time at the rated speed and half of the rated current. Instead, the AdNN2 is activated only once during each working cycle at the 10% of rated speed and at the rated current. The comparison of the operating conditions in which AdNN1 and AdNN2 are enabled shows how the algorithm is working in agreement with convergence condition in (A13).  Figure 9 shows the AdNNs activation signals Enable1, Enable2, and Enable3 in the first 12 s for the two motors. It is shown in which of the motors' steady states the rotor flux linkage, stator resistance, and stator inductance's AdNN estimators are working. In the first working cycle of the Bonfiglioli 65, the AdNN1 is activated for the first time at half of the rated speed and current and for the second time at the rated speed and at 10% of the rated current; instead, the AdNN2 is activated only at the 10% of rated speed and at the rated current. Instead, in the first working cycle of the Bonfiglioli 170, only the AdNN1 is activated: the first time at 10% of the rated speed and at the rated current, the second time at the rated speed and current, and the third time at the rated speed and half of the rated current. In both cases, the first activation of AdNN1 corresponds to the first estimation to initialize the multistep algorithm while the other two activations satisfy the other conditions expressed in Algorithm A1. The trend of the enable signals is periodical from 6s to the end of the parameter identification, i.e., the time instant in which the stop criterion is satisfied. In particular, for the Bonfiglioli 170, after 6s, the AdNN1 is activated twice during each working cycle: the first time at the rated speed and current and the second time at the rated speed and half of the rated current. Instead, the AdNN2 is activated only once during each working cycle at the 10% of rated speed and at the rated current. The comparison of the operating conditions in which AdNN1 and AdNN2 are enabled shows how the algorithm is working in agreement with convergence condition in (A13).  Figure 10 shows the results of the rotor flux linkage identification for the Bonfiglioli 65. The estimations are progressively improved and the percentage estimation error at the end of the identification process is 0.23%. This figure shows how changes its value only when the AdNN1 is disabled and is equal to * at the end of the AdNN1 operation. In particular, the zoomed part in Figure 10a shows that * is affected by perturbations before the AdNN1 is disabled, i.e., when is updated, while the signal * overcomes this issue thanks to the filtering operation of the moving average. Note that a progressive improvement of the estimations is achieved in the first two working cycles since the condition expressed by (A13) is satisfied. At the third working cycle, an accurate estimation has already been achieved and there are no substantial improvements until the stop criterion is satisfied. Figure 10b shows that R-statistic manages to detect only the first transient state of * , while the other transients are too small to be detected. The comparison between Figure 10a,b shows how all the estimations performed by the AdNN1 have reached the steady state and, thus, have been processed by the DMA to update .  Figure 10 shows the results of the rotor flux linkage identification for the Bonfiglioli 65. The estimations are progressively improved and the percentage estimation error at the end of the identification process is 0.23%. This figure shows howψ stop changes its value only when the AdNN1 is disabled and is equal toψ * ave at the end of the AdNN1 operation. In particular, the zoomed part in Figure 10a shows thatψ * is affected by perturbations before the AdNN1 is disabled, i.e., whenψ stop is updated, while the signal ψ * ave overcomes this issue thanks to the filtering operation of the moving average. Note that a progressive improvement of the estimations is achieved in the first two working cycles since the condition expressed by (A13) is satisfied. At the third working cycle, an accurate estimation has already been achieved and there are no substantial improvements until the stop criterion is satisfied. Figure 10b shows that R-statistic manages to detect only the first transient state ofψ * ave , while the other transients are too small to be detected. The comparison between Figure 10a,b shows how all the estimations performed by the AdNN1 have reached the steady state and, thus, have been processed by the DMA to updateψ stop . Figure 11 shows the results of the stator resistance identification for the Bonfiglioli 65. Moreover, in this case, there is a progressive improvement of the estimations and at the end of the identification process a percentage estimation error of 0.35% is achieved. Similar considerations already performed for the flux linkage estimations can be repeated for the updating ofR stop and the R-statistic analysis. Figure 12 shows how the stator inductance estimation accuracy slightly depends on the SPMSM operating condition in which the AdNN3 is activated. The condition in which the minimum estimation error is achieved is the one at the rated current and 10% of the rated speed while the condition in which the maximum error is achieved is the one at 10% of the rated current and at the rated speed. At the end of the parameter identification process, the relative percentage estimation error is 0.11%. Moreover, in this case, the figure shows the importance of filtering the AdNN estimation with the MovAvOut block in order to attenuate perturbations that affect the estimation accuracy. Figure 12b shows that, in this case, R-statistic managed to detect the transients of the estimations. However, all the estimations performed by the AdNN3 reached the steady state and have been processed by the DMA to updateL stop . The percentage estimation errors associated with the filtered parameter estimationsψ * ave ,R * s,ave, andL * s,ave for the entire simulation time of the Bonfiglioli 65 are shown in Figure 13.  Figure 11 shows the results of the stator resistance identification for the Bonfiglioli 65. Moreover, in this case, there is a progressive improvement of the estimations and at the end of the identification process a percentage estimation error of 0.35% is achieved. Similar considerations already performed for the flux linkage estimations can be repeated for the updating of and the R-statistic analysis. Figure 12 shows how the stator inductance estimation accuracy slightly depends on the SPMSM operating condition in which the AdNN3 is activated. The condition in which the minimum estimation error is achieved is the one at the rated current and 10% of the rated speed while the condition in which the maximum error is achieved is the one at 10% of the rated current and at the rated speed. At the end of the parameter identification process, the relative percentage estimation error is 0.11%. Moreover, in this case, the figure shows the importance of filtering the AdNN estimation with the MovAvOut block in order to attenuate perturbations that affect the estimation accuracy. Figure 12b shows that, in this case, R-statistic managed to detect the transients of the estimations. However, all the estimations performed by the AdNN3 reached the steady state and have been processed by the DMA to update . The percentage estimation errors associated with the filtered parameter estimations * , , , * and , * for the entire simulation time of the Bonfiglioli 65 are shown in Figure 13.        Figure 14 shows the parameters identification for the Bonfiglioli 170. The estimations of rotor flux linkage, stator resistance, and inductance have trends similar to those of Bonfiglioli 65. In particular, the rotor flux linkage and stator resistance estimation progressively improve and the percentage estimation errors at the end of the identification process are 0.27% and 3.18%, respectively. Moreover, in this case, the stator inductance estimation accuracy depends on the operating conditions of the SPMSM in which the AdNN3 is activated. The condition in which the minimum estimation error is achieved is the one at the rated current and 10% of the rated speed while the condition in which the maximum error is achieved is the one at 50% of the rated current and at the rated speed. At the end of the parameter identification process, the percentage estimation error is 0.16%. The percentage estimation errors associated with the filtered parameter estimationsψ * ave ,R * s,ave, andL * s,ave for the entire simulation time Bonfiglioli 170 are shown in Figure 15.

Conclusions
In this paper, a new multistep multiparametric identification method of SPMSMs highly suitable for cloud computing environment is presented. The results achieved show how the proposed solution overcame the rank deficiency problem without any signal injection and extra sensor, with the combined action of the stator inductance, resistance, and rotor PM flux linkage estimations. The R-static algorithm managed to identify the steady states of the SPMSM and of the parameter estimations. Moreover, the results of the parameter estimations confirm the validity of the convergence condition obtained by means of analytical studies and employed in the design of the DMA. We achieved final estimation errors of 0.23%, 0.35%, and 0.11% for the rotor flux linkage, stator resistance, and stator inductance, respectively, in a simulation environment. Finally, the criteria proposed for the automated tuning of the Box-Muller noise and of the AdNNs learning rates have also been validated by the achieved results.
We stress that the proposed method is highly suitable for the cloud computing implementation in large-scale applications since it requires a simple tuning of a few parameters. Moreover, note that the proposed algorithm can be partially reused for applications Finally, the algorithm runs on a cloud resource with 8 Intel Xeon Platinum 8259CL @ 2.50 GHz processors and 32 GB of RAM. Execution timings of multiple parallel runs of the identification procedure of the Bonfiglioli 170 are reported in Table 3. Note that the total execution time highlights the computational sustainability of the proposed scheme in large-scale applications using cloud computing resources. The achieved values are appropriate since for parallel runs less than the number of cores, the execution times are similar, i.e., approximately one core per run is dedicated. Instead, for parallel runs greater than the number of cores, the execution time proportionally increases as the number of parallel runs increases.

Conclusions
In this paper, a new multistep multiparametric identification method of SPMSMs highly suitable for cloud computing environment is presented. The results achieved show how the proposed solution overcame the rank deficiency problem without any signal injection and extra sensor, with the combined action of the stator inductance, resistance, and rotor PM flux linkage estimations. The R-static algorithm managed to identify the steady states of the SPMSM and of the parameter estimations. Moreover, the results of the parameter estimations confirm the validity of the convergence condition obtained by means of analytical studies and employed in the design of the DMA. We achieved final estimation errors of 0.23%, 0.35%, and 0.11% for the rotor flux linkage, stator resistance, and stator inductance, respectively, in a simulation environment. Finally, the criteria proposed for the automated tuning of the Box-Muller noise and of the AdNNs learning rates have also been validated by the achieved results.
We stress that the proposed method is highly suitable for the cloud computing implementation in large-scale applications since it requires a simple tuning of a few parameters. Moreover, note that the proposed algorithm can be partially reused for applications other than PMSMs. In particular, steady-state identification procedures of the motor and the estimated parameters can also be implemented for other motor typologies. Instead, the AdNNs and the convergence condition need to be modified since they are based on the mathematical model of the motor to be identified.
Future developments will include tests with motor data provided by real production plants, and the extension of the method to different machines, such as internal PMSMs, will be considered.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy issues.


Deep Unsupervised Embedding for Remotely Sensed Images Based on Spatially Augmented Momentum Contrast

Convolutional neural networks (CNNs) have achieved great success when characterizing remote sensing (RS) images. However, the lack of sufficient annotated data (together with the high complexity of the RS image domain) often makes supervised and transfer learning schemes limited from an operational perspective. Despite the fact that unsupervised methods can potentially relieve these limitations, they are frequently unable to effectively exploit relevant prior knowledge about the RS domain, which may eventually constrain their final performance. In order to address these challenges, this article presents a new unsupervised deep metric learning model, called spatially augmented momentum contrast (SauMoCo), which has been specially designed to characterize unlabeled RS scenes. Based on the first law of geography, the proposed approach defines spatial augmentation criteria to uncover semantic relationships among land cover tiles. Then, a queue of deep embeddings is constructed to enhance the semantic variety of RS tiles within the considered contrastive learning process, where an auxiliary CNN model serves as an updating mechanism. Our experimental comparison, including different state-of-the-art techniques and benchmark RS image archives, reveals that the proposed approach obtains remarkable performance gains when characterizing unlabeled scenes since it is able to substantially enhance the discrimination ability among complex land cover categories. The source codes of this article will be made available to the RS community for reproducible research.

semantic variety of RS tiles within the considered contrastive learning process, where an auxiliary CNN model serves as an updating mechanism. Our experimental comparison, including different state-of-theart techniques and benchmark RS image archives, reveals that the proposed approach obtains remarkable performance gains when characterizing unlabeled scenes, since it is able to substantially enhance the discrimination ability among complex land cover categories. The source codes of this paper will be made available to the RS community for reproducible research.

Index Terms
Deep learning, metric learning, self-supervised learning, unsupervised learning, remote sensing, scene characterization.

I. INTRODUCTION
With the growing development of deep learning (DL) technologies, these kinds of methods have achieved tremendous success in many important Remote Sensing (RS) applications [1], [2], such as scene classification [3]- [7], object localization [8]- [12] and change detection [13]- [15], owing to their prominent capabilities to uncover highly representative features from RS scenes [16]. In general, DL techniques aim at projecting the visual content of input images onto a particular label space, using a hierarchy of nonlinear layers to generate a high-level semantic abstraction that is very useful to characterize RS data. Most available DL-based image characterization methods in the RS field rely on a supervised learning scheme, in which a large amount of labeled scenes is required to properly train the models and prevent over-fitting [17].
However, the task of obtaining relevant annotations for vast volumes of RS data can be very difficult and time consuming. This may severely constrain the applicability and potential of the supervised DL paradigm in operational RS environments, especially under the most challenging conditions [18].
In order to mitigate the need for labeled RS data, different strategies have been effectively explored in the literature, e.g. [19]- [22]. One of the most popular schemes is based on the use of pre-trained convolutional neural networks (CNNs) [23], where different pre-defined CNN architectures (e.g. AlexNet [24], VGGNet [25], GoogleNet [26] and ResNet [27]), trained on large-scale computer vision datasets (e.g. ImageNet [28]) are directly used as feature extraction methods for RS data. Despite its remarkable success [29]- [32], the existing limitations on the number of spectral bands and the data complexity make this transfer learning scheme unable to fully exploit the advantages of RS imagery [33]. An attractive option to relieve these limitations consists of using unsupervised methods to characterize unlabeled RS scenes. As a result, different methods have been successfully proposed within and outside the DL field [34]- [37]. However, the general unsupervised framework is often unable to introduce appropriate prior knowledge about the RS image domain, eventually constraining the resulting performance. Although a recently developed deep metric learning method -the Tile2Vec [38]-is certainly able to obtain promising results by using the geospatial information as prior knowledge, the unprecedented availability of massive RS archives, together with the constant development of the acquisition technology, still make unsupervised DL-based image characterization a major challenge in RS.
Note that the integration of the unsupervised mode into the deep metric learning approach [39] is highly limited by the contrasting land cover types that can be sampled in a single batch, which may eventually reduce the capacity of the model to distinguish a broader range of complex RS categories, and also motivates the development of novel techniques useful to deal with the large-scale variance complexity of the RS image domain [40].
With all these considerations in mind, this paper proposes a new unsupervised deep metric learning approach, called spatially augmented momentum contrast (SauMoCo), which has been specially designed to characterize unlabeled RS scenes. Inspired by Tobler's first law of geography [41], the proposed approach provides a new perspective on unsupervised land cover characterization in which not only the semantic similarities among nearby scenes are exploited to learn the corresponding feature embeddings, but also the inherent diversity within RS semantic concepts. To achieve this goal, we define a spatial augmentation criteria to uncover enhanced semantic relationships among RS tiles for the embedding space. Then, we build a queue of deep embeddings, where the size of queue is forced to be larger than the batch size in order to further increase the semantic variety of contrasting land cover tiles during the training process. Moreover, we introduce an auxiliary CNN into our model to consistently update the deep embeddings of the RS tiles in the queue. With the objective of validating the proposed approach, we conduct a comprehensive experimental comparison, using two benchmark datasets and different state-ofthe-art characterization techniques, which demonstrates the superior performance of the presented method in the task of categorizing RS scenes without using any land cover class information.
In short, the main contributions of this paper can be summarized as follows: 1) We propose a new unsupervised deep metric learning model (SauMoCo) to characterize unlabeled RS images. The presented approach pursues to exploit not only the semantic similarities among nearby geospatial locations, but also the inherent diversity within land cover concepts, by using a newly defined spatial augmentation criteria with a contrastive loss formulation and a momentum update-based optimization.
2) We investigate how the proposed SauMoCo model performs with large-scale training data, which gives us important insights about the working mechanism and practical advantages of the proposed method with respect to other unsupervised RS image characterization techniques available in the literature. The codes of this work will be released for reproducible research inside the research community.
The rest of this paper is organized as follows. Section II reviews some related works on RS scene characterization while highlighting their main limitations. Section III details the proposed unsupervised deep metric learning model for RS images. Section IV presents the experimental part of the work. Finally, Section V concludes the paper with some remarks and hints at plausible future research lines.

II. RELATED WORK
Different strategies have been successfully adopted within the RS field to relieve the need for labeled data when characterizing aerial scenes. This section reviews some of the most relevant trends, including pre-trained (Section II-A), unsupervised (Section II-B) and deep metric learningbased (Section II-C) methods. Additionally, we also analyze their main limitations in the context of RS problems (Section II-D).

A. Pre-trained Methods
One of the most popular schemes to characterize RS data is based on the use of pre-trained convolutional neural networks (CNNs). In more details, these approaches make use of pre-defined CNN models, such as AlexNet [24], VGGNet [25], GoogleNet [26] or ResNet [27], which are pre-trained on large-scale computer vision datasets, such as the ImageNet [28] collection. In this way, the amount of labeled RS scenes can be substantially reduced by transferring the knowledge from the standard image domain to the RS field. For instance, Hu et al. define in [42] two different schemes to take advantage of the VGGNet model pre-trained on ImageNet. In the first scheme, the authors employ the last fully connected layers as image descriptors. In the second one, an additional encoding procedure is used to fuse the last convolutional feature maps.
In both cases, a support vector machine (SVM) is adopted to finally classify the RS images.
Precisely, Marmanis et al. analyze in [29] the effectiveness of classifying remotely sensed scenes using different CNN-based representations transferred from ImageNet. Similarly, Li et al. present in [31] a multi-layer feature fusion framework which integrates several pre-trained DL models for RS scene classification. Zheng et al. build in [43] a holistic representation of RS images using a multi-scale pooling over pre-trained features. Kang et al. also combine in [44] several pre-trained CNN architectures to define a building-instance level land-use classification framework. In [30], the authors propose using a sparse autoencoder (AE) over pre-trained features to generate the final representation of RS scenes.
Despite the remarkable performance achieved by these and other pre-trained models, there are still some important limitations that substantially reduce the applicability of such transfer learning strategies within the RS field. On the one hand, standard image collections, such as ImageNet, are often made up of RGB imagery, which makes existing pre-trained networks unable to take advantage of the additional spectral bands provided by air-borne and space-borne optical sensors [45]. Note that RS instruments are often designed to provide valuable information outside the visible spectrum, and these data are essential in may important applications, such as biophysical parameter analysis [46] and land cover material study [47]. On the other hand, standard images often contain natural object-centric photographs that hardly represent the complexity of RS scenes, comprising fully focused multi-band shots of the Earth surface with plenty of complex spatio-spectral details within the same acquisition frame [48]. Precisely, these important differences often make necessary to consider broader strategies than pre-tranined DL models.

B. Unsupervised Methods
A more general option to relieve the need for annotated RS data is based on using unsupervised image characterization methods. In more details, these techniques work for characterizing aerial scenes without using any class label information, which becomes particularly attractive in RS problems [18]. Consequently, different unsupervised models (both within and outside the DL field) have been proposed to learn informative representations from unlabeled RS scenes. For instance, Cheriyadat presented in [49] a feature learning approach for aerial scenes which adopts a sparse coding framework to generate unsupervised data representations based on a set of basis functions derived from low-level measurements. Following this idea, other authors proposed using different unsupervised decomposition frameworks instead. This is the case of the works presented in [37], [50], [51], which make use of probabilistic topic models to represent the RS data as probability distributions of feature patterns. In [52] [34] a greedy layer-wise unsupervised pre-training method for learning sparse features from aerial images. In the case of [36], the authors define a shallow weighted deconvolution network for extracting features from RS scenes by minimizing the Euclidean distance between the original and the reconstructed images. Alternatively, some works in the literature also show the utility of convolutional generative adversarial networks to characterize standard and remotely sensed imagery [53], [54].

C. Unsupervised Deep Metric Learning
Notwithstanding the positive results of these and other important unsupervised methods, all these works mainly rely on generic clustering, decomposition or encoding procedures that are often unable to introduce relevant prior knowledge about the RS domain without using supervised information. Among all the conducted research, one of the most promising trends to adequately characterize RS images is based on the so-called deep metric learning approach [39]. In particular, deep metric learning aims at learning a low-dimensional metric space based on CNN models, where the feature embeddings of similar images should be close, and those of dissimilar images should be separated. Despite its great potential in RS problems [55]- [58],

D. Current Limitations in RS
Certainly, the Tile2Vec algorithm sets a path for learning more informative CNN-based characterizations of RS data from a completely unsupervised perspective. However, the task of generating highly meaningful representations of aerial scenes without using any kind of class label information still remains a very important challenge in RS [59], [60]. The recent availability of massive data archives, together with the constant development of the airborne and space acquisition technology, are steadily increasing the complexity of RS data and, hence, their semantic understanding. Precisely, this growing complexity often produces a huge within-class diversity and between-class similarity that introduce important limitations within the aforementioned learning scheme [39]. When integrating the unsupervised mode into the deep metric learning approach [38], it is logically necessary to train the CNN model by sampling negative RS image tiles within each batch. However, this strategy significantly reduces the capacity of the model to distinguish between a broader range of contrasting land cover types, since the learning process is constrained by the tiles that can be sampled in a single batch. Note that this point may become particularly critical in complex large-scale archives, which stimulates the development of more advanced unsupervised characterization techniques within the RS field [18], and ultimately motivates the research conducted in this work.

E. Novelty of the Proposed Method
To address all these challenges, this paper proposes a new unsupervised deep metric learning model that jointly exploits two different aspects: a spatially augmented contrastive loss and a momentum update-based optimization. In contrast to Tile2Vec [38], the proposed approach integrates a new spatial augmentation criteria that allows considering not only semantic similarities among nearby RS scenes but also the inherent semantic diversity of land cover concepts when learning the corresponding metric space in an unsupervised fashion. Note that this within-class variability has not yet been exploited in the context of characterizing unlabelled RS scenes despite the fact it may become very useful to relieve the large-scale variance problem of RS data [40]. Using this methodological improvement, the proposed approach is able to avoid the triplet loss limitations with scalable data while also taking advantage of additional contrastive RS image pairs during training. In order to further improve such contrasting land cover variety, the proposed approach also utilizes a momentum update-based optimization [61]. The general idea behind the momentum update is based on managing a dynamic dictionary of encoders to enhance the contrastive learning process. Following this idea, we build a queue of deep embeddings of RS scenes in order to force the length of such queue to be larger than the mini-batch size.
Unlike the standard momentum scheme which shows limited results with large-scale data [61], the proposed end-to-end approach is designed to exploit vast unlabeled RS archives by using a CNN-based backbone architecture to jointly characterize land cover scenes and update the queue. Compared with different state-of-the-art methods to characterize unlabelled RS scenes, the proposed approach is able to achieve a better performance than the methods in [29], [38], [54], [62], which also reveals the novelty and advantages provided by this work for the RS community.

III. SPATIALLY AUGMENTED MOMENTUM CONTRAST
Our newly proposed end-to-end unsupervised deep metric learning model for characterizing unlabeled RS scenes (SauMoCo) can be summarized in the following three parts: • A backbone architecture (called anchor CNN) which is used to generate the corresponding feature embedding of the input RS scenes. Note that this CNN architecture can be defined according to a specific off-the-shelf topology, such as AlexNet [24], VGGNet [25], GoogleNet [26] and ResNet [27].
• A spatially augmented loss, based on the contrastive loss formulation and a newly defined spatial augmentation criteria, which exploits not only the semantic similarities among nearby RS scenes but also the inherent diversity within land cover semantic concepts.
• The corresponding optimization algorithm, which learns the proposed model parameters using a momentum contrast update. To achieve this goal, a queue of deep embeddings is constructed, and an additional CNN model (called momentum CNN) is introduced to update such queue. It is important to highlight that this network should be defined with the same architecture with regards to the one of the anchor CNN for a scalable training. Figure 1 illustrates in a graphical manner the proposed unsupervised deep metric learning framework. In the following sections, we detail the newly defined loss function and the considered optimization algorithm.

A. Spatially Augmented Loss
x i , an anchor patch x a i (located in its center) can be cropped with a certain size W × W . With a certain distance d of the anchor patch x a i , a neighborhood patch cropped from x i is defined as its spatial augmentation, which is x n i . If the distance is 100 pixels (in both vertical and horizontal directions), the center of x n i should be within 100 pixels with respect to the center of x a i . Let f a i ∈ R D denote the deep embedding of x a i obtained by a CNN model F(·; θ) on the unit sphere (i.e., f a i = F(x a i ; θ)/ F(x a i ; θ) 2 ), where D is its dimension and θ represents the parameters of the CNN model. We identify this model as anchor CNN.
As noted by the first law of geography [41], everything is related to everything else, but nearby things are more related than distant things. Following this rule, the proposed method relies on the assumption that images that are geographic neighbors should be semantically more similar than distant images [38]. Therefore, the embeddings of nearby images should be closer than those of distant images in the metric space. However, it is important to highlight that the proposed spatial augmentation criteria is different from the one considered in other works, such as in [38]. Specifically, we do not fix the position of the spatially augmented patches to a specific neighbour position, but to a neighbourhood region of the anchor patch. As a result, our augmentation criteria allows certain spatial variations on the cropped areas in order to increase the variety of spatially augmented patches that are extracted on-line. To achieve this in a scalable way, we adopt a contrastive learning mechanism [63], [64], where the contrastive loss of x a i can be defined as: In this equation, the inner product f a i , f n i measures the cosine similarity between the embedding f a i of the anchor patch x a i and the one f n i of its spatial augmented patch x n i . Besides, τ represents a temperature parameter controlling the concentration level of the sample distribution [65]. Intuitively, Equation (1) describes the log-likelihood of the spatial augmented patch, which can be classified as its anchor patch among all the anchor patches in X . Then, the corresponding contrastive loss over the whole dataset can be formally expressed as: By optimizing Equation (2), we can obtain the deep embeddings of X and the trained CNN model, which is useful for characterizing unlabeled RS scenes and conducting the corresponding downstreaming land-cover categorization tasks.

B. Optimization via Momentum Update
In order to sufficiently train the CNN model based on Equation (2) in an unsupervised fashion, a scalable dataset is logically required to be fed into the deep model. For scalable datasets, how to sufficiently sample the negative patches, (i.e., x a j ) with respect to x a i should be carefully defined, since the number of spatially augmented patches and their consistency are both critical aspects within the proposed contrastive unsupervised learning scheme.
One common strategy adopted in the literature is based on sampling the negative patches within each mini-batch [39]. However, this optimization mechanism has important limitations for training our deep model with scalable spatially augmented data. In more details, this minibatch sampling process assumes that each patch can be seen once during one epoch of training and hence x a i only exists in one mini-batch for the current iteration. Consequently, the CNN model is only able to see its corresponding negative patches x a j belonging this mini-batch, while other important samples outside the mini-batch cannot be considered. Precisely, this fact can substantially reduce the semantic variety of contrasting land cover types during training, which is a key factor to allow learning more informative RS image representations from an unsupervised perspective.
To solve this problem, we adopt the momentum update rule [61], [66] for training our newly proposed unsupervised RS image characterization model. Specifically, a queue of the deep embeddings of image patches x a j is constructed, where the size of the queue is forced to be larger than that of the mini-batch. In this way, the unsupervised learning process can be substantially enhanced by considering contrasting patches beyond a single batch. During the training phase, the embeddings of the current mini-batch are compared with the ones in the queue, as they are progressively replaced. The embeddings of the current mini-batch are enqueued and the oldest ones are dequeued. Moreover, in order to consistently update the deep embeddings in the queue, an auxiliary CNN model with parameter set θ aux is introduced. We identify the θ aux model as momentum CNN and it is updated as follows: where m ∈ [0, 1) is a momentum coefficient. It is worth noting that only the CNN with θ is updated by means of back-propagation. The momentum CNN with parameters θ aux can be evolved more smoothly than the CNN with θ. Then, the embeddings in the queue (encoded by the momentum CNN) are updated by:f In this expression,f i denotes the features generated by the momentum CNN. In other words, the embeddings in the queue are replaced by the ones encoded by the momentum CNN after each training epoch. To this end, the proposed optimization mechanism is detailed in Algorithm 1.

A. Dataset Description
In this work, we use two benchmark RS image datasets to validate the effectiveness of the proposed method. A detailed description of the datasets is provided below: Sample a mini-batch of x i .

4:
Within x i , randomly generate the center pixels of x n i .

5:
Crop the patches x n i and x a i online. Crop, Residential, River, and Sea Lake. The EuroSAT archive is also publicly available 2 .
These two RS archives have been selected to evaluate the performance of the unsupervised RS image characterization process from a single-source land cover acquisition perspective because  Figure 2(a). Then, we randomly select a total of 100,000 images (with a size of 250 × 250 pixels) from the downloaded tiles.
2) In the case of EuroSAT, we downloaded 100 Sentinel-2 Level-1C image products which have been globally sampled from the entire globe. Figure 2(b) illustrates the geo-locations of the downloaded Sentinel-2 products. Then, we select 100,000 random images (with a size of 264 × 264 pixels) from the downloaded products. centers is 100 pixels, according to the defined spatial augmentation criteria.

B. Experimental Setup
The proposed method is implemented in PyTorch [68].  [54], 2) MARTA GAN [62], 3) the ResNet18 model pre-trained on ImageNet while considering the most discriminating principal components (pretrained CNN+PCA) [29] and 4) the Tile2Vec [38]. In the case of the pre-trained CNN+PCA, it is important to highlight that we make use of the PCA method after extracting the pre-trained features to generate the corresponding deep embeddings with the same dimensionality as the other methods. All the experiments are conducted on an NVIDIA Tesla P100 graphics processing unit (GPU).
To measure the effectiveness of the proposed approach as compared with the other methods, we extract the deep embeddings for the NAIP and EuroSAT collections after training. Then, we use the available annotations to compute the corresponding classification results for each dataset.
In more details, we provide five different experiments for validating and analyzing the results from several perspectives: 1) Evaluation of deep embeddings based on random forest (RF) classification: We first utilize the random forest (RF) classifier to measure the classification performance based on the extracted feature embeddings of the two datasets obtained by the considered methods. For each dataset, we randomly select 80% images for training the classifier and evaluate its performance on the rest 20% images. In order to obtain a mean score of the overall accuracy, a total of 100 trials are conducted. Then, we calculate the mean and standard deviation values of the obtained accuracy scores.
2) Visualization of image retrieval: In this experiment, we conduct a retrieval test to explore, from a qualitative perspective, the performance of the considered characterization methods. In particular, we extract one query image patch from a complete NAIP scene. Then, we use the pretrained CNN+PCA, Tile2Vec and SauMoCo models to obtain the deep embeddings of the selected patch as well as the rest of the patches in the scene. Finally, we calculate their corresponding similarity maps with respect to the query and retrieve the 10 nearest neighbor patches within the whole scene.      images are cropped on-line, which also provides additional advantages as data augmentation strategy. By doing so, a higher semantic variety of similar images can be generated during the training phase (with respect to the anchor images). In comparison, the triplet set utilized in Tile2Vec is constructed beforehand and it does not exhibit any data augmentation capability during the training. Therefore, the proposed approach can also take advantage of the proposed spatial augmentation criteria. To analyze the differences between Tile2Vec and SauMoCo in more details, Table II provides the corresponding class-wise F1 scores obtained by the RF classifier on the EuroSAT dataset,

3) Evaluation of CNN model initialization
where the two best results are highlighted in bold and gray-shaded font. As it is possible to observe, the proposed approach obtains the best and second best performances in 2 and 5 landcover classes, respectively. Although Tile2Vec also exhibits positive results in 5 categories, the performance for the rest of the classes is rather limited, being even worse than that of the pre-trained CNN+PCA and MARTA GAN in some cases. Precisely, these important differences make SauMoCo more stable and accurate from a global perspective, indicating that the proposed approach is able to extract more relevant information about a wider range of semantic concepts.
2) Visualization of image retrieval: As shown in Figure 5, we extract one query image patch from a NAIP tile. Then, we obtain its deep embedding and the ones of the rest of the patches in the tile based on the considered methods. Subsequently, we calculate their similarities and obtain the corresponding heatmaps for SauMoCo (b), Tile2Vec (c) and the pre-trained CNN+PCA (d), where brighter colors denote a higher similarity in the embedding space. Additionally, the 10 nearest neighbor patches are illustrated in (e). As it is possible to observe in the heatmaps, the locations of the most similar patches with respect to the query can be more clearly identified in (b). Precisely, these results indicate that the semantic information of the RS scene is not properly encoded based on Tile2Vec and the pre-trained CNN+PCA model, since there are larger parts in the image that are considered to be similar to the query in the embedding space. Regarding the nearest neighbor results, it is possible to see that the image patches retrieved from the embedding space generated by SauMoCo are the most visually similar with regards to the query. That is, the proposed approach is able to model the semantic content of the query more accurately than the other methods, since all the retrieved images display similar land-cover patterns.
3) Evaluation of CNN model initialization: Figure 6 shows the learning curves of the ResNet18 model over the EuroSAT collection when using two different initialization strategies: with the parameters obtained by the proposed approach and with the pre-trained ImageNet parameters.
According to the displayed results, it can be seen that the classification accuracy can be slightly improved when the parameters of the CNN model are initialized via the proposed approach.
Although the pre-trained ImageNet initialization exhibits higher classification accuracies at the beginning of the training process, SauMoCo is able to consistently achieve better results after 30 epochs. This fact reveals that the proposed approach is able to capture richer semantic information in the corresponding embedding space, since a better minimum location of the loss function can be discovered by the parameters pre-trained via SauMoCo.  ResNet50 is slightly improved in comparison with ResNet18. As shown in Table III, compared with ResNet18, the classification accuracy based on the deep embeddings from ResNet50 can be improved with a score of 0.5% on the NAIP dataset.

V. CONCLUSIONS AND FUTURE LINES
This paper presents a new unsupervised deep metric learning framework (SauMoCo) to characterize unlabeled RS scenes. Specifically, the proposed approach initially defines a spatial augmentation criterion to uncover semantically similar RS images based on the the first law of geography. Then, a queue of deep embeddings is built such that the size of queue is forced to be substantially larger than the batch size, to improve the semantic variety of contrasting land cover types during the training. To achieve this goal, an auxiliary CNN model is also used to consistently update the deep embeddings in the queue. The experimental part of the work, conducted over two benchmark datasets and based on the use of different characterization methods, reveals that the proposed unsupervised deep metric learning model is able to provide competitive advantages with respect to other state-of-the-art techniques in the task of representing unlabeled RS images.
One of the main conclusions that arises from this work is the relevance of considering a broader variety of land cover types when learning unsupervised RS image characterizations.
In this regard, the proposed approach takes advantage of the defined spatial augmentation criteria and the considered queue of deep embeddings to enrich the semantic information of different RS categories during the contrastive learning process. Precisely, this feature allows our SauMoCo to enhance the global discrimination ability among unsupervised land-cover classes, and also to provide a more robust behavior with different datasets and settings. Owing to the remarkable performance achieved by the presented method, our future work will be directed towards adapting it to inter-sensor data and other important RS tasks, such as dimensionality reduction of hyperspectral imagery or fine-grained land-use categorization.