Application of the Singular Spectrum Analysis on electroluminescence images of thin-film photovoltaic modules

Singular Spectrum Analysis (SSA) is a model-free time series analysis method that belongs to the so-called subspace methods, Van Der Veen et al. (1993). In subspace methods, a signal estimation is performed by taking a certain linear subspace. SSA can also be considered as a low-rank approximation method, Markovsky (2018).

The general theory of SSA for one-dimensional time series is elaborated in Golyandina et al. (2001). Golyandina et al. (2018) provide more updated information on extensions with a strong focus on applications.

The output of SSA-like methods is a decomposition of an observed signal $x$ (e.g. a time series, a multivariate time series or an image) into a sum of identifiable components:

Among all SSA-like methods, the following four common algorithm steps can be isolated.

1. 1.

   Embedding. The original signal $x$ (e.g. time series or an image) is mapped into a matrix $\mathbf{X}$, that is called a trajectory matrix.

   The embedding is parametrised by a single parameter denoted by $L$ (a number or a vector).

2. 2.

   Decomposition. The second step consists of the decomposition of the trajectory matrix $\mathbf{X}$ into a sum of matrices of rank 1.

   Often the singular value decomposition (SVD) is used for this purpose, which is an optimal rank-one matrix decomposition in the Frobenius norm sense.

3. 3.

   Grouping. The third step is a grouping of the decomposition components. At the grouping step, the elementary rank-one matrices are grouped and summed within groups.

   The grouping algorithm step is often semi-automatic and depends on the type of data and application.

4. 4.

   Reconstruction. The grouped components are not necessarily valid trajectory matrices. Hence a projector operator to trajectory matrix space is applied on the grouped components, resulting in the final signal decomposition (1).

In this paper, we utilise the 2D-SSA variant of the algorithm, Golyandina and Usevich (2009, 2010); Golyandina et al. (2015). For this variant, the trajectory matrix is a Hankel-block-Hankel matrix and the corresponding embedding is parametrised by a two-dimensional parameter $L=(L_{1},L_{2})\in\mathbb{N}^{2}$. The decomposition is performed with SVD, and the reconstruction projection operator is a 2-step diagonal averaging procedure. The precise forms of the embedding, SVD, and reconstruction projection operator are provided in the Appendix.

The grouping step of the SSA algorithm determines the form of the final signal decomposition. A typical SSA decomposition of a signal is the decomposition into a slowly-varying trend, regular oscillations, and noise. An important notion in the SSA theory is the notion of the signal of finite rank. This notion allows us to classify and group together rank-one matrices into a trend, oscillations, and noise components.

Informally, finite rank signals are those that have a trajectory matrix of a finite rank (see a formal definition in the Appendix). 2D-SSA produces a class of signals of specific objects of finite rank. Those objects have the form described in the following theorem, Golyandina and Usevich (2009); Golyandina et al. (2015).

Estimation of Signal Parameters via Rotational Invariance Techniques (ESPRIT) is a method to estimate parameters of a mixture of one-dimensional, Roy and Kailath (1989), and two-dimensional amplitude-modulated sinusoids, Rouquette and Najim (2001); Wang et al. (2005), in background noise.

For the 2D-ESPRIT method, Rouquette and Najim (2001), the observed data $y$ is generated by the following additive model:

where $0\leq m\leq N_{x}-1$ and $0\leq n\leq N_{y}-1$, $\varepsilon$ is the zero-mean Gaussian noise with variance $\sigma^{2}$. The model for the signal $x$ is given by the sum of amplitude modulated two-dimensional sinusoid

where $\omega_{1k},\omega_{2k}$ are the normalised frequencies in different directions, $\alpha_{1k},\alpha_{2k}$ are the damping factors, $s_{k}$ amplitudes, and $\phi_{k}$ phases.

By Theorem 2.1 the signal (4) is of finite rank. The ESPRIT methods utilises the fact of the rank-deficiency of the trajectory matrix of the observed signal $y$, and a certain transformation matrix between sub-trajectory matrices is obtained. The frequency and damping factor parameters are computed from the argument and absolute values of the complex-valued eigenvalues of the obtained transformation matrix. See more details in Roy and Kailath (1989); Rouquette and Najim (2001).

The ESPRIT estimates the frequencies $\omega_{\cdot k}$ and the damping factors $\rho_{\cdot k}$ of the signal model (4). Here only the amplitudes $s_{k}$ and phases $\phi_{k}$ remain to be estimated. This problem can be reformulated as a linear regression model using the formula of the cosine of sums. With this formula, (4) can be rewritten as

where $\rho_{\cdot k},\omega_{\cdot k}$ are the parameters estimated from the ESPRIT, and $A_{k},B_{k}$ are the parameters of the linear model to be estimated.

