A method and apparatus for joint equalization of PMD and RSOP impairments

Abstract

The application provides a PMD and RSOP damage joint equalization method and device, the method comprises the following steps: constructing a polarization damage joint equalization model according to the characteristics of polarization damage in a fiber link under a thunderstorm scene, and initializing SCKF parameters; time updating SCKF to predict damage values of damaged signals collected in a current sliding window; measurement updating SCKF according to a measurement equation of the sliding window to obtain a measurement prediction value; combining the radius distribution of received symbols to perform a circular ring decision operation on the measurement prediction value based on maximum posterior probability to obtain an optimized radius of a corresponding constellation circular ring, calculate an innovation matrix, and fuse a prediction estimation value and a measurement estimation value to obtain a posterior estimation value; the sliding window is moved forward to equalize the damage of signals in a next window. The application has excellent performance in a low OSNR scene, and has the advantages of high tracking accuracy, large initial error tolerance and fast convergence speed.

Description

Technical Field

[0001] This invention belongs to the field of optical fiber communication technology, and particularly relates to a joint equalization method and apparatus for PMD and RSOP impairments. Background Technology

[0002] First, during thunderstorms, for fiber-optic composite overhead ground wire (OPGW) links widely used in power communication backbone networks, the instant of lightning strike causes a sharp increase in polarization state rotation (RSOP), reaching up to 5.1 Mrad / s. Second, due to the strong correlation between RSOP evolution and polarization mode dispersion (PMD) impairment, the faster the RSOP, the faster the temporal correlation between PMD at different times fades, making it difficult for classical polarization demultiplexing algorithms based on multiple-input multiple-output (MIMO) structures (such as constant-mode / multi-mode algorithms (CMA / MMA)) to balance this combined polarization impairment. Third, as a breakthrough technology, probabilistic constellation shaping (PCS) approaches the Shannon limit in optical communication systems by changing the probability distribution of constellation points and is currently being widely applied in long-distance coherent optical communication systems. However, because PCS technology causes constellation points of different amplitudes to no longer follow an equiprobable distribution, many digital signal processing (DSP) algorithms developed for standard quadrature amplitude modulation (QAM) signals will experience severe performance degradation or even failure in PCS systems. Meanwhile, since the polarization demultiplexing (PolDemux) algorithm is located relatively early in the DSP process of the entire probabilistic constellation shaping (PCS) system, its algorithm performance directly affects the accuracy of subsequent frequency offset estimation (FOE) and carrier phase recovery (CPR). Therefore, joint equalization of RSOP and PMD damage in the probabilistic constellation shaping (PCS) system under thunderstorm scenarios has important theoretical research and practical application value.

[0003] Existing polarization demultiplexing algorithms suffer from numerous problems, including unsuitability for probabilistic shaping systems, inability to handle equalization pressures under thunderstorm conditions, and weak adaptability to higher-order modulation formats. For example, the blind polarization demultiplexing method in existing probabilistic constellation-modulated coherent optical communication systems selects circular radius symbols that occupy a large proportion of the PCS signal and performs least-squares plane fitting on these symbols in Stoke space to achieve polarization demultiplexing. However, this method has poor dynamic RSOP equalization capability and is only applicable to short-range transmission systems. Another existing polarization demultiplexing algorithm based on maximum likelihood independent component analysis (ML-ICA) estimates the inverse Jones matrix W using ML-ICA to maximize the mutual information between the X-polarization and Y-polarization inputs and outputs, thereby achieving polarization demultiplexing. However, this method has a low performance ceiling and also cannot handle the polarization equalization problem of probabilistic shaping systems under thunderstorm conditions. Summary of the Invention

[0004] To address the aforementioned shortcomings in the prior art, this invention provides a method and apparatus for joint equalization of PMD and RSOP damage, which solves the problem of joint equalization of PMD and RSOP damage in PCS-64QAM systems under thunderstorm scenarios.

[0005] To achieve the above objectives, the technical solution adopted by this invention is: a joint equalization method for PMD and RSOP damage, comprising the following steps:

[0006] S1. Based on the characteristics of polarization impairment in optical fiber links under thunderstorm scenarios, construct a joint equalization model for polarization impairment and initialize the parameters of the square root capacitive Kalman filter (SCKF).

[0007] S2. Based on the prediction equation, the posterior estimate of the state parameter at the previous time step is used as the prior estimate of the state parameter at the next time step. The square root capacitive Kalman filter (SCKF) is updated in time to predict the damage value of the damaged signal collected by the current sliding window.

[0008] S3. Based on the measurement equation of the sliding window, update the measurement of the square root capacitive Kalman filter (SCKF) to obtain the measurement prediction value.

[0009] S4. Combining the radius distribution of the received symbols, perform a ring decision operation based on the maximum a posteriori probability on the measurement prediction value to obtain the optimized radius of the constellation ring, calculate the innovation matrix, and combine the damage prediction value and the measurement prediction value to calculate the a posteriori estimate of the state parameter.

[0010] S5. Determine whether the iteration is complete. If yes, end the process and complete the joint equilibrium of PMD and RSOP damage. Otherwise, send the posterior estimate of the state parameter and the square root coefficient of the error covariance matrix to step S2, update the iteration number k, and proceed to the next iteration.

[0011] The beneficial effects of this invention are as follows: Based on a joint damage equalization model of RSOP and PMD and a square root capacitive Kalman filter, combined with a Rice distribution model, this invention uses the maximum a posteriori probability method for circular ring decision and calculates the optimized radius. Compared with the traditional EKF scheme, SCKF effectively reduces linearization error. Furthermore, this invention avoids circular ring decision errors caused by uneven constellation point distribution and ASE noise, and can especially solve the problem of joint equalization of polarization damage in thunderstorm scenarios. The introduction of probabilistic sensing innovation in this invention effectively solves the noise amplification problem caused by multiplicative innovation, greatly improving the computational accuracy of the PDMPCS-64QAM system in low optical signal-to-noise ratio (OSNR) scenarios. This invention exhibits excellent performance in low OSNR scenarios and has the advantages of high tracking accuracy, high tolerance to initial errors, and fast convergence speed.

[0012] Further, step S1 includes the following steps:

[0013] S101. Based on the characteristics of polarization damage in optical fiber links under thunderstorm scenarios, a joint equalization model for polarization damage is constructed. Based on the joint equalization model, the state parameters of the square root capacitive Kalman filter (SCKF) are selected and normalized.

[0014] S102. Initialize the square root commutative Kalman filter (SCKF) according to its settings. The expression for the measurement noise covariance R of the square root commutative Kalman filter (SCKF) is as follows:

[0015] R = diag([1; 1])

[0016] Where diag(·) represents a diagonal matrix operation;

[0017] The expression for the prediction noise covariance Q of the square root capacitive Kalman filter (SCKF) is as follows:

[0018] Q=diag([1E-4; 1E-4; 1E-4; 1E-6; 1E-6; 1E-6])

[0019] Where E is expressed in scientific notation, 1E-4 is 0.0001, and 1E-6 is 0.000001;

[0020] The error covariance matrix S of the square root capacitive Kalman filter SCKF 0|0 The expression is as follows:

[0021] S 0|0 =diag[(1E-2; 1E-2; 1E-2; 1E-2; 1E-2; 1E-2)].

[0022] The beneficial effects of the above-mentioned further solutions are as follows: Based on the correct establishment of the joint equalization model for polarization impairment, this invention can maximize the performance of SCKF; the normalization operation can ensure that the state parameters are of the same order of magnitude, reducing errors. Simultaneously, initializing the filter parameters can maximize the equalization accuracy and upper limit of the equalization performance of the solution.

[0023] Furthermore, step S101 includes the following steps:

[0024] S1011. Based on the characteristics of polarization impairment in fiber optic links under thunderstorm scenarios, RSOP is tracked in the time domain and PMD is compensated in the frequency domain to construct a joint equalization model for polarization impairment, where the RSOP tracking matrix is ​​R. eq The PMD compensation matrix is ​​U comp ;

