Ultrasound-guided prostate percutaneous intervention robot system and calibration by informative particle swarm optimization

Applying a robot system in ultrasound-guided percutaneous intervention is an effective approach for prostate cancer diagnosis and treatment. The limited space for robot manipulation restricts structure volume and motion. In this paper, an 8-degree-of-freedom robot system is proposed for ultrasound probe manipulation, needle positioning, and needle insertion. A novel parallel structure is employed in the robot system for space saving, structural rigidity, and collision avoidance. The particle swarm optimization method based on informative value is proposed for kinematic parameter identification to calibrate the parallel structure accurately. The method identifies parameters in the modified kinematic model stepwise according to parameter discernibility. Verification experiments prove that the robot system can realize motions needed in targeting. By applying the calibration method, a reasonable, reliable forward kinematic model is built, and the average errors can be limited to 0.963 and 1.846 mm for insertion point and target point, respectively.


Introduction
Prostate cancer is a widespread disease that threatens men's health [1]. Transrectal ultrasonography (TRUS)guided biopsy and brachytherapy are effective approaches for prostate diagnosis and treatment [2]. To accomplish these procedures, feasible and reliable prostate percutaneous intervention is needed. Generally, the physician operates an ultrasound probe and inserts a needle simultaneously with the help of a tablet. The insertion trail is restricted in the ultrasound probe imaging plane to help the physician reduce insertion difficulty, but the restriction reduces insertion flexibility [3]. Moreover, the insertion procedure needs repetitive attempts and greatly depends on the operator's experience. To reduce experience dependence, manipulation difficulty, and patient trauma, and increase insertion accuracy and manipulation flexibility, image-guided percutaneous intervention robots are developed [4], and an automatic ultrasoundguided prostate intervention is the basis for further medical image fusion application [5,6].
Researchers have proposed different varieties of prostate robot systems. Robot "Apollo" was designed with three motors and three brakes to manipulate the endorectal ultrasound probe [7]. The structure employed transrectal access and could realize remote center motion. A 9-degree-of-freedom (9-DOF) structure for prostate brachytherapy was proposed to manipulate the ultrasound probe and adjust the tablet for needle attitude guidance [3]. The structure worked with a high flexibility, but the overlarge size limited its practical application. In Ref. [8], a 4-DOF hands-free probe manipulator for TRUS-guided prostate biopsy was designed. The needle driver employed transrectal access and was coupled with the probe motion. In addition to the customized manipulators, commercial robots were employed as assistance to manipulate the probe in prostate therapy [9]. However, commercial robots may not meet the manipulation or sterilization requirement.
Though several researchers tried to develop a control method for a system with unknown parameters [10], identifying robot system parameters by calibration is usually fundamental for effective, intuitive robot control [11]. Calibration contains two steps: kinematic modeling, and parameter identification and compensation [12]. Many mature methods are used for serial robot calibration, but the situations for parallel structure calibration are much more complicated [13]. Parallel structure models vary and are difficult to be summarized as a universally applicable model. The kinematic parameters of parallel structures are coupled nonlinearly, and the end pose is affected by the aggregation of manufacturing and assembly errors. Methods based on the Denavit-Hartenberg (DH) models are the most extensively used approaches for kinematic modeling in real application [14,15]. The modified DH and Gauss-Newton method were applied to calibrate the parallel manufacturing machine [16]. To distinguish errors from different sources, a generalized Jacobian method was proposed [17]. The error model of the parallel mechanism was built by screw theory in Ref. [18], and the dual-vector space method was applied to distinguish errors from different sources. However, this method was based on a one-order linearization. Reference [19] designed a camera calibration technique by collecting data of a flat pattern in several poses and applied the Levenberg-Marquardt algorithm as the optimization method to identify the parameters; however, this method may fall into local optimum. Reference [20] identified the structural parameters of a 3-DOF overconstrained parallel robot by the nonlinear least-squares method and the nongeometric parameters by the trained neural network method. Reference [21] developed a local convergence method and applied Tabu Search to calibrate Gough platform. Several novel approaches are flexible and effective for solving traditional problems in different regions [22,23]. Particle swarm optimization (PSO) is a nonlinear optimization method with rapid calculation speed [24] and can be applied to solve optimization problems in parameter identification [25,26]. A hybrid algorithm based on neural network and PSO was applied for industrial robot kinematic parameter identification [27]. Reference [28] deployed a ball-plate-based calibration approach and applied PSO for parameter identification. The PSO method was intuitive and effective, but the calculation cannot achieve real-time supervision and was stochastic and easily fell into local minimum [29]. Calibration methods based on parameter characteristics were proposed to increase identification accuracy. The global error transformation index was related to the global maximum pose error, and genetic algorithm based on the index was applied for 3-DOF parallel robot parameter optimization [30]. Sensitivity-based parameter calibration was proposed for limited observations, and the results showed better parameter estimation and prediction accuracy compared with the least squares calibration method and the Bayesian calibration method [31]. Reference [32] performed optimization and calibrationbased model validation in a sequence of small domains in the entire design space and applied Gaussian distribution and maximum likelihood estimation to calibrate the prediction model. Reference [33] proposed a calibration method focused on the least error sensitive regions of parallel kinematic machines.
In this paper, the main contribution is the design of an innovative 8-DOF robot system for transrectal ultrasoundimaging-guided prostate percutaneous intervention and a novel optimization method for kinematic model calibration, named informative particle swarm optimization (InfoPSO). The robot system is compact for transrectal ultrasound probe manipulation, needle positioning, and needle insertion. A parallel structure is employed in the robot system to increase structural rigidity and avoid the potential risk of confliction between the robot body and the patient. To model the 4-DOF parallel structure accurately, the structure kinematic model is modified to contain manufacturing and assembly errors. The error parameters are identified by the proposed InfoPSO method. In the method, the error parameters are grouped based on the order of magnitudes of parameter discernibility, and PSO is carried out stepwise. Then, targeting experiments are carried out to verify the structure function and mechanical accuracy.
The paper is organized as follows. The structure design of the 8-DOF percutaneous intervention robot system is introduced in Section 2. The kinematic calibration method InfoPSO is proposed to identify the parallel structure parameters in Section 3. Experiments are carried out for structure calibration and targeting accuracy verification in Section 4. The paper is concluded in the final section. The calculation base of parameter discernibility is introduced in Appendix A.