Note that the dependent variable in the linear regression model are the values of $x(n,m)$, the finite-rank signal extracted from $y$. In terms of the SSA algorithm this corresponds to the series reconstructed from the selected components.

The amplitude and phase are given by:

For our application, we use an R-package “Rssa”, Golyandina et al. (2018); Korobeynikov (2010); Golyandina and Korobeynikov (2014); Golyandina et al. (2015), where all the required functionality including the ESPRIT method is implemented. It should be noted that the trajectory matrix for 2D-SSA is large and has $O(N^{2})$ number of elements, where $N$ is a number of pixels in the signal image. However, the structure of the Hankel-block-Hankel matrix allows implementing of SVD with Lanczos algorithm efficiently in time and memory by computing product of a matrix and a vector with Fast Fourier Transform, Korobeynikov (2010); Lanczos (1950).

For computing amplitude and phase parameters, we utilise a standard R-function “lm”, Chambers (1992).

The form of the parametric model of the SSA signal suggests that Fourier analysis can be used to obtain similar results. However, in order to obtain a compact representation of a signal a small Fourier coefficient should be discarded. For that purpose a sparse Fourier analysis approach can be used, Hassanieh et al. (2012).

However, the parametric model of the SSA signal is more flexible, as every periodic has an amplitude modulation. Furthermore, the Fourier analysis is a low-resolution type method, as in the context of time series, frequency can be estimated only up to $1/N$, where $N$ is the length of the time series. Whereas ESPRIT is regarded as a high-resolution method.

The grouping step of the SSA algorithm allows combining of decomposed rank-one components into groups. We define these groups that result in EL image decomposition onto 3 components: global intensity variation, cell, and aperiodic components.

Figure 2 shows a close up image of the module, that clearly identifies a set of 10 parallel cells separated by interconnection lines. The EL intensity varies systematically over the cell width. This is the result of the series resistance of the electrodes. As the transparent zinc-oxide front electrode exhibits a much larger sheet resistance than the molybdenum back contact, the voltage over the diode junction drops from one side to the other Helbig et al. (2010); Pieters and Rau (2015). As the EL intensity depends primarily on the cell voltage this leads to a clear intensity gradient over each\linelabelphysical reason for periodic cell.

To achieve good separability of the periodic components, the parameter $L$ of the SSA embedding is selected to equal approximately half of the dimensions of the input image, Golyandina et al. (2001, 2018). The ESPRIT estimates frequencies and damping factors of the components, which are used in the decision of the grouping step.

A thin-film module consists of a fixed number of interconnected cells, therefore, the cell components can be identified as periodic components that have a period smaller than $\frac{\text{image width}}{150}$. The other components are grouped together as the global variation and the residual image captures information of non-low rank signals, such as aperiodic features of an EL image.

The image decomposition algorithm steps are described in the Algorithm 1.

The choice of the 50 computed components is arbitrary, as 10 components incorporate $99.9\%$ of the Frobenius norm of the trajectory matrix $\mathbf{X}$. The RMSE of the linear regression model in the ESPRIT indicates the accuracy of the model (4). For a set of 50 EL images, the mean RMSE equals 0.033.

Figure 3 depicts the non-cell components, the cell components, and the residual image. It can be argued that different components have a different physical origin. By its nature, the global variation component describes changes in a material which results in large losses that spread over large portions of a module. The cell component variation is influenced by the EL measurement conditions. Lastly, the aperiodic component captures effects caused by non-regular changes in material like shunts, or droplets (see Sovetkin et al. (2021)).\linelabelphysical interpretation of components

Shunts are characterized by a more conductive connection between the front and back electrodes than the normal solar cell structure (i.e. the solar cell structure is damaged or missing). There are many causes for shunts. Commonly shunts originate from debris of the copper evaporation source or pinholes in the CIGS absorber Misic et al. (2015); Misic (2017). Shunts are generally relevant to the solar module performance, in particular under low light conditions Weber et al. (2011).

In addition to shunts we noticed the CIGS modules often exhibit “droplets” in the EL images. The appearance of droplets resembles water stains and thus we speculate these structures originate from the chemical bath deposition. At this point it is unknown what the impact of droplets is on the module performance, however, the bright appearance imply a local change in quantum efficiency according to the reciprocity relations between luminescence and quantum efficiency Rau (2007).

We remark that decomposition can be applied to an image (additive model assumption\linelabeladditive and multiplicative models) as well as to the logarithm of an image (multiplicative model assumption). The logarithm of an image corresponds to the internal voltage (see (8), Section 4.3).

Furthermore, we remark that it is important to correct any perspective distortion present in images, as a periodic image distorted in such a way is no longer a signal with a finite rank in the settings of the 2D-SSA methodology.

In order to identify interconnection lines, it is sufficient to locate the global minima for every level in the normal direction of the interconnection lines.