[0025] The expression for tracking RSOP in the time domain is as follows:

[0026]

[0027] Among them, R eq Represents the equilibrium operator of the three-parameter RSOP, e (·) denoted by exponentiation, j represents the imaginary unit, ζ and η both represent the phase rotation angle of the RSOP damage to be tracked, and κ represents the azimuth rotation angle of the RSOP damage to be tracked.

[0028] The expression for compensating PMD in the frequency domain is as follows:

[0029]

[0030]

[0031]

[0032] Among them, U comp (ω) represents the compensation matrix for the first-order PMD, ω represents the optical angular frequency, Δτ represents the value of the differential group delay (DGD), and I represents the identity matrix. The vector representing PMD, Let denote the Pauli matrix, and τ1, τ2, τ3 represent the three components of the PMD vector in Stokes space. T Indicates the transpose operation;

[0033] S1012. Based on the joint equalization model, select and normalize the state parameters of the square root capacitive Kalman filter (SCKF):

[0034]

[0035] Among them, X norm T represents the state parameter of the square root capacitive Kalman filter (SCKF). s This indicates the symbol period of the current PCS-64QAM system.

[0036] The beneficial effects of the above-mentioned further scheme are: by performing equalization processing on RSOP in the time domain and compensating for PMD in the frequency domain, the present invention can obtain the six parameters tracked by the square root capacitive Kalman filter to achieve polarization equalization. At the same time, normalizing the state parameters can effectively reduce errors.

[0037] Furthermore, step S2 includes the following steps:

[0038] S201. Calculate the 2n column volume points X of the time update state of the square root volume Kalman filter (SCKF). i,k-1|k-1;

[0039] S202, Update the state of the 2n column volume points X over time. i,k-1|k-1 It propagates through the prediction equation;

[0040] S203, Predicted volume point values ​​after propagation using the prediction equation. By performing an averaging operation, the damage prediction value of the square root capacitive Kalman filter (SCKF) state parameter is obtained.

[0041] S204. Based on the damage prediction value of the damaged signal collected by the current sliding window. Calculate the square root factor S of the prediction error covariance matrix. k|k-1 The estimate is used to predict the damage value of the damaged signal collected in the current k-th sliding window.

[0042] Furthermore, the square root factor S k|k-1 The expression is as follows:

[0043]

[0044]

[0045]

[0046]

[0047]

[0048]

[0049] Where Tria(·) denotes the orthogonal triangular QR decomposition operation, Let S represent the center-weighted matrix. Q,k-1 Q represents the prediction noise covariance matrix. k-1 The square root factor, and All represent the predicted values ​​of the volume points after propagation through the prediction equation, where n represents the number of state parameters, i represents the sequence number of the volume points, i = 1, 2, ..., 2n, and S k-1|k-1 ξ represents the posterior estimate of the square root factor of the error covariance matrix at the previous time step. i Indicates the volume point. This represents the posterior estimate of the tracked state parameter at the previous time step. [1] i It represents the i-th column of point set 1.

[0050] The beneficial effect of the above-mentioned further scheme is that, during transmission, polarization impairments in the system are random, but their autocorrelation function characteristics indicate that polarization impairments have high correlation over a short period of time. Therefore, the establishment of the above process equations can better adapt to the constantly changing channel environment.

[0051] Furthermore, step S3 includes the following steps:

[0052] S301, based on the square root factor S k|k-1 Volume point ξ i and predicted value The volume point of the measurement update state of the square root capacitive Kalman filter (SCKF) is calculated using the following formula:

[0053]

[0054] S302, Based on the updated volume point X i,k|k-1 The k-th output sliding window signal sequence after polarization damage joint equalization is calculated.

[0055] S303. Based on the equalized k-th output sliding window signal sequence and the measurement equation of the fused sliding window, the measurement prediction value of the square root capacitive Kalman filter (SCKF) is calculated.

[0056]

[0057]

[0058] Where n represents the number of state parameters, i represents the sequence number of the volume point, i = 1, 2, ..., 2n, Z i,k|k-1 This represents the predicted volume point value, and abs(·) represents the modulo operation of a complex number. and Let represent the middle elements of the i-th group of X and Y polarization volume point sequences in the k-th sliding window, respectively. This indicates the rounding up operation.

[0059] The beneficial effect of the above scheme is that by performing a modulo operation on the sign of the intermediate sequence of the sliding window, the noise value can be reflected to the greatest extent.

[0060] Furthermore, step S302 includes the following steps:

[0061] S3021. Capture a segment of signal using a sliding window. Among them, S x,k With S y,k These represent the damaged X and Y polarization signal sequences extracted in the k-th sliding window, respectively.

[0062] S3022. Perform a Fourier transform on the intercepted signal to convert the signal to the frequency domain;

[0063] S3023. Based on the signal converted to the frequency domain, the volume points after the update state are propagated using the PMD compensation matrix to obtain the measurement prediction volume points after PMD compensation.

[0064] The expression for the PMD compensation matrix is ​​as follows:

[0065]

[0066]

[0067]

[0068] Among them, U comp,i (ω) represents the PMD compensation matrix, Δτ i U represents the PMD compensation matrix that uses the volume points updated from the i-th group of states as parameters. comp,i (ω) is the calculated differential group delay DGD. τ represents the vector of PMD obtained based on the i-th set of volume points as parameters. i,1,k|k-1 Let τ represent the value of the Stokes component τ1 in the i-th column volume point of the measurement update state. i,2,k|k-1 Let τ represent the value of the Stokes component τ2 at the i-th column volume point in the measurement update state. i,3,k|k-1 The value of the Stokes component τ3 in the i-th column volume point of the measurement update state;

[0069] S3024. Transform the measurement prediction volume points after PMD compensation to the time domain through inverse fast Fourier transform.

[0070] S3025. The measurement prediction volume points after PMD compensation converted to the time domain are propagated using the RSOP tracking matrix to obtain the measurement prediction volume points after PSOP tracking.

[0071] The expression for the RSOP tracking matrix is ​​as follows:

[0072]

[0073] Among them, R eq,i This indicates that the state parameters corresponding to the volume points in the i-th column are used as the tracking matrix of RSOP, ζ i,k|k-1 η represents the value of the phase rotation angle ζ in the i-th column volume point of the measurement update state. i,k|k-1 κ represents the value of the phase rotation angle η in the i-th column of the volume point in the measurement update state. i,k|k-1This represents the value of the azimuth rotation angle κ in the i-th column of the volume point in the measurement update state;

[0074] S3026. Average the measured predicted volume points after PSOP tracking to obtain the kth output sliding window signal sequence after polarization damage joint equalization.

[0075] The expression for the average calculation operation is as follows:

[0076]

[0077] in, and Let n and i represent the output sequences of the X and Y polarization sliding windows after the final joint equalization, respectively. Let n represent the number of state parameters and i represent the index of the volume point, i = 1, 2, ..., 2n. and They represent respectively through U comp,i and RSOP tracking matrix R eq,i The combined equalized X-polarization and Y-polarization volume point sequence, U comp,i This indicates the volume point X used to update the state using the measurement in column i. i,k|k-1 The first three parameters are used as the matrix of parameters for the PMD compensation matrix.

[0078] The beneficial effects of the above scheme are as follows: the measurement equation is propagated through volume points, and the state mean and covariance of the nonlinear system with additive white Gaussian noise are approximated by a set of volume points. It is theoretically the closest approximation algorithm to Bayesian filtering at present. Even in high-dimensional state space, the filter weights are always positive, and it has higher robustness and filtering accuracy than the unscented Kalman (UKF). RSOP occurs in the time domain and PMD occurs in the frequency domain. By using time-frequency conversion, RSOP is tracked in the time domain and PMD is compensated in the frequency domain, which conforms to the physical model constraints and can maximize the performance upper limit of polarization-damped joint equalization.

[0079] Furthermore, step S4 includes the following steps:

[0080] S401. Estimate the square root coefficient S of the new information covariance matrix. zz,k|k-1 :

[0081]

[0082]

[0083]