Percutaneous intervention robot system design
Robot-assisted percutaneous intervention is an effective approach to realize needle insertion. According to humanity anatomic structure and surgery requirement, the workspace should at least contain the following region: The needle insertion point can reach a region of circle with 50 mm diameter; the needle insertion depth is larger than 70 mm from the insertion point; and the needle insertion angle should reach a circular cone of apex angle 50°. According to intervention requirements, the targeting accuracy of a prototype should be higher than 2 mm for a successful surgical procedure [34].
In this paper, an 8-DOF robot system is designed for ultrasound-imaging-guided prostate percutaneous intervention, and the scenario is shown in Fig. 1. The patient is in a lithotomy position, and the robot system uses the space between the two legs. The robot system is designed to realize hands-free activities of probe manipulation, needle positioning, and needle insertion. The decoupled motions allow the physician carry out more complicated and flexible activities. To simplify the robot system structure, the horizontal and vertical adjustments of the robot system base are accomplished by a manual lift table. As shown in Fig. 1, the robot system base is fixed on a lift table. By adjusting vertical and horizontal positions of the lift table, the ultrasound probe axis is aligned with the patient's anus. Then, the lift table is locked on the bed, and the lift table and the robot system base are kept fixed during the whole scanning and intervention process.
The robot structure is shown in Fig. 2. The robot dimension is 254 mm × 420 mm × 488 mm (the length 420 mm is measured excluding the two devices because the total length changes as the needle and the probe move). The robot contains a 2-DOF structure for ultrasound probe movement, as shown in Fig. 2(a), a 4-DOF parallel structure for needle positioning, as shown in Fig. 2(b), and a 2-DOF structure for needle insertion, as shown in Fig. 2(c). The bases of the ultrasound probe manipulator and the needle positioning manipulator are both fixed on the lift table, as shown in Fig. 1.