The ESPRIT model satisfies equation (4), which is a sum of amplitude-modulated cosine functions. A derivative of such function is again a sum of polynomial amplitude-modulated cosine functions, similar to the general form of the signal of finite rank (2). Such derivatives can be computed symbolically. Hence all the global minima in the normal direction of the interconnection lines can be identified by evaluating precise values of the derivatives, and filtering out points that satisfy a minima extreme point requirement.

Algorithm 2 describes the steps of the interconnection line identification. Symbolic derivatives are computed using the “Deriv” R-package, Clausen and Sokol (2019).

We remark that the resulting expression for the cell component signal, as well as its derivatives, can be computed on a finer grid $P$ than the original pixel coordinates. Hence the interconnection line identification is performed with sub-pixel accuracy.

Figure 4 visually demonstrates the steps of the method. Figure 4a displays the estimated cell component. Figure 4b shows a slice of a pixel value intensities in Figure 4a, where computed local minima identified with red dots. Lastly, Figure 4c shows the module with its interconnection lines identified by red lines.

Another application of the obtained parametric expression for a module is an estimation of the inverse characteristic length $\lambda$ Helbig et al. (2010); Pieters and Rau (2015).

In a defect-free PV module, there are no currents in the vertical direction (along the interconnection line direction), and thus current satisfies the following linearised 1D Poisson equation, Helbig et al. (2010); Augarten et al. (2016).

where $\lambda=\sqrt{\frac{R_{\mathrm{sheet}}}{r_{\textrm{j}}}}$ is the inverse characteristic length, where $R_{\mathrm{sheet}}$ is the sheet resistance and $r_{\textrm{j}}$ is the local differential junction resistance. Note, that (7) is an equation between two fields, where $V(x,y)$ and $\lambda(x,y)$ depend on the position $(x,y)$ in a PV module.

The luminescence intensity related\linelabelintensity and voltage to the internal voltage via the following relation

where the constant $c_{0}=\frac{q}{kT}$ is the thermal voltage (with the elementary charge, $q$, Boltzmann’s constant, $k$, and the temperature, $T$), $c$ is a parameter that describes the optical system from the photon generation in the solar cell absorber material to the photon detection within the camera. As such the parameter $c$ depends on the quantum efficiency of both the solar cell and the used camera including all optical components, and the spectral photon density of a black body Rau (2007). Generally, the constant $c$ may vary over the solar cell area due to the camera optics and variations in the solar cell properties. However, in general, the variations in $c$ are small compared to the exponential voltage-dependency Abou-Ras et al. (2016). In our analysis, we assume that $c$ is constant over the whole module area, and is accessible to a researcher.

The decomposition of a signal onto trend and cell components provides us an EL image signal without small aperiodic defects like shunts (small dark areas within cells boundaries). The global intensity and cell components without small aperiodic defects can be considered as a module without defects. Hence, we consider a module without defects to be an image $G+S$, of the output of Algorithm 1.

Hence combining (7) and (8), allows us to express inverse characteristic length as a function of intensity and its derivative:

and therefore,

All derivatives are computed symbolically from the estimated cell component signal.

In order to achieve higher image resolution, several EL images can be stitched together. However, the image aligning can be imperfect, as shown in Figure 5a. These misaligned image patches and the resulting stitch line can be attributed to incorrect perspective as well as radial distortions. The latter distortion leads to a non-linear transformation that is needed to be applied to an image for the stitch line correction.

Algorithm 3 proposes an approach to correct such distortion. The basic idea of the algorithm is the estimation of phase shifts in neighbouring locations, where each shift is estimated by application of MSSA, Golyandina et al. (2018), in a direction perpendicular to the interconnection lines.

The result of the algorithm is a displacement map that defines shifts in the horizontal direction. Note that shifts have sub-pixel accuracy.

Figure 5 shows the result of the algorithm. Figure 5a shows a patch of the original EL image with a stitched part. Figure 5b depicts the shifted image, resulting from applying the displacement map. As each displacement is estimated locally between neighbours, the resulting transformation is a non-linear transformation.

To evaluate the accuracy of the phase shift estimation we model the following two time series

where $x\in\{1,2,\ldots,1000\}$, periodics with period $20$ and $30$ correspond to the signal (models the cell component), periodics with period $50$ and $70$ slowly varying trend, and $\varepsilon_{1},\varepsilon_{2}$ are two independent Gaussian iid processes with zero-mean and unit variance.

The signal part of the series $s_{1}$ is shifted by $7$ units relative to the signal of the series $s_{2}$. Figure 6 shows a part of those series on the interval $[1,100]$.

Table 1 shows the accuracy of the shift estimation of the selected signal components. That simulation was performed using 100 repetitions.