[0084] Where Tria(·) denotes the orthogonal triangular QR decomposition operation, Let S represent the weighted centrality matrix. R,k-1Represents the noise covariance matrix R k-1 The square root factor, n represents the number of state parameters, and i represents the index of the volume point, i = 1, 2, ..., 2n. Z represents the predicted measurement value. 2n,k|k-1 S represents the predicted volume point value. R,k-1 Represents the noise covariance matrix R k-1 The square root factor,

[0085] S402, Estimating the cross-covariance matrix P xz,k|k-1 :

[0086]

[0087]

[0088] Where, χ k|k-1 and Both represent center-weighted matrices, X 2n,k|k-1 The volume point represents the measurement update state. This represents the damage prediction value of the state parameter;

[0089] S403, based on the square root coefficient S zz,k|k-1 and cross-covariance matrix P xz,k|k-1 The Kalman gain W was calculated. k :

[0090]

[0091] S404. Combining the radius distribution of the received symbols, i.e., the Rice distribution model, to predict the measurement value. Perform a ring decision operation based on maximum a posteriori probability to obtain the optimized radius of the ring belonging to the constellation, and calculate the innovation matrix;

[0092] S405. Based on the innovation matrix, combined with damage prediction values Measurement prediction value and Kalman gain W k The posterior estimates of the state parameters are calculated.

[0093]

[0094] Among them, Z k This represents the target observation value of the k-th sliding window. This describes the process of calculating the innovation matrix.

[0095] S406. Update the mean root coefficient S of the error covariance matrix according to the following formula. k|k :

[0096]

[0097] Where, χ k|k-1 Represents the center-weighted matrix. Let S represent the center-weighted matrix. R,k-1 Represents the noise covariance matrix R k-1 The square root factor,

[0098] The beneficial effect of the above scheme is that it can estimate the cross-covariance matrix P by calculating it. xz,k|k-1 The square root coefficient S of the new information covariance matrix zz,k|k-1 Only then can the Kalman gain be calculated. The Kalman gain is an important parameter for measuring the degree of confidence in the observations and directly affects the posterior estimate of the state parameters. The mean root coefficient S of the error covariance matrix... k|k The calculations prepare for the next Kalman filter iteration.

[0099] Furthermore, the innovation matrix calculation process in step S404 includes the following steps:

[0100] S4041. Combining the Rice distribution model, the measurement prediction value is calculated based on the maximum a posteriori probability. Optimized radius of the constellation ring with the highest probability of belonging Among them, the predicted measurement value The process for determining which ring belongs to which ring is as follows:

[0101]

[0102]

[0103] in, This indicates the value of m when the function reaches its maximum value, where m represents each of the nine rings in the ideal PCS-64QAM constellation diagram. This indicates that the received signal amplitude is Under the condition that the amplitude of the transmitted signal is The probability is given by p(·), which represents the probability calculation operation. This indicates the radius of the constellation ring of the PCS-64QAM signal at the transmitting end. σ represents the predicted value of X or Y polarization in the k-th sliding window of the square root capacitive Kalman filter (SCKF). 2 I0(·) represents the noise variance, and I0(·) represents the 0th-order modified Bessel function. Represents different radii of the transmitter after probabilistic integer coding. The probability of occurrence;

[0104] S4042. Calculate the optimal radius using the average value of the Rice distribution. Among them, the optimized radius The expression is as follows:

[0105]

[0106] Where I1(x) represents the first-order modified Bessel function, exp(·) represents the exponential operation, and ln(·) represents the logarithmic operation;

[0107] S4043, optimize the radius Z of the annulus with the highest maximum posterior probability. k Compared with the predicted value By subtraction, we obtain the innovation matrix e:

[0108]

[0109]

[0110] in, This represents the optimized radius of the ring to which the predicted X-polarization values ​​of the PCS-64QAM system belong. This represents the optimized radius of the ring to which the predicted Y-polarization values ​​of the PCS-64QAM system belong. The measured predicted value representing X-polarization is... The measured predicted value representing Y-polarization is... L represents the length of the sliding window, S x,k With S y,k These represent the damaged X or Y polarization signal sequences captured by the k-th sliding window, respectively.

[0111] The beneficial effects of the above scheme are as follows: In PCS systems, due to the uneven distribution of constellation points and the influence of amplifier spontaneous emission noise (ASE), the traditional Euclidean distance-based decision method is suboptimal. By introducing maximum a posteriori probability estimation, the decision error caused by the traditional Euclidean distance decision threshold can be effectively solved, the probability of circular decision errors in low OSNR scenarios can be effectively reduced, and the filtering accuracy in low OSNR scenarios can be improved.

[0112] This invention provides a joint equalization device for a joint equalization method of PMD and RSOP damage, comprising:

[0113] The first processing module constructs a joint equalization model for polarization damage based on the characteristics of polarization damage in fiber optic links during thunderstorms, and initializes the parameters of the square root commutative Kalman filter (SCKF). The second processing module updates the SCKF by using the posterior estimates of the state parameters from the previous time step as prior estimates for the state parameters at the next time step, based on the prediction equation, to predict the damage value of the damaged signal acquired in the current sliding window. The third processing module updates the SCKF by measurement based on the measurement equation of the fused sliding window. The first processing module obtains the measurement prediction value; the second processing module, in conjunction with the radius distribution of the received symbols, performs a ring decision operation based on the maximum a posteriori probability on the measurement prediction value to obtain the optimized radius of the ring of the constellation to which it belongs, calculates the innovation matrix, and calculates the posterior estimate of the state parameter by combining the damage prediction value and the measurement prediction value; the third processing module is used to determine whether the iteration is complete. If so, the process ends and the joint equilibrium of PMD and RSOP damage is completed. Otherwise, the posterior estimate of the state parameter and the square root coefficient of the error covariance matrix are sent to the second processing module, and the iteration number k is updated for the next iteration.

[0114] The beneficial effects of this invention are as follows: Based on the joint damage equalization model of RSOP and PMD and SCKF, this invention integrates a probabilistic sensing model with maximum a posteriori probability, effectively reducing linearization error compared to the traditional EKF scheme. Furthermore, this invention avoids circular decision errors caused by uneven constellation point distribution and ASE noise, and particularly solves the problem of joint equalization of polarization damage in thunderstorm scenarios. The introduction of probabilistic sensing innovation in this invention effectively solves the noise amplification problem caused by multiplicative innovation, greatly improving the computational accuracy of the PDM PCS-64QAM system in low OSNR scenarios. This invention exhibits excellent performance in low OSNR scenarios and has the advantages of high tracking accuracy, high tolerance to initial errors, and fast convergence speed. Attached Figure Description

[0115] Figure 1 This is a flowchart of the method of the present invention.

[0116] Figure 2 This is a schematic diagram of the simulation system of the present invention.

[0117] Figure 3 This diagram illustrates the upper limit of RSOP performance under different schemes and different differential group delays (DGD) for OSNR=25dB in this embodiment.

[0118] Figure 4 This is a schematic diagram of the relationship between OSNR and NGMI under different polarization damage conditions in this embodiment.

[0119] Figure 5This is a schematic diagram showing the relationship between source entropy H and normalized generalized mutual information (NGMI) under different polarization impairments in this embodiment.

[0120] Figure 6 This is a schematic diagram of the device structure of the present invention. Detailed Implementation

[0121] The specific embodiments of the present invention are described below to enable those skilled in the art to understand the present invention. However, it should be understood that the present invention is not limited to the scope of the specific embodiments. For those skilled in the art, various changes are obvious as long as they are within the spirit and scope of the present invention as defined and determined by the appended claims. All inventions utilizing the concept of the present invention are protected.

[0122] Example 1

[0123] Before describing this invention, the parameters k and i are explained as follows:

[0124] The kth iteration is the kth sliding window;

[0125] The volume points of the state parameters in the i-th column are used as parameters of the RSOP and PMD compensation matrices, and the compensated sequence is the volume points of the i-th group of signal sequences.