2-DOF ultrasound probe manipulator
The 2-DOF ultrasound probe manipulator is used to rotate and move the transrectal ultrasound probe along its axis, as shown in Fig. 3. The two DOFs can realize the motions needed for ultrasound probe scanning. Motors 1 and 2 are used to realize the rotation and the linear motion of the probe, respectively. A rubber block is wrapped around the probe for adjustment and locking. The rubber block is a quick-changing part that is easy to be sterilized.

4-DOF needle positioning manipulator
The 4-DOF needle positioning structure is used to orient and position the 2-DOF needle driver, as shown in Fig. 4. This part utilizes a parallel structure to increase rigidity and avoid collision between the moving parts and the patient. The structure mainly contains two similar stages assembled on the base, and each stage contains a scissor mechanism as end effector, as shown in Fig. 4(a).

Single-stage structure and driving system
The single-stage structure is shown in Fig. 4(b) and mainly contains two separate transmission systems (labeled by purple and blue colors separately in the figure). Each transmission is composed of motor module, timing belt, and disc module. The motor module contains motor, relative encoder, reducer, coupler, small pulley, related shafts, and mounting parts (motor: DCX16L, planetary gearhead: GPX16, sensor: ENX16, Maxon, Switzerland). The disc module contains a large pulley, which is driven by timing belt transmission.
Bearings are used to position the rotating disc modules. As shown in Fig. 4(c), positioning slide bearing is used to change timing belt shape and avoid confliction between timing belt 2 and the base part. Small guide slots are reserved on the supporting board to allow position adjustments of the motor modules, such that pretightening force can be exerted on the timing belts. Two stages have the same transmission structures, where four motor modules are used to drive four rotating disc modules.

Scissor mechanism end effector
The two stages contain similar end effectors, and each end effector is a scissor mechanism formed by four hinged arms. Two arms are hinged to the two rotating  disc modules of the same stage, as shown in Fig. 5(a). The relative motion between two hinged points impels the scissor mechanism end to extend or retract, as shown in Fig. 5(b). When two hinged points rotate together, the end trail can cover all the space inside the rotating disc. The needle driver passes through the two scissor mechanism ends of the top and bottom stages. As shown in Fig. 6, the end structure of the bottom stage is a customized Hooke joint. This structure restricts axial rotation and sliding motion, but allows rotation along the two other axes. The end structure of the top stage is a centripetal joint bearing, where the needle driver passes through the central hole, and can rotate freely in all directions and slide along the needle driver axis. The structure design allows distance changing between the two end effectors when the needle driver moves.

Position measurement devices
Limited by the mechanical structure, the absolute position of parallel structure is acquired by the combination of the relative encoder and the limit switch. The parallel structure contains four timing belt-pulley transmission. For each transmission, a relative encoder is assembled in motor module, and an optoelectronic switch (PM-L25, Panasonic, Japan) is assembled near the disc module, as shown in Fig. 7. A flag is mounted on the disc module, and the optoelectronic switch is mounted on the supporting board. To calibrate the system parameter, an optical tracking system is used, and the positions for spherically mounted retroreflectors (SMRs) are reserved on the rotating discs.