Following notations of Golyandina and Usevich (2010), let $x$ be an 2D-array with dimensions $(N_{x},N_{y})\in\mathbb{N}^{2}$:

The 2D-SSA embedding is defined by the window size vector $(L_{x},L_{y})$, which is restricted by $1\leq L_{x}\leq N_{x},1\leq L_{y}\leq N_{y}$, and $1<L_{x}L_{y}<N_{x}N_{y}$. Let $K_{x}\coloneqq N_{x}-L_{x}+1$ and $K_{y}\coloneqq N_{y}-L_{y}+1$, then the trajectory matrix of 2D-SSA is given by the following Hankel-block-Hankel matrix:

where each block

By construction there is a one-to-one correspondence between 2D-arrays of size $N_{x}\times N_{y}$ and the Hankel-block-Hankel matrices.

Let $\mathbf{Z}$ be an arbitrary matrix with a block-structure, where each block $\mathbf{Z}_{i}$ has the same dimension as the matrix $\mathbf{X}_{i}$

Then a projection of $\mathbf{Z}$ to Hankel-block-Hankel matrix can be computed in two steps. Firstly, diagonal averaging is performed within blocks $\mathbf{Z}_{i},i\in 0,\ldots,N_{y}-1$. Secondly, the blocks of the matrix $Z$ are averaged between themselves. The projection can be applied in the reverse order as well.

Let $\mathbf{S}=\mathbf{X}\mathbf{X}^{T}$, $\lambda_{1}\geq\ldots\geq\lambda_{d}>0$ be non-zero eigenvalues of the matrix $\mathbf{S}$, $U_{1},\ldots,U_{d}$ be the corresponding eigenvectors, and $V_{i}\coloneqq\mathbf{X}^{T}U_{i}/\sqrt{\lambda_{i}},i=1,\ldots,d$.

Then SVD of the matrix $\mathbf{X}$ can be written as

The values $\sqrt{\lambda_{i}}$ are the singular values of $\mathbf{X}$.

An image $x$ has rank $r$ if the rank of trajectory matrix $\mathbf{X}$ equals $r<\min(L_{x},L_{y},K_{y},K_{y})$. In other words, the trajectory matrix is rank-deficient. If rank $r$ does not depend on the choice of $L$ for any sufficiently large dimensions of $x$, then $x$ is called to have a finite rank.

Objects of finite rank are closely related to the linear recurrent sequences, Kurakin et al. (1995). Linear recurrent formulae are used to build forecast based on the SSA signal.

computational complexity

From the computational point of view, the hardest steps of the proposed algorithms are the singular value decomposition (SVD), ESPRIT (Section 2.2), and the phase and amplitude least-squares fit (Section 2.3).

In the context of the SSA application, the computational complexity of the SVD of a Hankel-block-Hankel matrix is $O(kN\log N+k^{2}N)$, where $N$ is the number of pixels in an image and $k$ is the number of computed eigentriples, Korobeynikov (2010). The linear equations solved for ESPRIT and the phase and amplitude parameters estimations are solved using QR-decomposition, and require $O(N^{3})$ operations. Our numerical experiments show that the SVD decomposition dominates the computational time for the proposed algorithms for image sizes with width and height less than 4000 pixels.

To measure the required time and memory of the proposed algorithms we utilise a machine with Intel(R) Xeon(R) CPU E5-1620 3.5GHz processor and 31GB of RAM. The time measurements are performed on a single CPU.

Figure 7 shows the time required to perform steps of Algorithm 1 for square images (width equals height) for different image widths. The red points correspond to the measured time in seconds, and the black line is a parabola fitted to the points. Table 2 shows the amount of memory required for the steps of Algorithm 1.

Algorithm 2 utilises results of the Algorithm 1 and requires additional symbolical differentiation. Such computation is also required for the inverse characteristic length (Section 4.3). The running time of the symbolical differentiation depends on the number of terms of the signal (4). Table 3 demonstrates the relationship between the number of terms in the signal and the running time required for the symbolic differentiation. Algorithm 2 and the inverse characteristic length computation works with the cell component that typically consists of 5–7 terms.

Algorithm 3 utilises a different version of SSA, that requires less time and memory, however, the algorithm computes MSSA multiple times for several pairs of neighbour rows. A single iteration of the algorithm loop requires 0.7 seconds and 0.78 of RAM for an image with width of 4000 pixels. Different iterations of the loop can be run in parallel. Note that most of the used memory is occupied by an image itself, and hence the memory can be shared by multiple processes. Algorithm 3 usually requires running about 100–200 iterations, as an approximate location of the stitching lines is known. Hence, the total running time an 8-core processor can be as little as 20 seconds.

Lastly, we remark the proposed algorithms run in a deterministic amount of time. Hence, the methods can be run in real-time applications.\linelabelreal time remark