[0126] To address the aforementioned shortcomings in existing technologies, this invention provides a joint equalization method for RSOP and PMD damage in a PDM-PCS-64QAM system under thunderstorm conditions. This method can jointly equalize RSOP, which can reach several Mrad / s, and PMD, which can last for several symbol periods, under thunderstorm conditions. Figure 1 As shown, a joint equalization method for PMD and RSOP damage is implemented as follows:

[0127] S1. Based on the characteristics of polarization impairment in optical fiber links under thunderstorm scenarios, a joint equalization model for polarization impairment is constructed, and the parameters of the square root capacitive Kalman filter (SCKF) are initialized. The implementation method is as follows:

[0128] S101. Based on the characteristics of polarization impairment in optical fiber links under thunderstorm scenarios, a joint equalization model for polarization impairment is constructed. Then, based on the joint equalization model, the state parameters of the square root capacitive Kalman filter (SCKF) are selected and normalized. The implementation method is as follows:

[0129] S1011. Based on the characteristics of polarization damage in fiber optic links under thunderstorm scenarios, RSOP is tracked in the time domain and PMD is compensated in the frequency domain to construct a joint equalization model for polarization damage; S1012. Based on the joint equalization model, the state parameters of the square root capacitive Kalman filter SCKF are selected and normalized.

[0130] S102. Initialize the square root capacitive Kalman filter (SCKF) according to its settings.

[0131] S2. Based on the prediction equation, the posterior estimate of the state parameter at the previous time step is used as the prior estimate of the state parameter at the next time step. The square root capillary Kalman filter (SCKF) is updated over time to predict the damage value of the damaged signal acquired in the current sliding window. The implementation method is as follows:

[0132] S201. Calculate the 2n column volume points X of the time update state of the square root volume Kalman filter (SCKF). i,k-1|k-1 S202, Update the volume points X of the 2n columns of the state over time. i,k-1|k-1 Propagation is performed through the prediction equation; S203, the volume point prediction values ​​after propagation using the prediction equation are... By performing an averaging operation, the damage prediction value of the square root capacitive Kalman filter (SCKF) state parameter is obtained. S204. Based on the damage prediction value of the damaged signal collected by the current sliding window. Calculate the square root factor S of the prediction error covariance matrix. k|k-1 The estimate is used to predict the damage value of the damaged signal collected in the current k-th sliding window;

[0133] S3. Based on the measurement equation of the fused sliding window, the square root capacitive Kalman filter (SCKF) is updated to obtain the measurement prediction value. The implementation method is as follows:

[0134] S301, based on the square root factor S k|k-1 Volume point ξ i and predicted value The volume point of the measurement update state of the square root capacitive Kalman filter (SCKF) is calculated using the following formula:

[0135]

[0136] S302, Based on the updated volume point X i,k|k-1 The k-th output sliding window signal sequence after polarization-damped joint equalization is calculated, and its implementation method is as follows:

[0137] S3021. Capture a segment of signal using a sliding window. Among them, S x,k With S y,kS3022: Perform a Fourier transform on the intercepted signal to convert it to the frequency domain; S3023: Based on the signal converted to the frequency domain, propagate the updated volume points using the PMD compensation matrix to obtain the PMD-compensated measurement prediction volume points; S3024: Convert the PMD-compensated measurement prediction volume points to the time domain using an inverse fast Fourier transform; S3025: Propagate the PMD-compensated measurement prediction volume points converted to the time domain using the RSOP tracking matrix to obtain the PSOP-tracked measurement prediction volume points; S3026: Average the PSOP-tracked measurement prediction volume points to obtain the k-th output sliding window signal sequence after polarization damage joint equalization.

[0138] S303. Based on the equalized k-th output sliding window signal sequence and the measurement equation of the fused sliding window, the measurement prediction value of the square root capacitive Kalman filter (SCKF) is calculated.

[0139] S4. Combining the radius distribution of the received symbols, perform a ring decision operation based on the maximum a posteriori probability on the measurement prediction values ​​to obtain the optimized radius of the constellation ring, calculate the innovation matrix, and combine the damage prediction values ​​and the measurement prediction values ​​to calculate the posterior estimates of the state parameters. The implementation method is as follows:

[0140] S401. Estimate the square root coefficient S of the new information covariance matrix. zz,k|k-1 ;

[0141] S402, Estimating the cross-covariance matrix P xz,k|k-1 ;

[0142] S403, based on the square root coefficient S zz,k|k-1 and cross-covariance matrix P xz,k|k-1 The Kalman gain W was calculated. k ;

[0143] S404, Predicted Measurement Values Perform a ring decision operation based on the maximum a posteriori probability to obtain the standard radius of the ring, and calculate the innovation matrix. The implementation method is as follows:

[0144] S4041. Combining the Rice distribution model, the measurement prediction value is calculated based on the maximum a posteriori probability. Optimized radius of the constellation ring with the highest probability of belonging S4042. Calculate the optimized radius using the average value of the Rice distribution. S4043, optimize the radius Z of the annulus with the highest maximum posterior probability. k Compared with the predicted value By subtracting the components, we obtain the innovation matrix e;

[0145] S405. Based on the innovation matrix, combined with damage prediction values Measurement prediction value and Kalman gain W k The posterior estimates of the state parameters are calculated.

[0146] S406, Update the root coefficients of the mean of the error covariance matrix S k|k ;

[0147] S5. Determine whether the iteration is complete. If yes, end the process and complete the joint equilibrium of PMD and RSOP damage. Otherwise, send the posterior estimate of the state parameter and the square root coefficient of the error covariance matrix to step S2, update k, and proceed to the next iteration, where k represents the kth iteration.

[0148] In this embodiment, a joint equalization model is constructed, and state parameters are selected: Based on the damage characteristics, RSOP is equalized in the time domain using formula (1), and PMD is compensated in the frequency domain using formula (2).

[0149]

[0150]

[0151] Among them, R eq The equilibrium quantum of the three-parameter PSOP is represented by e, which represents the exponentiation operation, j, which represents the imaginary number, ζ and η both represent the phase rotation angle, κ represents the azimuth rotation angle, and U comp (ω) represents the compensation matrix of the first-order PMD, ω represents the optical angular frequency, and the PMD vector is described as... The superscript T indicates the transpose operation, τ1, τ2, τ3 represent the three components in the Stokes space, and the DGD value... I represents the identity matrix. This represents the Pauli matrix.

[0152] Therefore, τ1, τ2, τ3, κ, η, ζ become the six parameters tracked by the square root capacitive Kalman filter to achieve polarization equalization. To reduce errors and normalize the state parameters, the final state parameters of the SCKF are expressed as:

[0153]

[0154] Among them, X norm T represents the state parameter of the square root capacitive Kalman filter (SCKF). s This indicates the symbol period of the current PCS-64QAM system.

[0155] In this embodiment, the correct establishment of the joint equalization model for polarization damage can maximize the performance of SCKF, and the normalization operation can ensure that the state parameters are at the same order of magnitude, thus reducing errors.

[0156] In this embodiment, the square root commutative Kalman filter is initialized according to the settings of the square root commutative Kalman filter.

[0157] Measurement noise covariance R:

[0158] R = diag([1; 1]);

[0159] Here, diag(·) represents a diagonal matrix operation.

[0160] Predicted noise covariance Q:

[0161] Q=diag([1E-4; 1E-4; 1E-4; 1E-6; 1E-6; 1E-6]);

[0162] Where E is expressed in scientific notation, 1E-4 is 0.0001, and 1E-6 is 0.000001;

[0163] Error covariance matrix S 0|0 :

[0164] S 0|0 =diag[(1E-2; 1E-2; 1E-2; 1E-2; 1E-2; 1E-2)].

[0165] In this embodiment, the filter parameters are initialized to maximize the equalization accuracy and upper limit of the equalization performance of the scheme.

[0166] In this embodiment, the 2n column volume points X of the time update state are calculated according to formula (4). i,k-1|k-1 .

[0167]