2-DOF needle driver
The 2-DOF needle driver realizes the rotation and the linear motion of the needle along its axis. The structure uses a compact design to decrease structure dimension    and avoid overturning. As shown in Fig. 8, the needle driver passes through the top and bottom scissor mechanism ends. As the two ends are positioned, the pose of the needle driver is oriented. Fixture 1 is hinged to the Hooke joint of the bottom stage; Fixture 2 passes through the centripetal joint bearing of the top stage; and the two fixtures are fixed by screws.
The needle driver contains two identical motor modules (motor: DCX10S, planetary gearhead: GPX10, sensor: ENX10, Maxon, Switzerland). For the linear joint, the mounting plate moves relative to Fixture 2 by gear rack pair transmission. For the rotary joint, the insert needle is fixed with a bevel gear; driven by motor module 2 and bevel gear pair transmission, the needle can realize selfrotation. To separate the needle and other parts of the needle driver, a sterilizable long tube is applied as a needle guide. The needle guide is fixed at the end of Fixture 1 by screws, in case of undesired hurts caused by needle guide motions during positioning.
Same as the 4-DOF needle positioning manipulator, the 2-DOF needle driver uses the combination of relative encoder and limit switch-flag pair to acquire joint absolute position. As shown in Fig. 9, the two limit switches are assembled on Fixture 2 and the mounting plate, and flags are assembled on the moving parts in each joint.
To sum up, the 8-DOF robot system is designed for ultrasound probe manipulation, needle positioning, and insertion. In the robot system, the 2-DOF probe manipulator is a reliable structure for probe axial rotation and insertion. The other 6-DOF structure is placed above the probe manipulator to decrease the system volume. The 4-DOF parallel structure can realize needle positioning without motor or cable movement. This design increases system rigidity and decreases the potential risk of undesired collision between the robot body and the patient. The 2-DOF needle insertion structure uses compact distributions to reduce overturning torque.

Kinematic calibration and parameter identification
In the robot system, the 2-DOF probe manipulator and the 2-DOF needle driver are serial structures with coaxial linear and rotation joints, and calibration can be accomplished by the combination of the DH method and the least squares algorithm. The 4-DOF needle positioning manipulator is a parallel structure and needs an effective method to calibrate error parameters. In this paper, the structure model is built first, and a calibration method for the 4-DOF needle positioning manipulator is carried out. The inputs of kinematics are the rotating angles of scissor mechanism arms. The geometrical relationships show that the position of end point T satisfies the following equation: where P T is the position vector of point T, f(·) is the position function of vector P T , θ 1 and θ 2 are the rotating angles of points A and B, respectively, d 1 and d 2 are the lengths of the scissor mechanism arms, l AC = l BC = d 1 , l CD = l DT = l TE = l EC = d 2 , and r is the radius of the circle where the hinged points slide. In Fig. 10(b), the geometric model of the two-stage scissor mechanism is presented, and both stages share the same models. The coordinate system is placed above the top stage, and the xOy plane is defined by the SMR assembly plane. The subscripts t, b and f indicate that the parameters belong to top stage, bottom stage and needle tip, respectively. Assuming a needle passes through top end T 1 and bottom end T 2 , the needle tip is denoted as T f , and then the position vector of needle tip T f can be calculated as follows: l, (2) P Ti where is the position vector of point T i , i ∈{1, 2, f}, and l is the distance between needle end T f and hinged point T 2 .

