Chinese Optics Letters, Volume. 23, Issue 9, 091201(2025)

High-accuracy distributed forward state-of-polarization reconstruction in optical fibers enabled by polarization-sensitive optical frequency domain reflectometry

Yidai Zhu, Xinyu Fan*, Yangyang Wan, and Zuyuan He
Author Affiliations
  • State Key Laboratory of Photonics and Communications, Shanghai Jiao Tong University, Shanghai 200240, China
  • show less

    We propose a method for reconstructing the distributed forward state of polarization (SOP) in single-mode fibers (SMFs) to solve the problem of an unpredictable blind box of SOP evolution in many applications such as fiber-optic parametric amplification systems. Using polarization-sensitive optical frequency domain reflectometry (POFDR) and a quaternion-based model to describe polarization changes, our approach achieves high spatial resolution and precision. By an improved iterative fitting algorithm, the mean square error (MSE) of forward SOP reconstruction for approximately 100 consecutive measurement points was reduced to below 0.1%. This method enables visualization of SOP dynamics along the fiber, offering critical insights for polarization-dependent systems.

    Keywords

    1. Introduction

    For optical fibers, normally, the state of polarization (SOP) can only be measured at the output end via a polarization analyzer or related techniques. Due to unevenness and residual stress during the drawing process, irregular random birefringence effects occur inside optical fibers, and the SOP of the transmitted lightwave in optical fibers will undergo random changes along the fiber due to environmental factors and mechanical vibrations. However, the evolution process of SOP in optical fibers is a very important problem for many applications such as optical parametric amplification, as the gain is strongly polarization dependent because of an intrinsic angular-momentum conservation[1]. Numerous studies were proposed to control the SOP, such as nonlinear polarization pulling (NLPP)[2] and active polarization control[3]. The traditional method considered the SOP changes in optical fibers as an unpredictable black box, and used a nonlinear-optical polarizer to control the SOP, such as nonlinear lossless polarizers (NLPs) or Raman polarizers[4]. Therefore, making such a black box visualized, which means obtaining the evolution process of SOP along the fiber, can better assist the parametric amplification research and control.

    Note that in the fiber-optic sensing field, the Rayleigh backscattered (RBS) lightwave can preserve the SOP information of the scattering points, which shows a potential distribution method for measuring the SOP. By measuring the SOP of RBS signals, the first polarization-sensitive distributed fiber-optic sensors (DFOSs) based on an optical time domain reflectometry (OTDR) system were proposed[5]. The backward SOP can be collected and calculated accurately using the above-mentioned polarization-sensitive DFOSs. Afterwards, the local linear birefringence vector, which indeed reflects the influence of optical fibers on polarization, can be calculated[6,7]. The forward SOP can be reconstructed through a certain algorithm[8]. However, a hypothesis must be pre-set, that within a certain distance, normally three times the distance of spatial resolution, the linear birefringence and the twist rate are evenly distributed[9]. As for the OTDR system, this will bring significant errors and the hypothesis cannot be easily satisfied.

    In this paper, we propose a method to reconstruct the SOP of forward lightwaves using a polarization-sensitive optical frequency domain reflectometry (POFDR) system. Compared with the OTDR system, the optical frequency domain reflectometry (OFDR) system can provide a much higher spatial resolution, which helps to calculate the evolution of forward SOP more accurately. Using a simple polarization-diversity scheme to realize the OFDR system, and improving the calculation and iterative fitting algorithm, the mean squared error (MSE) of forward SOP for approximately 100 consecutive measurement points has been reduced to a fairly low rate of less than 0.1%.

    2. Theory

    In this section, we will introduce the basic theory of the matrix-based analysis method and how to apply it to the OFDR system to reconstruct the forward SOP. The matrix-based analysis method was first proposed by Ellison and Siddiqui, who used it for polarization-sensitive optical time domain reflectometry (POTDR)[9]. The Stokes vector and Mueller matrix were used to describe the transmission process of polarized light in optical fibers. Normally, the Stokes vector is described by four-dimensional parameters, S=[S0,S1,S2,S3], where S0 represents the intensity and S1,2,3 describes the SOP. After the normalization of the Stokes vector and Mueller matrix, the SOP parameters are the essential elements to be used for the calculations. Ignoring the variations of intensity, the polarization-related transmission can be written as[10][S1S2S3]=[m11m12m13m21m22m23m31m32m33][S1S2S3],where mij is the element of the rotation matrix M.

    In an OFDR system, we assume that the probe lightwave has an initial SOP of Sin and the received RBS lightwave has an SOP of Sout. The forward-return transfer function of distance z can be written as[9]Sout(z)=Rreturn(z)Rfiber[τ(z),βL(z),Δl]Rforward(z)Sin,where Rreturn and Rforward represent the Mueller matrices of the return path and forward path till distance z, and Rfiber describes the influence on SOP of the fiber with length Δl at distance z. Here, τ is the twist rate and βL is the linear birefringence. During the round-trip transmission, the fiber with length Δl rotates the SOP twice around the vector βF with a magnitude of |βF|Δl. When the fiber is reciprocal and geometrically asymmetric, the vector βF and transmission matrix Rfiber can be written as[11]βF=[βLx(z),βLy(z),gτ(z)]T,Rfiber[τ(z),βL(z),Δl]=RT[βF(z),Δl]RMR[βF(z),Δl],where g is the elasto-optic coefficient, R is the above-mentioned rotation matrix, and RM is a mirror matrix representing the reflection point. Combining Eq. (2) and Eq. (4), the Mueller matrix M in Eq. (1) can be factored into the product of multiple rotation matrices and solved numerically when Sin and Sout are known. According to Eq. (1), the transmission function is a 3×3 rotation matrix, providing three degrees of freedom, which requires three sets of linearly independent Sin and Sout to obtain the numerical solution. Moreover, there are nine unknown parameters in Eq. (2) (three degrees of Rreturn, three of Rforward, and three of βF), which means that three different Mueller matrices are required to perform the calculation. Assuming that τ and βL are stable within a range of 3Δl, three 3×3 matrix equations can be acquired by expanding the range. By combining three matrix equations and analyzing three sets of SOPs (Sin and Sout) within the range of 3Δl, we achieved a nine-order nonlinear equation system [Eq. (5)], which can be solved using the quasi-Newton method: {M1Sin=Rreturn(z)Rfiber[τ(z),βL(z),Δl]Rforward(z)SinM2Sin=Rreturn(z)Rfiber[τ(z),βL(z),2Δl]Rforward(z)SinM3Sin=Rreturn(z)Rfiber[τ(z),βL(z),3Δl]Rforward(z)Sin.

    Through the above calculation, the whole fiber has been divided into several transmission units with length Δl and its corresponding transmission matrix Ri(βFi,Δl). The forward SOP at distance z can be written as S(z)=n=N1Rn[βF(n*Δl),Δl]Rforward(0)Sin,where N=z/Δl is the number of transmission units and Rforward(0) is calculated by the first part of the fiber, which describes the influence of the system optical path on SOP before this part of the fiber. Note that the product operation starts from N and ends at one.

    Obviously, the quasi-Newton iteration of the nine-order nonlinear equation system is a massive and complex computational task. We make a simulation test of barely 500 points, which costs over 10 min. Performing iteration fitting on each element of Eq. (5) brings a large amount of redundant computing. Actually, the truly useful iterative operations should be performed only on the three required parameters in Eq. (3). The problem is that matrix M describing the rotation process is not concise enough. Nine elements only provide three valid degrees of freedom information: two of rotation vector orientation and one of rotation angle size. Therefore, we consider introducing another description method, called quaternion, to describe the transmission process. Using this description method, we can reduce the required parameters for a single fitting process by three times and greatly reduce the computation complexity. Besides, it is unable to perform averaging processes on the transmission matrix to reduce the errors but this reduction can be performed if we adopt quaternions. In Sec. 3, we proved the efficiency and accuracy of this method in reducing fitting errors. The accuracy of the forward SOP is improved based on the proposed description method and the improved algorithm.

    In a two-dimensional plane, a regular complex number of unit length z=eiθ can describe a rotation process around the origin. Similarly, the “extended complex number”, quaternion of unit length q=e(vxi+vyj+vzk)θ/2, can encode rotations in the three-dimensional space. Normally, quaternions have three imaginary unit numbers (i,j,k) that can be written as[12]q=qω+qxi+qyj+qzk=qω,qv,i2=j2=k2=1,ij=ji=k,ki=ik=j,jk=kj=i.When qω equals 0, p (=0,pv) can be called a pure quaternion and obviously can be used to describe the SOP. The conjugate of a quaternion q* and the norm of a quaternion q are defined by[12]q*=qω,qv,q=qq*=qω2+qx2+qy2+qz2,where represents the quaternion product, and the product of any two quaternions is shown as pq=[pωqωpvTqvpωqv+qωpv+pv×qv].

    The norm of unit quaternions equals one, and each encodes a rotation around the vector u with a magnitude of θ, which can be written as q=cosθ+usinθ,where u=uxi+uyj+uzk is a unit vector. When this rotation is applied to p, it can be described through a quaternion product: q=cos(θ/2)+usin(θ/2),p1=qpq*.

    We substitute the quaternion into Eq. (6), and the forward SOP can be written in the form of a quaternion product: S(z)=n=N1qnqforward(0)Sinqforward*(0)n=1Nqn*,where qn=[cos(θn/2)Bnsin(θn/2)],θn=|βF(nΔl)|Δl,Bn=βF(nΔl)|βFi(nΔl)|.

    The forward-return transfer function in Eq. (2) can be obtained through a similar quaternion product: Sout(z)=qreturn(zΔl)qNqNqforward(zΔl)Sinqforward*(zΔl)qN*qN*qreturn*(zΔl),where q=qω,Dqv, D=diag(1,1,1).

    The fitting element is directly correlated to the required parameters, instead of numerical changes in the rotation matrix. The parameters during the iteration process decreased from nine in the third-order matrix to three in the quaternions. Moreover, this is helpful for the subsequent discussion on the accuracy loss caused by fitting errors, which is shown in Sec. 3.

    Normally, a Stokes analyzer is a useful equipment to collect the SOP of RBS lightwaves, which is widely used in the POTDR system[13]. However, the OFDR system requires receiving devices with a much higher sampling rate, far beyond that of the conventional commercial Stokes analyzers. Actually, the phase information can be obtained during the time-frequency domain conversion of the OFDR system, which can help calculate the SOP using the following calculation method after adopting a polarization diversity receiver (PDR)[14]: [S0(z)S1(z)S2(z)S3(z)]=[|Ap(z)|2+|As(z)|2|Ap(z)|2|As(z)|22|Ap(z)||As(z)|cos(δ)2|Ap(z)||As(z)|sin(δ)],δ=arg[As(z)Ap*(z)],where Ap(z) and As(z) represent the amplitudes at position z of the two polarization directions received by PDR.

    3. Experiment and Discussion

    The experimental system’s structure is depicted in Fig. 1. It is a conventional OFDR system using a polarization diversity scheme. The SOP of the probe lightwave is controlled by an electrical polarization controller (EPC). The reference path (Lref) with a length of around 100 m is used for phase noise compensation[14], and the length of FUT is around 50 m. Figure 2(a) shows the frequency domain information collected from two polarization directions. The quasi-Newton fitting method is used to calculate the mathematical solution of Rfiber. The forward SOP along the fiber is then obtained after using Eq. (13), and part of the results is shown by the red line in Fig. 2(b).

    Structure of the POFDR system. SSB: single-side band; RF, radio frequency; EPC, electrical polarization controller; PDR, polarization diversity receiver; BPD, balanced photodetector; FUT, fiber under test; Lref, reference fiber.

    Figure 1.Structure of the POFDR system. SSB: single-side band; RF, radio frequency; EPC, electrical polarization controller; PDR, polarization diversity receiver; BPD, balanced photodetector; FUT, fiber under test; Lref, reference fiber.

    Experimental results. (a) POFDR curves for two polarization directions. (b) Comparison of forward SOP results calculated by POFDR and POTDR systems when the birefringence changes abruptly within a range of 4.11 m. POTDR: black points; POFDR: red line.

    Figure 2.Experimental results. (a) POFDR curves for two polarization directions. (b) Comparison of forward SOP results calculated by POFDR and POTDR systems when the birefringence changes abruptly within a range of 4.11 m. POTDR: black points; POFDR: red line.

    We also constructed a traditional POTDR system to measure the linear birefringence as a comparison. The result is shown in Fig. 3. It is clear that the numerical values obtained from the two systems are very close, which validated the accuracy of the algorithm. However, due to the difference in spatial resolution, POTDR can only describe a rough range and changing process, while ignoring the short-range changes of linear birefringence over a large scale. We select part of the forward SOP obtained by POTDR and POFDR where the birefringence has abrupt changes within a range of 4.11 m, and the result is shown in Fig. 2(b). The red line is the forward SOP calculated by POFDR while the black points are collected via POTDR. Although they maintained a similar change trend, the results of POTDR have significant errors, due to its neglect of sudden changes in birefringence. Therefore, it is necessary to improve the measurement spatial resolution for the reconstruction of forward SOP with higher precision. According to Fig. 2, the spatial resolution of forward SOP calculated by POFDR is around 3.74 cm, far exceeding the spatial resolution precision of POTDR in meters.

    Comparison of birefringence measurement results between POFDR and POTDR systems.

    Figure 3.Comparison of birefringence measurement results between POFDR and POTDR systems.

    During the experiment and calculation process, the backward SOP is obtained and the forward SOP is reconstructed. Then we proposed a method to assess the accuracy of the results. First, we define the deviation of SOP (DSOP) at a certain point as the modulus of the difference vector between the measured and calculated SOP. The total deviation, which is also called mean squared DSOP (MS-DSOP), is the mean square of all the forward DSOP along the fiber. Unfortunately, only the SOP of the emitted lightwave at the fiber end can be measured. Therefore, we use the backward DSOP to indirectly assess the accuracy. The calculation process is divided into the following steps: Input three sets of linearly independent correlation SOPs to acquire the parameters τ, βL, and transmission matrix;Adjust the EPC to set the input SOP to Sin1 and collect the output backward SOP sequence Sout1;Calculate the forward SOP sequence S and the backward SOP sequence Sout1 using Eqs. (13) and (15);Obtain the DSOP of Sout1 and Sout1 to indirectly assess the accuracy.

    Since the forward SOP is calculated by the parameters τ and βL, the assessment of the accuracy of the forward SOP is essentially the assessment of parameters τ and βL, which can also be revealed by the calculated backward SOP. Therefore, the accuracy of forward SOP can be represented by the backward DSOP.

    According to Eq. (13), the forward SOP at every point is accumulated from all the previous transmission matrix segments. The stop condition of the quasi-Newton iteration is set as the DSOP between the backward SOP Sout calculated via Eq. (15) and the collected true Sout to be less than 0.1%. The selection of this numerical accuracy is based on the precision that can be achieved using EPC to regulate the SOP in this experiment. Although the DSOP of each unit can be reduced to 0.1%, the total deviation can be very large. We simulate the process of deviation accumulation, and the result is shown in Fig. 4(a). Four thousand transmission units are used for simulation and the process is repeated 100000 times. Actually, the total deviation of N transmission units obeys the Wiener process with a mean value of 0 and a variance of N2σ, where σ=0.01%. This means that there is a 99.7% probability that the accumulated deviation of N units will fall in (3σN,3σN). Figure 4(b) shows the histogram of the DSOP at the fiber end with 100000 times repetitive simulations. Taking N equal to 4000 as an example, the range will be around ±18.97%, which is an unacceptable value. To solve the problem, another iteration termination condition is proposed to be added to the original quasi-Newton method. The forward transmission quaternion at a distance of z can be written as follows according to Eqs. (13) and (15): qforward(z)=n=N1qnqforward(0),where N=zΔl. The left part of Eq. (18) can be calculated by the backward SOP collected at a distance of z, and the right part is simulated by the previous iteration results. By adding a deviation of less than 1% between two terms as another termination condition to the iteration process, the accuracy of calculation can be improved greatly.

    Simulation results. (a) Process of deviation accumulation of forward SOP. (b) Histogram of DSOP at the fiber end.

    Figure 4.Simulation results. (a) Process of deviation accumulation of forward SOP. (b) Histogram of DSOP at the fiber end.

    Then we discuss the accuracy loss caused by fitting errors. Suppose the error between the normalized calculated birefringence vector B and the real value B is B=B+Δ.

    Within a distance of Δl, considering that the unit quaternion can be rewritten in exponential form through Euler’s equation, we can rewrite Eq. (13) as S(z)=qNS(zΔl)qN*=eθ2BS(zΔl)eθ2(B).

    We use vector u to replace the previous forward SOP S(zΔl) and expand Eq. (20): S(z)=eθ2Bueθ2(B)=eθ2(B+Δ)ueθ2(BΔ)=(eθ2B+sinθ2Δ)u(eθ2(B)sinθ2Δ)=eθ2Bueθ2(B)+sinθ2(Δueθ2(B)eθ2(B)uΔ)sin2θ2ΔuΔ.

    The first term of Eq. (20) is the true value of forward SOP we need, and the remaining terms are the errors. The Δ and u are pure quaternions, and their product is Δu=[ΔT·uΔ×u],uΔ=[uT·Δu×Δ].

    So Δu=(uΔ)*, and the norm of the fourth term equals sin2θ2|Δ|2|u|sin2θ2|Δ|2, which is the higher-order term of error, and can be ignored. The second and third terms of Eq. (21), ignoring the term sinθ2, can be rewritten as follows through substitution of p=Δu,q=eθ2B: Term2+Term3=pq*qp*=pq*(pq*)*=[02(pq*)v]=2[pω(qv)+qωpv+pv×(qv)]=2[(ΔT·u)sinθ2B+cosθ2Δ×usinθ2Δ×u×B].

    The error Δ is a tiny vector whose direction follows a uniform distribution, so that Δi=0. Therefore, we can eliminate the second and third terms of Eq. (23) by averaging multiple fitting results: {1NiN(cosθ2Δi×u)=cosθ21N(Δi)×u=01NiN(sinθ2Δi×u×B)=sinθ21N(Δi)×u×B=0.

    By an averaging operation, the first term of Eq. (23) turns into 1NiN(2ΔT·usinθ2B)=2sinθ2B1NiN(ΔT·u)=2sinθ2B·1NiN(Δixux+Δiyuy+Δizuz)=2sinθ2B·1N(Δixux+Δiyuy+Δizuz)=0.

    Therefore, by averaging multiple fitting results, further improvements can be efficiently achieved on the forward SOP.

    Figure 5 shows the contrast of part of the backward DSOP, which assesses the accuracy of the forward SOP using the original iteration process and the improved algorithm with an averaging operation. The red line represents the result of the original algorithm and the blue line represents that of the improved algorithm. In the original fitting algorithm (the red line), due to limitations in computational complexity, the stop condition is set as the DSOP between the calculated backward SOP Sout and the collected true backward SOP Sout to be less than 1%. Without changing the accuracy of the fitting termination condition, we incorporated the calculation accuracy of transmission quaternions [Eq. (18)] into the iteration conditions. Meanwhile, according to Eqs. (21)–(25), we averaged 100 sets of iteration results under different initial iteration inputs. The proposed algorithm has improved SOP accuracy by at least one order of magnitude, reaching a level of 0.1%, which can be observed from the blue line. We calculated the MS-DSOP represented by different algorithms, and it decreases sharply from 3.85% (red lines) to 0.0475% (blue lines). Due to the accuracy limitations of the EPC and Stokes analyzers used in our laboratory, the SOP accuracy we obtained from the laboratory can only reach the level of 0.1%. This limited the final accuracy. However, compared to the 7% error level calculated by previous researchers[15], significant improvement was achieved using the proposed algorithm.

    Contrast of absolute DSOP value for approximately 100 points. DSOP1: DSOP between the collected backward SOP and the calculated SOP using the improved algorithm. DSOP2: DSOP of the backward SOP and the calculated SOP using the original algorithm.

    Figure 5.Contrast of absolute DSOP value for approximately 100 points. DSOP1: DSOP between the collected backward SOP and the calculated SOP using the improved algorithm. DSOP2: DSOP of the backward SOP and the calculated SOP using the original algorithm.

    4. Conclusion

    We proposed a method to calculate the SOP of forward lightwaves in optical fibers with high accuracy. Using the backscattered SOP measured by the POFDR system, the SOP of forward lightwaves in the fiber can be reconstructed and visualized, which has important research significance for the fiber parametric amplification technology. With the improved iterative fitting algorithm, the MS-DSOP of a forward SOP with over 100 points decreased to a low rate of 0.0475%, which is a tenfold improvement over OTDR-based methods, demonstrating the superiority of POFDR in high-accuracy SOP reconstruction.

    [12] J. Sola. Quaternion kinematics for the error-state Kalman filter(2017).

    [15] A. Galtarossa, L. Palmieri. Mapping of intense magnetic fields based on polarization sensitive reflectometry in single mode optical fibers. AFRICON Conference, 1(2013).

    Tools

    Get Citation

    Copy Citation Text

    Yidai Zhu, Xinyu Fan, Yangyang Wan, Zuyuan He, "High-accuracy distributed forward state-of-polarization reconstruction in optical fibers enabled by polarization-sensitive optical frequency domain reflectometry," Chin. Opt. Lett. 23, 091201 (2025)

    Download Citation

    EndNote(RIS)BibTexPlain Text
    Save article for my favorites
    Paper Information

    Category: Instrumentation, Measurement, and Optical Sensing

    Received: Feb. 28, 2025

    Accepted: May. 9, 2025

    Published Online: Aug. 21, 2025

    The Author Email: Xinyu Fan (fan.xinyu@sjtu.edu.cn)

    DOI:10.3788/COL202523.091201

    CSTR:32184.14.COL202523.091201

    Topics