[0168] in, [1] i Let S represent the i-th column of point set 1 (i = 1, 2, 3, ..., n), where n is the number of state parameters, and the subscript k represents the iteration number index of the filter, i.e., the number of sliding windows. k-1|k-1 This represents the posterior estimate of the square root factor of the error covariance matrix at the previous time step. This represents the posterior estimate of the tracked state parameter at the previous moment. It should be noted that the number of columns of the volume points is 2n, which is determined according to the radial volume rule of the third-order sphere.

[0169] In this embodiment, the 2n column volume points of the time-updated state are propagated through a prediction equation. Regarding the construction process of the prediction equation in this invention, due to the high randomness of polarization damage, these damages cannot be accurately predicted. Therefore, the predicted value of the current k-th iteration is assumed to be the posterior estimate of the (k-1)-th iteration. That is, the propagation process of the time-updated state volume points is calculated according to the following formula:

[0170]

[0171] in, X represents the volumetric point prediction value after propagation through the prediction equation. k-1|k-1 This represents the posterior estimate of the state parameters at the previous time step, where n is the number of state parameters and i represents the index of the volume point.

[0172] In this embodiment, the predicted value of the state parameter is calculated. Volume point prediction value after propagation through the prediction equation By averaging the values, we can obtain the predicted values ​​of the state parameters. The calculation formula is:

[0173]

[0174] In this embodiment, the square root factor S of the prediction error covariance matrix is ​​estimated according to the following formula. k|k-1 :

[0175]

[0176] Where Tria(·) denotes the orthogonal triangular QR decomposition operation, Let S represent the center-weighted matrix. Q,k-1 Q represents the prediction noise covariance matrix. k-1 The square root factor, T represents the transpose operation. and All represent the predicted values ​​of the volume points after propagation through the prediction equation, where n represents the number of state parameters, i represents the sequence number of the volume points, i = 1, 2, ..., 2n, and S k-1|k-1 ξ represents the posterior estimate of the square root factor of the error covariance matrix at the previous time step. i Indicates the volume point. This represents the posterior estimate of the tracked state parameter at the previous time step. [1] i It represents the i-th column of point set 1.

[0177] In this embodiment, polarization impairments in the system are random during transmission, but the autocorrelation function characteristics indicate that polarization impairments have a high correlation in a short period of time. Therefore, the establishment of the above process equations can better adapt to the ever-changing channel environment.

[0178] In this embodiment, based on the square root factor S obtained in step S2 k|k-1 Volume point ξ i and predicted value The volume points (i = 1, 2, ..., 2n) of the measurement update state are calculated using the following formula:

[0179]

[0180] In this embodiment, the signal sequence after joint equalization for polarization impairment is calculated according to the measurement equation:

[0181] Extracting a segment of the signal using a sliding window:

[0182]

[0183] Among them, S x,k With S y,k Let X and Y represent the damaged X and Y polarization signal sequences intercepted by the k-th sliding window, respectively. Performing a Fourier transform on the intercepted signals converts them to the frequency domain, yielding:

[0184]

[0185] Where FFT{·} denotes Fast Fourier Transform.

[0186] The signal converted to the frequency domain is propagated using formula (11) of PMD compensation matrix for volume point propagation:

[0187]

[0188] Among them, U comp,i (ω) represents the PMD compensation matrix, Δτ i U represents the PMD compensation matrix that uses the volume points updated from the i-th group of states as parameters. comp,i (ω) is the calculated differential group delay DGD. τ represents the vector of PMD obtained based on the i-th set of volume points as parameters. i,1,k|k-1 Let τ represent the value of the Stokes component τ1 at the volume point corresponding to the measurement update state in the i-th column. i,2,k|k-1 Let τ represent the value of the Stokes component τ2 at the volume point corresponding to the measurement update state in the i-th column. i,3,k|k-1 The value of the Stokes component τ3 in the volume point of the measurement update state in the i-th column;

[0189] It should be noted that in actual operation, the propagation process of the volume point requires attention to the position of the state parameters. In formula (3), the position of the tracked state parameters has been specified. Taking the volume point with i=1, i.e., the state update of the first column, as an example, the corresponding state parameter order is: Where, τ 1,1,k|k-1 This represents the value of the Stokes component τ1 at the volume point of the measurement update state in column 1. 2,1,k|k-1 This represents the value of the Stokes component τ2 at the volume point of the measurement update state in the first column, and so on. Therefore, U comp,i It uses the volume point X of the i-th column to measure and update the state. i,k|k-1 The first three parameters are used as the matrix of parameters for the PMD compensation matrix, namely:

[0190]

[0191] Since there are 12 sets of volume points, 12 PMD compensation matrices U can be obtained. comp,i After compensating for the damaged signals, 12 sets of PMD-compensated measurement prediction volume points were obtained:

[0192]

[0193] Then, the 12 sets of PMD-compensated measurement prediction volume points were transformed to the time domain using inverse fast Fourier transform (IFFT):

[0194]

[0195] The measurement prediction volume point after PMD compensation (i.e., formula (13)) converted to the time domain is then propagated according to the RSOP tracking matrix (14):

[0196]

[0197] Among them, R eq,i This indicates that the state parameters corresponding to the volume points in the i-th column are used as the tracking matrix of RSOP, ζ i,kk-1 η represents the value of the phase rotation angle ζ at the volume point corresponding to the measurement update state in the i-th column. i,kk-1 κ represents the value of the phase rotation angle η at the volume point corresponding to the measurement update state in the i-th column. i,k|k-1 This represents the value of the azimuth rotation angle κ in the volume point corresponding to the measurement update state in the i-th column;

[0198] The parameters of the RSOP tracking matrix are determined based on the volume points of the measurement update state. Since there are 12 sets of volume points, there are 12 RSOP tracking matrices R. eq,iRSOP tracking was performed on the 12 sets of time-domain measurement prediction volume points after PMD compensation obtained by formula (13) to obtain the 12 sets of measurement prediction volume points after RSOP tracking:

[0199]

[0200] For ease of description, the sequence after joint equilibrium of formula (15) is written as and They represent respectively through U comp,i and R eq,i The X-polarized and Y-polarized volume point sequence after joint equalization.

[0201] The volume point sequence obtained from formula (15) There are 12 groups. To obtain the final symbol sequence after joint equalization of polarization damage, the volume point sequence needs to be averaged using the following formula:

[0202]

[0203] in, and Let n and i represent the output sequences of the X and Y polarization sliding windows after the final joint equalization, respectively. Let n represent the number of state parameters and i represent the index of the volume point, i = 1, 2, ..., 2n. and They represent respectively through U comp,i and RSOP tracking matrix R eq,i The combined equalized X-polarization and Y-polarization volume point sequence, U comp,i This indicates the volume point X used to update the state using the measurement in column i. i,k|k-1 The first three parameters are used as the matrix of parameters for the PMD compensation matrix.

[0204] Finally, the predicted measurement value is calculated based on the measurement equation.

[0205] In this embodiment, the volume point sequence obtained according to formula (15) There are 12 groups. The predicted volume point value Z is calculated according to formula (17). i,k|k-1 We take the middle element of each volume point sequence and perform a modulo operation. Since the sequence length may be odd or even, we round it up to prevent decimals.

[0206]

[0207] The 12 sets of measured predicted volume point values ​​Z obtained according to formula (17) i,k|k-1 The average value is calculated to obtain the measurement prediction value.

[0208]

[0209] Among them, Z i,k|k-1 This represents the predicted volume point value, and abs(·) represents the modulo operation of a complex number. and Let represent the middle elements of the i-th group of X and Y polarization volume point sequences in the k-th sliding window, respectively. This indicates the rounding up operation.

[0210] In this embodiment, the square root coefficient S of the new information covariance matrix zz,k|k-1 Calculate according to formula (19):

[0211]

[0212]

[0213]

[0214] Where Tria(·) denotes the orthogonal triangular QR decomposition operation, Let S represent the weighted centrality matrix. R,k-1 Represents the noise covariance matrix R k-1 The square root factor, n represents the number of state parameters, and i represents the index of the volume point, i = 1, 2, ..., 2n. Z represents the predicted measurement value. 2n,k|k-1 S represents the predicted volume point value. R,k-1 Represents the noise covariance matrix R k-1 The square root factor,