Robot kinematic model modification
To calibrate the system kinematic parameters, the kinematic model is reconstructed with manufacturing and assembly errors.
Single stages are considered first. The values of θ 1 and θ 2 are obtained by the sum of relative encoder readings and home positions. The home position can be obtained from the computer aided design model, and a more precise value should be identified by calibration. To present the manufacturing and assembly errors on the system precision, these error items are integrated into the single-stage kinematic model Eq. (1), as shown in the following equation: where is the modified position vector of point T, (·) is the modified position function of vector , d 1 , d 2 , r, θ 1 , θ 2 , and f(·) have same meanings as in Eq. (1), δj is the error of corresponding variable, where j∈{d 1 , d 2 , r, θ 1 , θ 2 }, and δx and δy are disc center errors with respect to the system base coordinate. For the single stage, seven error parameters, namely, {δθ 1 , δθ 2 , δd 1 , δd 2 , δr, δx, δy} need to be calibrated.
When the two stages are assembled as in Fig. 10(b), extra errors are introduced to the system. The errors mainly include position and coaxial errors, and the needle end tip equation is modified as Eq. (4): where is the modified position vector of point T i , i∈{1, 2, f}, and are the modified kinematic position functions of the top and bottom stage end vectors, respectively, {θ t1 , θ t2 } and {θ b1 , θ b2 } are the angle inputs for the top and bottom stages, respectively, and {δx t , δy t , δz t } and {δx b , δy b , δz b } are the scissor mechanism end position errors of the top and bottom stages, respectively. The distances between the xOy plane and top stage, bottom stage, and needle tip are denoted as z t , z b , and z f , respectively. ||·|| is the norm of the corresponding vector. The parameters for two-stage system calibration are {δx t , δy t , δz t , δx b , δy b , δz b }.

Calibration overview
Calibration is applied to identify the unknown parameters in the modified kinematic models. Considering the nonlinearity of the parallel robot kinematic model, PSO is applied in this paper during parameter identification for its rapid calculation speed. However, during iteration, the PSO method is not under a real-time supervision, and the final results may be converted to unreasonable values. To solve this problem, the PSO method is expected to be combined with the parameter characteristics.
In this paper, the informative value is combined with the PSO, and a method named InfoPSO is proposed. The informative value is a nondimensional value that reflects the discernibility of each parameter, and the calculation detail is explained in Section 3.2.3.
The scheme of InfoPSO is shown in Fig. 11. Each single-stage kinematic model has seven parameters that need to be identified. Combining the informative values with PSO, the modified single-stage model can be acquired. Based on the modified top-stage model, bottom-stage model, informative values of two-stage model parameters, and PSO, the modified two-stage model can be acquired.

PSO approach review
PSO is a group intelligence optimization algorithm imitating the process of birds searching for food based on information exchange. The main process is as follows: Generate primary particles. In the D-dimension solution space, N particles are defined randomly, and each particle is a candidate solution for the problem. The particle swarm is . For each particle, , where "0" represents the first generation, i is the particle number, and means the value of ith particle at jth dimension. The primary velocity of the ith particle is defined as .  For the searching, the fitness values are used to guide the particle swarm finding the target. The swarm fitness value gbest k is the particle that fits the objective function best in the particle swarm until the current generation k; for each particle i, the value fits the objective function best in its history is the individual fitness value . c) Update the particle swarm velocity. Referring to the fitness values of individual and swarm, the particle values are updated by corresponding velocity. For the ith particle, the velocity is where is the inertia weight reflecting the influence of current velocity, c 1 and c 2 are particles' cognitive and social factors, respectively, and µ and η are random values in (0, 1). The velocity of dimension j is restricted by a maximum velocity . If velocity , then ; if , then . d) Update particle value. Based on the current value and updated velocity, the next generation of particles can be obtained as e) Update the fitness values. The current fitness of each particle in the new generation is calculated. For each particle, the individual fitness value is replaced with the new fitness value if the latter is better. Among all the fitness values in the current generation, the best fitness value is searched and compared with the swarm fitness, and the better one is selected as the swarm fitness value.
Steps c)-e) are repeated until the swarm fitness value satisfies the task requirement or the cycle counter reaches the maximum value.

Informative value calculation
In PSO, particles are randomly generated, and each particle is a candidate for kinematic modeling. During searching, the fitness value is the only evaluation standard, and all the parameters are identified simultaneously. Though a large difference exists between two particle candidates, the fitness values may locate close to each other. The phenomenon results from that different parameters have varied influence rates or different "parameter discernibility." The changes led by parameters with less discernibility may be omitted by those changes led by parameters with more discernibility. To realize the identification accurately, parameters with different discernibility magnitude orders should be identified separately. The derivation of parameter discernibility is shown in Appendix A. In general, the derivation includes two steps: First, the gradients of all parameters belonging to {δθ 1 , δθ 2 , δd 1 , δd 2 , δr, δx, δy} are calculated; then, the high-order terms in the gradient are neglected for a convenient calculation. When high-order terms are ignored, the constant term coefficient is remarkable enough to show the discernibility of one parameter because the identified error parameters are all near 0. The constant term is defined as the "informative value" to reflect the parameter discernibility.
The informative values of all the parameters are shown in Table A1. For each parameter, the x and y components are functions of θ 1 and θ 2 , that is, the informative values of the error parameters may vary as the function inputs change. In theory, the definition domain of θ 1 and θ 2 are (−180°, 180°). Considering the restriction of parallel structure, θ 1 changes freely in (−180°, 180°), whereas θ 2 is limited in the range of (θ 1 +30°, θ 1 +50°). These input restrictions hold for all parameter informative value calculations.
The informative values of all the parameters are illustrated in Fig. 12. The informative values are labelled as I(·) for short, and for each parameter, the informative value contains components on the x and y axis and both are shown in Fig. 12. For a clear comparation, the maximum informative values of the parameters are listed in Table 1. Figure 12 and Table 1 show that the informative values vary greatly. The maximum values of I(δθ 1 ) and I(δθ 2 ) are much larger than the other parameters, which means the value fluctuations of δθ 1 and δθ 2 greatly influence the values of other parameters during identification. The informative values of δd 1 and δr are also larger than δd 2 , δx, and δy. If all the parameters are explored simultaneously, the parameters with larger informative values harmfully affect other parameter identification processes. The optimization should consider informative values and carry out the process stepwise.
As shown in Eq. (4), the two-stage system contains offset parameters {δx t , δy t , δz t , δx b , δy b , δz b }, where {δx t , δy t , δz t } and {δx b , δy b , δz b } are offsets for the top end and the bottom end, respectively. The error parameters are homogenous in Eq. (4), and their informative values are all constantly 1.0 (similar to δx and δy in Fig. 12). These six parameters can be identified in a single-step process.