[0215] The cross-covariance matrix P is calculated according to formula (20). xz,k|k-1 Estimate:

[0216]

[0217] in, χ k|k-1 and Both represent center-weighted matrices, X 2n,k|k-1 The volume point represents the measurement update state. This represents the damage prediction value of the state parameter, and the superscript T indicates the transpose operation.

[0218] Calculate the Kalman gain W k The calculation of Kalman gain depends on the cross-covariance matrix P. xz,k|k-1 And the square root coefficient S of the new information covariance matrix zz,k|k-1The calculation formula is shown in formula (21):

[0219]

[0220] Calculate the posterior estimate This process involves fusing the prediction results and measurement results to obtain a more accurate posterior estimate of the state parameters. The calculation formula is as follows:

[0221]

[0222] Among them, Z k This represents the target observation value of the k-th sliding window. The innovation calculation process is an important part of the SCKF design, as it can measure the reliability of measurement results.

[0223] In this embodiment, the innovation calculation method based on maximum a posteriori probability sensing is as follows: Because the measurement prediction value is affected by probabilistic constellation shaping and amplifier spontaneous emission noise (ASE), the traditional Euclidean distance decision scheme is suboptimal for probabilistic constellation shaping systems. To adapt to the PDM PCS-64QAM system, this invention proposes a decision scheme based on maximum a posteriori probability sensing, which can obtain more accurate target observation values ​​Z. k To assist in information calculation. Due to the presence of frequency offset phase noise in the system, the constellation points of the PDM PCS-64QAM system will exhibit multiple rings after polarization equalization. This invention first combines the radius distribution of the received symbols, i.e., the Ricean distribution model, to predict the measurement values. Perform a circular decision operation based on the maximum a posteriori probability, and then optimize the radius Z with the highest probability. k Compared with the predicted value Divide the difference and calculate the new information e.

[0224] In this embodiment, the specific determination process for the constellation ring to which the measured predicted value belongs can be represented as follows:

[0225]

[0226]

[0227] in, This indicates the value of m when the function reaches its maximum value, where m represents each of the nine rings in the ideal PCS-64QAM constellation diagram. This indicates that the received signal amplitude is Under the condition that the amplitude of the transmitted signal is The probability is given by p(·), which represents the probability calculation operation. This indicates the radius of the constellation ring of the PCS-64QAM signal at the transmitting end. σ represents the predicted value of X or Y polarization in the k-th sliding window of the square root capacitive Kalman filter (SCKF). 2 Represents the noise variance, and I0(·) is the 0th-order modified Bessel function. It is a transmitter with different radii via probabilistic integer coding. The probability of occurrence;

[0228] Optimized radius The value is calculated using the average of the Rice distribution. An optimized radius is employed. Z, the target observation in SCKF k To calculate more accurate information. Optimize radius. The calculation formula is:

[0229]

[0230] Where I1(x) represents the first-order modified Bessel function, exp(·) represents the exponential operation, and ln(·) represents the logarithmic operation.

[0231] At this point, the resulting innovation matrix e can be written as:

[0232]

[0233]

[0234] in, This represents the optimized radius of the ring to which the predicted X-polarization values ​​of the PCS-64QAM system belong. This represents the optimized radius of the ring to which the predicted Y-polarization values ​​of the PCS-64QAM system belong. The measured predicted value representing X-polarization is... The measured predicted value representing Y-polarization is... L represents the length of the sliding window, S x,k With S y,k These represent the damaged X or Y polarization signal sequences captured by the k-th sliding window, respectively. In this embodiment, the square root coefficient S of the error covariance matrix is ​​updated according to formula (27). k|k .

[0235] S k|k =Tria([χ k|k-1 -W k z k|k-1 W k S R,k (27)

[0236] In this embodiment, the sliding window slides forward with a step size Δl = 1, and the signal of the next window is subjected to damage equalization. Steps S2, S3, and S4 are repeated until the polarization damage of all signals is jointly equalized.

[0237] In this embodiment, the present invention uses a recursive iteration method to continuously move the sliding window, which can effectively reduce the filter delay. The iteration operation is repeated until the polarization impairment of all signals is equalized.

[0238] In this embodiment, to verify the effectiveness of the present invention, the feasibility of the proposed scheme was verified using the optical communication simulation software VPItransmission Makers 11.2 and MATLAB. A schematic diagram of the simulation system is shown below. Figure 2 As shown. To obtain the PCS-64QAM signal, the transmitter first generates the required bit sequence, then generates the required probability distribution through CCDM, and finally generates a 28Gbaud PDM-PCS-64QAM constellation symbol using a symbol mapping step. The transmitter uses a square root raised cosine filter with a roll-off factor of 0.1 for pulse shaping. The laser center wavelength is 1550nm. The signal output by the transmitter is first superimposed with ASE noise, and then passes through a three-stage polarization impairment simulator consisting of RSOP1, PMD, and RSOP2. In the simulator, the rate of change of PSPs is the same as the rotation speed of RSOP, and the azimuth and phase angles of RSOP increase linearly with time. Before entering the receiver, it passes through a 10th-order Bessel filter with a bandwidth of 33.6GHz for noise suppression. At the receiver, the linewidth and frequency offset of the local oscillator laser are set to 100kHz and 100MHz, respectively. After entering the coherent receiver, the DSP process first performs clock recovery, then the signal is resampled to 2 samples / symbol, and normalized and standardized. Next, polarization equalization is performed. After polarization equalization, the signal is downsampled to 1 symbol / sample, followed by frequency offset estimation for frequency offset recovery, and a blind phase search (BPS) scheme is used for carrier phase recovery. Finally, the normalized generalized mutual information (NGMI) of the signal is calculated. SCKF-PA will be used as an abbreviation for this invention. In contrast, the EKF scheme using multiplicative innovation is abbreviated as EKF-CM, and the traditional multimode algorithm is abbreviated as MMA.

[0239] In this embodiment, to test the performance ceiling of the proposed scheme, the system performance ceiling was tested under the conditions of 25dB OSNR and source entropy H = 4.5bit / symbol, respectively. Figure 3As shown, it can be observed that when DGD gradually increases from 30 ps to 90 ps, ​​the MMA and EKF-CM schemes used in comparison are almost incapable of handling the combined polarization damage under any RSOP. However, the SCKF-PA scheme of this invention exhibits superior performance. Under strong shaping conditions of H = 4.5 bits / symbol and an NGMI threshold of 0.9, this invention can jointly balance the maximum RSOP of 8.34 Mrad / s and the DGD damage of 90 ps.

[0240] In this embodiment, the present invention verifies the OSNR cost under the condition of source entropy of 4.5 bits / symbol, such as... Figure 4 As shown, when the OSNR test range increases from 17dB to 27dB in 1dB increments, joint damage recovery is performed under scenarios of DGD=35ps, RSOP=3Mrad / s and DGD=70ps, RSOP=6Mrad / s. Simulation results show that for a damage-free curve, approximately 18dB OSNR is required to reach the NGMI threshold. The EKF-CM and MMA schemes used for comparison cannot balance such high polarization damage. Under the conditions of DGD=35ps, RSOP=3Mrad / s, compared to the damage-free curve, this invention achieves an NGMI threshold of 0.9 with only 0.79dB OSNR cost. Even in the more extreme scenario of DGD=70ps, RSOP=6Mrad / s, the OSNR cost of this invention is only 0.9dB. Therefore, this invention has a lower OSNR cost and superior performance.

[0241] In this embodiment, to illustrate the applicability and flexibility of the SCKF-PA scheme for the PCS-64QAM system, Figure 5 The NGMI curves for three schemes under different shaping intensities and different polarization damage conditions at OSNR=23dB are presented. Figure 5 The horizontal axis represents the source entropy H, ranging from 4.1 bit / symbol to 5.9 bit / symbol in intervals of 0.1 bit / symbol. The vertical axis represents NGMI. It can be seen that, without any polarization impairment, when the OSNR is 23 dB, due to the influence of ASE noise and probabilistic constellation shaping, as the source entropy increases from 4.1 bit / symbol to 5.69 bit / symbol, the NGMI value of the PCS-64QAM system remains above the threshold of 0.9. As a comparison, the MMA and EKF-CM schemes... Figure 5All of the conditions listed cannot be balanced. However, when the channel DGD = 35ps and RSOP = 3Mrad / s, this invention can achieve an NGMI threshold of 0.9 or higher within a source entropy range of 4.1 bit / symbol to 5.56 bit / symbol. When the impairment further increases to DGD = 70ps and RSOP = 6Mrad / s, the maximum source entropy range that this invention can handle is 4.1 bit / symbol to 5.3 bit / symbol.

[0242] Example 2

[0243] like Figure 6 As shown, the present invention provides a joint equalization device for a joint equalization method of PMD and RSOP damage, comprising:

[0244] The first processing module constructs a joint equalization model for polarization damage based on the characteristics of polarization damage in fiber optic links during thunderstorms, and initializes the parameters of the square root commutative Kalman filter (SCKF). The second processing module updates the SCKF by using the posterior estimates of the state parameters from the previous time step as prior estimates for the state parameters at the next time step, based on the prediction equation, to predict the damage value of the damaged signal acquired in the current sliding window. The third processing module updates the SCKF by measurement based on the measurement equation of the fused sliding window. The first processing module obtains the measurement prediction value; the second processing module, in conjunction with the radius distribution of the received symbols, performs a ring decision operation based on the maximum a posteriori probability on the measurement prediction value to obtain the optimized radius of the ring of the constellation to which it belongs, calculates the innovation matrix, and calculates the posterior estimate of the state parameter by combining the damage prediction value and the measurement prediction value; the third processing module is used to determine whether the iteration is complete. If so, the process ends and the joint equilibrium of PMD and RSOP damage is completed. Otherwise, the posterior estimate of the state parameter and the square root coefficient of the error covariance matrix are sent to the second processing module, and the iteration number k is updated for the next iteration.

[0245] like Figure 6 The joint equalization registration device provided in the embodiment shown can execute the technical solution shown in the joint equalization registration method of the above method embodiment. Its implementation principle and beneficial effects are similar, and will not be repeated here.

[0246] In this embodiment, the functional units can be divided according to the joint equalization registration method. For example, each function can be divided into its own functional units, or two or more functions can be integrated into one processing unit. The integrated unit can be implemented in hardware or as a software functional unit. It should be noted that the unit division in this invention is illustrative and represents only a logical division; in actual implementation, other division methods may be used.

[0247] In this embodiment of the invention, the joint equalization registration device, in order to realize the principle and beneficial effects of the joint equalization registration method, includes hardware structures and / or software modules corresponding to each function. Those skilled in the art should readily recognize that, in conjunction with the illustrative units and algorithm steps described in the embodiments disclosed herein, the invention can be implemented in hardware and / or a combination of hardware and computer software. Whether a function is executed by hardware or computer software depends on the specific application and design constraints of the technical solution. Different methods can be used to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0248] In this embodiment, the present invention is based on a joint damage equalization model of RSOP and PMD and SCKF, fusing a probabilistic sensing model with maximum a posteriori probability, which effectively reduces linearization error compared to the traditional EKF scheme. Furthermore, the present invention can avoid circular decision errors caused by uneven distribution of constellation points and ASE noise, and can particularly solve the problem of joint equalization of polarization damage in thunderstorm scenarios. The introduction of probabilistic sensing innovation in the present invention effectively solves the noise amplification problem caused by multiplicative innovation, greatly improving the computational accuracy of the PDM PCS-64QAM system in low OSNR scenarios. The present invention exhibits excellent performance in low OSNR scenarios and has the advantages of high tracking accuracy, high tolerance to initial errors, and fast convergence speed.

Claims

1. A joint equalization method for PMD and RSOP damage, characterized in that, Includes the following steps: S1. Based on the characteristics of polarization damage in the optical fiber link of the PCS-64QAM system under thunderstorm scenarios, construct a joint equalization model for polarization damage and initialize the parameters of the square root capacitive Kalman filter (SCKF). S2. Based on the prediction equation, the posterior estimate of the state parameter at the previous time step is used as the prior estimate of the state parameter at the next time step. The square root capacitive Kalman filter (SCKF) is updated in time to predict the damage value of the damaged signal collected by the current sliding window. S3. Based on the measurement equation of the sliding window, update the measurement of the square root capacitive Kalman filter (SCKF) to obtain the measurement prediction value. S4. Combining the radius distribution of the received symbols, perform a ring decision operation based on the maximum a posteriori probability on the measurement prediction value to obtain the optimized radius of the constellation ring, calculate the innovation matrix, and combine the damage prediction value and the measurement prediction value to calculate the a posteriori estimate of the state parameter. S5. Determine if the iteration is complete. If yes, end the process and complete the joint equilibrium of PMD and RSOP damages. Otherwise, send the posterior estimates of the state parameters and the square root coefficients of the error covariance matrix to step S2, and calculate the iteration count. k Update and proceed with the next iteration.

2. The joint equalization method for PMD and RSOP damage according to claim 1, characterized in that, Step S1 includes the following steps: S101. Based on the characteristics of polarization damage in optical fiber links under thunderstorm scenarios, a joint equalization model for polarization damage is constructed. Based on the joint equalization model, the state parameters of the square root capacitive Kalman filter (SCKF) are selected and normalized. S102. Initialize the square root capacitive Kalman filter (SCKF) according to its settings, where the measurement noise covariance of the SCKF is... The expression is as follows: in, Indicates diagonal matrix operations; Predictive noise covariance of the square root capacitive Kalman filter (SCKF) The expression is as follows: in, In scientific notation, It is 0.0001. It is 0.000001; Error covariance matrix of square root capacitive Kalman filter (SCKF) The expression is as follows: 。 3. The method for jointly balancing PMD and RSOP damage according to claim 2, characterized in that, Step S101 includes the following steps: S1011. Based on the characteristics of polarization impairment in fiber optic links under thunderstorm scenarios, RSOP is tracked in the time domain and PMD is compensated in the frequency domain to construct a joint equalization model for polarization impairment. The RSOP tracking matrix is: The PMD compensation matrix is ; The expression for tracking RSOP in the time domain is as follows: in, This represents the tracking matrix of the three-parameter RSOP. This indicates exponentiation. Represents the imaginary unit. and All represent the phase rotation angle of the RSOP damage to be tracked. Indicates the azimuth rotation angle of the RSOP damage to be tracked; The expression for compensating PMD in the frequency domain is as follows: in, This represents the compensation matrix for first-order PMD. Indicates the optical angular frequency. This represents the value of the differential group delay (DGD). Represents the identity matrix. The vector representing PMD, Represents the Pauli matrix. The three components of the PMD vector in Stokes space are represented. Indicates the transpose operation; S1012. Based on the joint equalization model, select and normalize the state parameters of the square root capacitive Kalman filter (SCKF): in, This represents the state parameters of the square root capacitive Kalman filter (SCKF). This indicates the symbol period of the current PCS-64QAM system.

4. The joint equalization method for PMD and RSOP damage according to claim 1, characterized in that, Step S2 includes the following steps: S201, Calculate the time update state of the square root volume Kalman filter (SCKF) 2. n Column volume points ; S202, Update status 2 for time. n Column volume points It propagates through the prediction equation; S203. Predicted volume point values ​​after propagation using the prediction equation. By performing an averaging operation, the damage prediction value of the square root capacitive Kalman filter (SCKF) state parameter is obtained. ; S204. Based on the damage prediction value of the damaged signal collected by the current sliding window. To calculate the square root factor of the prediction error covariance matrix. The estimate, in order to assess the current number k Damage values ​​are predicted by collecting damaged signals from a sliding window.

5. The joint equalization method for PMD and RSOP damage according to claim 4, characterized in that, The square root factor The expression is as follows: in, This represents the orthogonal triangular QR decomposition operation. Represents the center-weighted matrix. Represents the prediction noise covariance matrix The square root factor, , and All represent the volumetric point prediction values ​​after propagation through the prediction equation. Indicates the number of state parameters. The sequence number representing the volume point. , This represents the posterior estimate of the square root factor of the error covariance matrix at the previous time step. Indicates the volume point. This represents the posterior estimate of the tracked state parameter at the previous time step. Represents the first point set 1 List.

6. The joint equalization method for PMD and RSOP damage according to claim 5, characterized in that, Step S3 includes the following steps: S301, Based on the square root factor Volume point and predicted value The volume point of the measurement update state of the square root capacitive Kalman filter (SCKF) is calculated using the following formula: ； S302, Based on the updated volume points The calculation yields the first polarization-damped joint equalization result. k A sequence of output sliding window signals; S303, according to the balanced... k The measurement prediction values ​​of the square root capacitive Kalman filter (SCKF) are calculated using the output sliding window signal sequence and the measurement equation of the fused sliding window. : in, Indicates the number of state parameters. The sequence number representing the volume point. , Indicates the predicted volume point value. This represents the modulo operation of complex numbers. and They represent the first k The first sliding window Group X and Y The intermediate element of the polarization volume point sequence, This indicates the rounding up operation.

7. The joint equalization method for PMD and RSOP damage according to claim 6, characterized in that, Step S302 includes the following steps: S3021. Capture a segment of signal using a sliding window. ,in, and They represent the first Damaged portion captured within a sliding window and Polarization signal sequence; S3022. Perform a Fourier transform on the intercepted signal to convert the signal to the frequency domain; S3023. Based on the signal converted to the frequency domain, the volume points after the update state are propagated using the PMD compensation matrix to obtain the measurement prediction volume points after PMD compensation. The expression for the PMD compensation matrix is ​​as follows: in, Represents the PMD compensation matrix. Indicates the use of the first The PMD compensation matrix with the volume points of the group state update as parameters The calculated differential group delay (DGD) Indicates according to the first i The vector of PMD obtained by using the group volume points as parameters. Indicates Stokes components In the measurement update state of the i The value corresponding to the volume point in the column. Indicates Stokes components In the measurement update state of the i The value corresponding to the volume point in the column. Stokes components In the measurement update state of the i The value corresponding to the volume point in the column; S3024. Transform the measurement prediction volume points after PMD compensation to the time domain through inverse fast Fourier transform. S3025. The measurement prediction volume points after PMD compensation converted to the time domain are propagated using the RSOP tracking matrix to obtain the measurement prediction volume points after PSOP tracking. The expression for the RSOP tracking matrix is ​​as follows: in, Indicates the first i The state parameters corresponding to the column volume points are used as the tracking matrix of RSOP. Indicates the phase rotation angle In the measurement update state of the i The value corresponding to the volume point in the column. Indicates the phase rotation angle In the measurement update state of the i The value corresponding to the volume point in the column. Indicates azimuth rotation angle In the measurement update state of the i The value corresponding to the volume point in the column; S3026. Calculate the average value of the predicted volume points after PSOP tracking to obtain the first value after joint equalization of polarization damage. k A sequence of output sliding window signals; The expression for the average calculation operation is as follows: in, and They represent the final joint equilibrium, respectively. X and Y The output sequence of the polarization sliding window Indicates the number of state parameters. Indicates the sequence number of the volume point. , and They represent passing through and RSOP tracking matrix After joint equilibrium X polarization and Y Polarization volume point sequence, Indicates the use of the first Volume points of column measurement update status The first three parameters are used as the matrix of parameters for the PMD compensation matrix.

8. The joint equalization method for PMD and RSOP damage according to claim 7, characterized in that, Step S4 includes the following steps: S401. Estimate the square root coefficients of the new information covariance matrix. : in, This represents the orthogonal triangular QR decomposition operation. Represents the weighted center matrix, Represents the noise covariance matrix The square root factor, Indicates the number of state parameters. Indicates the sequence number of the volume point. , Indicates the predicted measurement value. Indicates the predicted volume point value. Represents the noise covariance matrix The square root factor, ; S402, Estimating the cross-covariance matrix : in, and Both represent center-weighted matrices. The volume point represents the measurement update state. This represents the damage prediction value of the state parameter; S403, Based on the square root coefficient and cross-covariance matrix The Kalman gain was calculated. : S404. Combining the radius distribution of the received symbols, i.e., the Rice distribution model, to predict the measurement value. Perform a ring decision operation based on maximum a posteriori probability to obtain the optimized radius of the ring belonging to the constellation, and calculate the innovation matrix; S405. Based on the innovation matrix, combined with damage prediction values Predicted measurement values and Kalman gain The posterior estimates of the state parameters are calculated. : in, Indicates the first k The target observation values ​​of a sliding window, This describes the process of calculating the innovation matrix. S406. Update the root mean coefficient of the error covariance matrix according to the following formula. : in, Represents the center-weighted matrix. This represents the center-weighted matrix.

9. The joint equalization method for PMD and RSOP damage according to claim 8, characterized in that, The innovation matrix calculation process in step S404 includes the following steps: S4041. Combining the Rice distribution model, the measurement prediction value is calculated based on the maximum a posteriori probability. Optimized radius of the constellation ring with the highest probability of belonging Among them, the predicted measurement value The process for determining which ring belongs to which ring is as follows: in, This indicates that when searching for the maximum value of a function, m The value of , This represents each of the nine rings in the ideal PCS-64QAM constellation diagram. This indicates that the received signal amplitude is Under the condition that the amplitude of the transmitted signal is The probability, This indicates a probability calculation operation. This indicates the radius of the constellation ring of the PCS-64QAM signal at the transmitting end. The square root capacitive Kalman filter SCKF represents the first... A sliding window or Predicted values ​​of polarization measurements Indicates the noise variance. Represents the 0th order modified Bessel function. Represents different radii of the transmitter after probabilistic integer coding. The probability of occurrence; S4042. Calculate the optimized radius using the average value of the Rice distribution. Among them, the optimized radius The expression is as follows: in, Represents the first-order modified Bessel function. express Exponentiation, Represents logarithmic operations; S4043, optimize the radius of the annulus with the highest posterior probability. Compared with the predicted value By subtraction, we obtain the innovation matrix e: in, Indicates PCS-64QAM system X The optimized radius of the annulus to which the polarization measurement prediction value belongs. Indicates PCS-64QAM system Y The optimized radius of the annulus to which the polarization measurement prediction value belongs. express The predicted value of polarization measurement is... , express The predicted value of polarization measurement is... , Indicates the length of the sliding window. and They represent the first k Damaged portion captured by a sliding window X or Y Polarization signal sequence.

10. A joint equalization apparatus for implementing the joint equalization method for PMD and RSOP damage according to any one of claims 1-9, characterized in that, include: The first processing module is used to construct a joint equalization model of polarization damage based on the characteristics of polarization damage in optical fiber links under thunderstorm scenarios, and initialize the parameters of the square root capacitive Kalman filter (SCKF). The second processing module is used to update the square root capacitive Kalman filter (SCKF) by taking the posterior estimate of the state parameter at the previous time step as the prior estimate of the state parameter at the next time step according to the prediction equation, so as to predict the damage value of the damaged signal collected by the current sliding window. The third processing module is used to update the square root capacitive Kalman filter (SCKF) according to the measurement equation of the fused sliding window, and obtain the measurement prediction value. The fourth processing module, in conjunction with the radius distribution of the received symbols, performs a ring decision operation based on the maximum a posteriori probability on the measurement prediction value to obtain the optimized radius of the constellation ring, calculates the innovation matrix, and combines the damage prediction value and the measurement prediction value to calculate the a posteriori estimate of the state parameter. The fifth processing module determines whether the iteration is complete. If so, the process ends, completing the joint equilibrium of PMD and RSOP damages. Otherwise, the posterior estimates of the state parameters and the square root coefficients of the error covariance matrix are sent to the second processing module, and the iteration count is calculated. k Update and proceed with the next iteration.

Images