PSO based on informative value
The 4-DOF needle positioning manipulator calibration is carried out, as shown in Fig. 13, and the two single stages are calibrated separately first. In this paper, InfoPSO is applied as the identification method. The informative values prove that different parameters have distinct discernibility, and the seven parameters are classified into three groups according to the magnitude orders of informative values, namely, {δθ 1 , δθ 2 }, {δd 1 , δr}, and {δd 2 , δx, δy}. Three groups of parameters are identified in a three-step identification sequence. When one group of parameters is identified, the other parameters are set to be 0 (not identified yet) or result values (already identified). The constants r, d 1 , and d 2 use designed values. θ 1 and θ 2 are obtained by the sum of theoretical home positions and real-time readings of relative encoders. The real end effector position P T is collected for each pair of {θ 1 , θ 2 }, and M groups of {θ 1 , θ 2 , P T } are collected and applied for PSO. The fitness value ε for each particle is defined as the quadratic criterion error: f * where (·) is the modified model as Eq. (3).
When one group of parameters changes, other parameters may be disturbed because of parameter interaction. Thus, a repetition of the three-step identification is needed until all of the parameters remain unchanged.
After the single stages are identified, the two-stage system is assembled, and {δx t , δy t , δz t , δx b , δy b , δz b } are parameters for identification. Equation (4) shows that the six parameters share the same informative value (all  Fig. 10(b). N groups of {θ t1 , θ t2 , θ b1 , θ b2 , , } are collected for PSO, and then the fitness value for one particle is defined as where and are defined as Eq. (4), is approximately regarded as the insertion point, α 1 and α 2 are the weight factors of target accuracy and insertion point accuracy, respectively, and α 1 + α 2 = 1.

Experiment
The robot system is assembled, and the platform is built for calibration, as shown in Fig. 14. The motors are driven by actuators (G-SOLWHI2.5/100EE, Elmo), and the controller is an industrial personal computer from Beckhoff. To collect position data, an optical tracking system is used (AT960, Leica, Switzerland).

Single-stage structure calibration
First, the two stages are analyzed separately. A single stage is assembled, and an equivalent weight is adhered to the top stage end effector to simulate the force exerted by the needle driver assembly. After the system powers on, the two rotating discs rotate together until both hit the limit switches. Then, the home positions are determined according to the limit switch positions and the predefined offsets of each joint. The end effector control is based on one-stage kinematics, and optical tracking system is used to collect position data of mechanism end, as shown in Fig. 15(a). Twelve traces are selected uniformly, and the positions of eight points are collected along each trace, as shown in Fig. 15(b). The process is repeated three times, and 288 groups of data are collected. The collected data are applied to the InfoPSO in Fig. 13, and the single-stage parameter errors can be identified.
The positions of targeting experiments at the top and bottom stages are shown in Fig. 16. The single-stage parameter identification results calculated by the traditional PSO and InfoPSO are shown in Tables 2 and  3, and the designed values are listed. To compare the results intuitively, the square root of the quadratic criterion error in Eq. (7) is defined as "average error" to judge the results. The average error is focused on the mechanical error, and no ultrasound image error is considered. Thus, the data collected by the optical tracking system can be used to calibrate the kinematic parameters.  In Tables 2 and 3, initially, the parameters are all set to be 0, and the average errors are 2.778 and 1.543 mm for the top stage and the bottom stage, respectively. By applying InfoPSO, the average errors are reduced to 1.499 and 0.975 mm for the top stage and the bottom stage, respectively. The traditional PSO and InfoPSO contain similar results for δθ 1 and δθ 2 (parameters with high informative values), whereas the other parameters are remarkably different. Though the average error of traditional PSO is smaller than that of InfoPSO, the length errors δr, δd 1 , and δd 2 identified by the traditional PSO are larger than 10 mm. For example, the δr of the top stage is approximately 59 mm, which means that the manufactured disc radius is 59 mm larger than the designed model. The unreasonable results suggest that the identified results of the traditional PSO are only fit for the collected data group. If another data group is collected in additional experiments, the large error parameters may lead to an unreliable end position. As an alternative choice, the results of InfoPSO are more trustworthy and applicable for robot kinematics. The reasonability is also proven by the fact that the same manufacturing errors of the two stages are close to one another. The average error of the top stage is larger than that of the bottom one because the added weight is exerted on the top stage (simulating the needle driver gravity), and the force leads to a larger deviation to the end effector position.

Two-stage structure calibration
The two stages and the needle driver are assembled together. As mentioned above, six parameters {δx t , δy t , δz t , δx b , δy b , δz b } are identified simultaneously. The coordinate system is built above the top stage, as shown in Fig. 10, and 21 points are selected uniformly from both stages, as shown in Fig. 17. For each point from the top stage, the neighbor points in the bottom stage are selected as point pairs. A total of 141 groups of point pairs are selected. For each pair, the needle driver is oriented by scissor mechanism end effectors, and {θ t1 , θ t2 , θ b1 , θ b2 , , } are collected. As in Eq. (8), the fitness value is defined as the average of the quadratic criterion error at the bottom stage point (near body insertion point) and the needle tip point.
The identified parameters and their designed values are listed in Table 4. The average error is a judgement of accuracy at the insertion point and the needle tip, and is defined as the square root of ε in Eq. (8). The weight To sum up, the parameters for kinematics modification are listed in Table 5. δθ 1 and δθ 2 are disc module offsets based on theoretical home positions, is the modified hinged point circle radius, and are the modified scissor mechanism arm lengths, δx and δy are the disc module center errors with respect to the system coordinate system, and is the z component of the end effector position with respect to the system coordinate system. All the modified parameters are substituted into system kinematics, and a verification experiment is carried out as follows.

System verification experiment
The identified parameters in Table 5 Fig. 18. The black and red points show the distribution of results calculated by original kinematics and modified kinematics, respectively, in contrast to the collected points with blue color. In Fig. 18(b), grey circles are applied to indicate the data groups.
The figures show that the points obtained by modified kinematics (red) are located closer to the measured data (blue). In Fig. 18(a), the average error of original kinematics is approximately 1.497 mm, and modified kinematics shows an error of 0.963 mm. In Fig. 18(b), the average error of the original kinematics is approximately 2.654 mm, and modified kinematics shows an error of 1.846 mm. The error of the needle tip is reduced by 43.8%, which proves that the kinematic modification substantially improves accuracy. Compared with similar robots for prostate percutaneous intervention [34], reaching a targeting accuracy smaller than 2 mm is acceptable for the parallel structure prototype.
For a further comparison, the boxplots of the system targeting experiments are shown in Fig. 19. The  Fig. 19(a) shows that the results of the bottom stage are better than those of the top stages because the needle driver exerts more force on the top stage, and gravity error (which is not considered in this paper) has more influence on the top stage, especially near the hinged circle center.

Conclusions
In this paper, an 8-DOF ultrasound-guided prostate percutaneous intervention robot system is designed to Fig. 17 Coordinate system distribution and needle driver orientation selection for two-stage targeting accuracy experiment.  assist physicians in realizing hands-free ultrasound probe scanning, needle pose adjustment, and needle insertion. The structure is novel with high rigidity and compactness, and can help the probe and the needle realize the required movements as well as decrease the potential risk of confliction between the robot body and the patient. To calibrate the parallel structure and identify the error parameters with different discernibility, InfoPSO is proposed. In stepwise InfoPSO, the single-and two-stage parameters are identified to modify the kinematics. The targeting experiments show that the identified parameters for the kinematic model modification are reasonable. By applying the modified kinematic model, the parallel structure can reach average errors of 0.963 and 1.846 mm at needle insert position and target position, respectively. The error of the needle tip is reduced by 43.8% compared with that of the initial kinematic model. The improvement proves that the calibration method is practicable, and InfoPSO provides an available, trustworthy solution for robot kinematic parameter identification.

Appendices
Appendix A Kinematic parameter discernibility analysis for single-stage parallel structure The geometric model of the single-stage parallel structure is shown in Fig. 10(a), and the position vector of end point P T is presented as Eq. (1). Considering the error parameters {δθ 1 , δθ 2 , δd 1 , δd 2 , δr, δx, δy}, the revised single-stage kinematic equation is presented as Eq. (3). Without loss of generality, δθ 1 is taken as an example. The partial derivative with respect to δθ 1 is taken, the six other errors {δθ 2 , δd 1 , δd 2 , δr, δx, δy} are set to be 0 for calculation simplification. Taylor series expansion is applied to the partial derivative at point δθ 1 = 0. When high-order terms are ignored, the constant term coefficient is remarkable enough to reflect the influence rate of parameter δθ 1 in calculating end position vector P T because the identified error parameters are all near 0. As explained in the paper, the constant term is the discernibility of the respective parameter or called  informative value, which is denoted as I(δθ 1 ). The calculation is repeated, and the informative value for all parameters can be acquired, as shown in Table A1. The x and y components of one parameter are shown separately.
Laboratory of Robotics and System (Harbin Institute of Technology, China) (Grant No. SKLRS202009B). No conflicts of interest exist in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution, and reproduction in any medium or format as long as appropriate credit is given to the original author(s) and source, a link to the Creative Commons license is provided, and indicate if changes were made.
The images or other third-party material in this article are included in the ) +