Kinematics for Spacecrafts and Satellites

Chapter 5

Abstract

1. Introduction

87

Kinematic Absolute Positioning

with Quad-Constellation GNSS

Lin Pan, Changsheng Cai, Jianjun Zhu and Xianqiang Cui

The absolute positioning technique is based on a point positioning mode with a single Global Navigation Satellite System (GNSS) receiver, which has been widely used in many fields such as vehicle navigation and kinematic surveying. For a long period, this positioning technique mainly relies on a single GPS system. With the revitalization of Global Navigation Satellite System (GLONASS) constellation and two newly emerging constellations of BeiDou Navigation Satellite System (BDS) and Galileo, it is now feasible to carry out the absolute positioning with quad-constellation of GPS, GLONASS, BDS, and Galileo. A combination of multi-constellation observations can offer improved reliability, availability, and accuracy for position solutions. In this chapter, combined GPS/GLONASS/BDS/Galileo point positioning models for both traditional single point positioning (SPP) and precise point positioning (PPP) are presented, including their functional and stochastic components. The traditional SPP technique has a positioning accuracy at a meter level, whereas the PPP technique can reach an accuracy of a centimeter level. However, the later relies on the availability of precise ephemeris and needs a long convergence time. Experiments were carried out to assess the kinematic positioning performance in the two different modes. The positioning results are compared among different constellation combina-

tions to demonstrate the advantages of quad-constellation GNSS.

Keywords: kinematic positioning, Global Navigation Satellite System,

multi-constellation combination, single point positioning, precise point positioning

Position services have become an inevitable demand for the human activities. Advanced technologies of the position services can significantly improve human's manufacturing efficiency, life quality, and resource utilization. Along with the development of human society, there is an increasing need of kinematic position services, such as automatic drive, intelligent transportation, precision agriculture, and so on. The Global Navigation Satellite System (GNSS), which rose in the 1980s of the last century, is an optimal infrastructure to realize the outdoor kinematic position services. The GNSS-based absolute positioning technologies have many advantages, such as no restriction by the inter-station distance, low cost, and simple data processing. The kinematic positions can be derived globally in all weather and any time based on the GNSS absolute positioning technique with a single receiver. According to different performance demands, two kinds of kinematic absolute positioning technologies can be employed, namely, single point positioning (SPP)

### Chapter 5

## Kinematic Absolute Positioning with Quad-Constellation GNSS

Lin Pan, Changsheng Cai, Jianjun Zhu and Xianqiang Cui

### Abstract

The absolute positioning technique is based on a point positioning mode with a single Global Navigation Satellite System (GNSS) receiver, which has been widely used in many fields such as vehicle navigation and kinematic surveying. For a long period, this positioning technique mainly relies on a single GPS system. With the revitalization of Global Navigation Satellite System (GLONASS) constellation and two newly emerging constellations of BeiDou Navigation Satellite System (BDS) and Galileo, it is now feasible to carry out the absolute positioning with quad-constellation of GPS, GLONASS, BDS, and Galileo. A combination of multi-constellation observations can offer improved reliability, availability, and accuracy for position solutions. In this chapter, combined GPS/GLONASS/BDS/Galileo point positioning models for both traditional single point positioning (SPP) and precise point positioning (PPP) are presented, including their functional and stochastic components. The traditional SPP technique has a positioning accuracy at a meter level, whereas the PPP technique can reach an accuracy of a centimeter level. However, the later relies on the availability of precise ephemeris and needs a long convergence time. Experiments were carried out to assess the kinematic positioning performance in the two different modes. The positioning results are compared among different constellation combinations to demonstrate the advantages of quad-constellation GNSS.

Keywords: kinematic positioning, Global Navigation Satellite System, multi-constellation combination, single point positioning, precise point positioning

#### 1. Introduction

Position services have become an inevitable demand for the human activities. Advanced technologies of the position services can significantly improve human's manufacturing efficiency, life quality, and resource utilization. Along with the development of human society, there is an increasing need of kinematic position services, such as automatic drive, intelligent transportation, precision agriculture, and so on. The Global Navigation Satellite System (GNSS), which rose in the 1980s of the last century, is an optimal infrastructure to realize the outdoor kinematic position services. The GNSS-based absolute positioning technologies have many advantages, such as no restriction by the inter-station distance, low cost, and simple data processing. The kinematic positions can be derived globally in all weather and any time based on the GNSS absolute positioning technique with a single receiver.

According to different performance demands, two kinds of kinematic absolute positioning technologies can be employed, namely, single point positioning (SPP)

and precise point positioning (PPP). The SPP technology can provide meter-level positioning, while the PPP technology has a positioning accuracy at a centimeter level. As satellite-based positioning technologies, the performances of the SPP and PPP are quite dependent on the observed satellites. For a long period, the kinematic absolute positioning technologies are mainly based on a single GPS system. With the recent revitalization of the Global Navigation Satellite System (GLONASS) constellation and two newly emerging constellations of BeiDou Navigation Satellite System (BDS) and Galileo, the quad-constellation integrated absolute positioning has become feasible. Multi-constellation combination is expected to improve the reliability, availability, and accuracy of the SPP and PPP solutions due to the increased measurement redundancy and enhanced satellite geometry, especially when they are performed in areas with GNSS signal blockages.

Following the deployment timeline of BDS, its implementation has been performed in three steps: BeiDou navigation demonstration system (BDS-1) by 2000, regional BDS (BDS-2) by 2012, and global BDS (BDS-3) by 2020. Although the BDS-3 is under construction based on the "three-step" strategy, it has been providing basic services of positioning, navigation, and timing (PNT) for global customers since 27 December 2018. The nominal space constellation of BDS-2 consists of five geostationary earth orbit (GEO) satellites, five inclined geosynchronous orbit (IGSO) satellites, and four medium earth orbit (MEO) satellites. The nominal space constellation of BDS-3 consists of 3 GEO satellites, 3 IGSO satellites, and 24 MEO satellites. The five BDS-2 GEO satellites are located at 58.75°E, 80°E, 110.5°E, 140°E, and 160°E respectively, while the three BDS-3 GEO satellites are placed at 80°E, 110.5°E, and 140°E above the earth's equator, respectively. Actually, spare satellites may be deployed in orbit, according to actual situation. An altitude of 35,786 km with a period of revolution of 23 h 56 min is adopted for the operation of the GEO satellites. The GEO satellites exhibit a non-zero inclination of 0.71.7°, as they are actively controlled in longitudes rather than latitudes. The altitude and orbital period of the IGSO satellites are the same as those of the GEO satellites, and they have an inclination of 55°. Each MEO satellite operates in a nearly circular orbit at an orbit inclination of 55° and an altitude of 21,528 km. The MEO satellites are arranged into three orbital planes, and an angle of 120° is used for the spacing between ascending nodes of different orbital planes. For MEO satellites, the period of the revolution is 12 h 53 min. As of March 2019, there are 15 BDS-2 satellites (5 GEO/7 IGSO/3 MEO), 5 BDS-3 demonstration system (BDS-3S) satellites (2 IGSO/3 MEO), and 19 BDS-3 satellites (1 GEO/18 MEO) in orbit. All BDS-3S satellites and the BDS-3 GEO satellite C59 are in the flight test phase, while the other BDS

Kinematic Absolute Positioning with Quad-Constellation GNSS

DOI: http://dx.doi.org/10.5772/intechopen.86368

When fully deployed, the Galileo constellation will contain 30 satellites in 3 orbital

Figure 1 shows a 24-h ground track of quad-constellation GNSS satellites on 18

March 2019. The quad-constellation mixed precise satellite orbit file is used to derive the satellite coordinates for all GNSS satellites. Different satellites are identified by different colors. The coordinate transformation from reference frame "IGS08" (three-dimensional positions x, y, and z), namely, the realization of international terrestrial reference frame 2008, to geodetic coordinate system (geodetic

planes. There are one inactive spare satellite and nine equally spaced operational satellites in each plane. The ascending nodes of the three planes are equally separated by 120°, and all of them are inclined at an angle of 56°. With a period of about 14 h 7 min and a semimajor axis of 29,600 km, all Galileo satellites are in a nearly circular orbit. The current Galileo space segment is composed of 26 satellites of 2 different generations. Two respective dual launches of four in-orbit validation (IOV) satellites in 2011 and 2012 initiated the buildup of the operational Galileo constellation. A permanent failure of the E5 and E6 signal transmission in May 2014 happened for the IOV-4 satellite due to a sudden power loss. Since then, the IOV-4 satellite can only transmit E1 signal. On 22 August 2014, the first pair of Full Operational Capability (FOC) satellites was launched. However, there was an "orbital injection anomaly" for the two FOC satellites 1 day later, which results in an elliptical orbit with an inclination roughly 5° smaller than planned. Because of the lack of broadcast ephemerides and single-frequency transmission, the IOV-4 satellite currently has the status "not available." As of June 2016, the other three IOV satellites are declared "available," namely, providing broadcast ephemerides and healthy signals that can be used in real-time navigation. Although the two FOC satellites FOC-1 and FOC-2 are not listed in the constellation status and in the eccentric orbit, they are generally transmitting broadcast ephemerides and navigation signals. The other 20 FOC satellites were

successively declared "available" over the past 4 years [4].

satellites are fully operational [3].

89

### 2. Global Navigation Satellite System

A rapid development has been undergone for the satellite-based global navigation systems in recent years. The GNSS family has expanded from a single GPS constellation to four constellations of Galileo, BDS, GLONASS, and GPS. An overview of the four GNSS systems is conducted in terms of their space segment status and navigation signals.

#### 2.1 Space segment status

A nearly circular orbit with an altitude of about 20,200 km is employed for GPS satellites. The GPS satellites pass a same place twice a day. All GPS satellites are located on six equally spaced orbital planes surrounding the earth. For each plane, there are four slots occupied by baseline satellites. The ascending nodes of the orbital planes are equally spaced 60° apart, and they are inclined at 55°. From virtually any point on the earth, users can view at least four GPS satellites, attributing to the 24-slot arrangement. Currently, a 27-slot constellation with improved coverage in most parts of the world is effectively operated for the GPS after expanding the constellation. The United States is committed to maintaining the availability of at least 24 operational GPS satellites 95% of the time. For the past few years, a total of 31 operational GPS satellites have been flying so as to ensure this commitment. The GPS constellation is a mix of old and new satellites, including 1 Block IIA, 11 Block IIR, 7 Block IIR-M, and 12 Block IIF satellites, as of March 2019 [1].

A complete revitalization of GLONASS with a full constellation including 24 operational satellites arranged into 3 orbital planes has been achieved since 2012. A nearly circular orbit is operated for each GLONASS satellite at an altitude of about 19,100 km, and approximately 11 h 16 min is needed for the GLONASS satellites to complete the orbit. The ascending nodes of orbital planes are separated by 120°, and each orbital plane has an inclination angle of 64.8°. The satellites within the same orbital plane are equally spaced 45° apart, while the difference in argument of latitude for satellites in equivalent slots in two different orbital planes is 15°. A continuous global navigation and position service can be provided due to the reasonable spacing of GLONASS satellites. Currently, there are 26 GLONASS satellites in orbit, but only 23 of them are in full operation, including 1 GLONASS-K1 and 22 GLONASS-M satellites. The GLONASS-M satellite SVN 716 is spare, and the SVN 720 satellite of the same series is in maintenance. The GLONASS-K1 satellite SVN 701 K is in the phase of flight tests [2].

#### Kinematic Absolute Positioning with Quad-Constellation GNSS DOI: http://dx.doi.org/10.5772/intechopen.86368

and precise point positioning (PPP). The SPP technology can provide meter-level positioning, while the PPP technology has a positioning accuracy at a centimeter level. As satellite-based positioning technologies, the performances of the SPP and PPP are quite dependent on the observed satellites. For a long period, the kinematic absolute positioning technologies are mainly based on a single GPS system. With the recent revitalization of the Global Navigation Satellite System (GLONASS) constellation and two newly emerging constellations of BeiDou Navigation Satellite System (BDS) and Galileo, the quad-constellation integrated absolute positioning has become feasible. Multi-constellation combination is expected to improve the reliability, availability, and accuracy of the SPP and PPP solutions due to the increased measurement redundancy and enhanced satellite geometry, especially

A rapid development has been undergone for the satellite-based global navigation systems in recent years. The GNSS family has expanded from a single GPS constellation to four constellations of Galileo, BDS, GLONASS, and GPS. An overview of the four GNSS systems is conducted in terms of their space segment status

A nearly circular orbit with an altitude of about 20,200 km is employed for GPS satellites. The GPS satellites pass a same place twice a day. All GPS satellites are located on six equally spaced orbital planes surrounding the earth. For each plane, there are four slots occupied by baseline satellites. The ascending nodes of the orbital planes are equally spaced 60° apart, and they are inclined at 55°. From virtually any point on the earth, users can view at least four GPS satellites, attributing to the 24-slot arrangement. Currently, a 27-slot constellation with improved coverage in most parts of the world is effectively operated for the GPS after expanding the constellation. The United States is committed to maintaining the availability of at least 24 operational GPS satellites 95% of the time. For the past few years, a total of 31 operational GPS satellites have been flying so as to ensure this commitment. The GPS constellation is a mix of old and new satellites, including

1 Block IIA, 11 Block IIR, 7 Block IIR-M, and 12 Block IIF satellites, as of

A complete revitalization of GLONASS with a full constellation including 24 operational satellites arranged into 3 orbital planes has been achieved since 2012. A nearly circular orbit is operated for each GLONASS satellite at an altitude of about 19,100 km, and approximately 11 h 16 min is needed for the GLONASS satellites to complete the orbit. The ascending nodes of orbital planes are separated by 120°, and each orbital plane has an inclination angle of 64.8°. The satellites within the same orbital plane are equally spaced 45° apart, while the difference in argument of latitude for satellites in equivalent slots in two different orbital planes is 15°. A continuous global navigation and position service can be provided due to the reasonable spacing of GLONASS satellites. Currently, there are 26 GLONASS satellites in orbit, but only 23 of them are in full operation, including 1 GLONASS-K1 and 22 GLONASS-M satellites. The GLONASS-M satellite SVN 716 is spare, and the SVN 720 satellite of the same series is in maintenance. The GLONASS-K1 satellite SVN

when they are performed in areas with GNSS signal blockages.

2. Global Navigation Satellite System

Kinematics - Analysis and Applications

and navigation signals.

March 2019 [1].

88

701 K is in the phase of flight tests [2].

2.1 Space segment status

Following the deployment timeline of BDS, its implementation has been performed in three steps: BeiDou navigation demonstration system (BDS-1) by 2000, regional BDS (BDS-2) by 2012, and global BDS (BDS-3) by 2020. Although the BDS-3 is under construction based on the "three-step" strategy, it has been providing basic services of positioning, navigation, and timing (PNT) for global customers since 27 December 2018. The nominal space constellation of BDS-2 consists of five geostationary earth orbit (GEO) satellites, five inclined geosynchronous orbit (IGSO) satellites, and four medium earth orbit (MEO) satellites. The nominal space constellation of BDS-3 consists of 3 GEO satellites, 3 IGSO satellites, and 24 MEO satellites. The five BDS-2 GEO satellites are located at 58.75°E, 80°E, 110.5°E, 140°E, and 160°E respectively, while the three BDS-3 GEO satellites are placed at 80°E, 110.5°E, and 140°E above the earth's equator, respectively. Actually, spare satellites may be deployed in orbit, according to actual situation. An altitude of 35,786 km with a period of revolution of 23 h 56 min is adopted for the operation of the GEO satellites. The GEO satellites exhibit a non-zero inclination of 0.71.7°, as they are actively controlled in longitudes rather than latitudes. The altitude and orbital period of the IGSO satellites are the same as those of the GEO satellites, and they have an inclination of 55°. Each MEO satellite operates in a nearly circular orbit at an orbit inclination of 55° and an altitude of 21,528 km. The MEO satellites are arranged into three orbital planes, and an angle of 120° is used for the spacing between ascending nodes of different orbital planes. For MEO satellites, the period of the revolution is 12 h 53 min. As of March 2019, there are 15 BDS-2 satellites (5 GEO/7 IGSO/3 MEO), 5 BDS-3 demonstration system (BDS-3S) satellites (2 IGSO/3 MEO), and 19 BDS-3 satellites (1 GEO/18 MEO) in orbit. All BDS-3S satellites and the BDS-3 GEO satellite C59 are in the flight test phase, while the other BDS satellites are fully operational [3].

When fully deployed, the Galileo constellation will contain 30 satellites in 3 orbital planes. There are one inactive spare satellite and nine equally spaced operational satellites in each plane. The ascending nodes of the three planes are equally separated by 120°, and all of them are inclined at an angle of 56°. With a period of about 14 h 7 min and a semimajor axis of 29,600 km, all Galileo satellites are in a nearly circular orbit. The current Galileo space segment is composed of 26 satellites of 2 different generations. Two respective dual launches of four in-orbit validation (IOV) satellites in 2011 and 2012 initiated the buildup of the operational Galileo constellation. A permanent failure of the E5 and E6 signal transmission in May 2014 happened for the IOV-4 satellite due to a sudden power loss. Since then, the IOV-4 satellite can only transmit E1 signal. On 22 August 2014, the first pair of Full Operational Capability (FOC) satellites was launched. However, there was an "orbital injection anomaly" for the two FOC satellites 1 day later, which results in an elliptical orbit with an inclination roughly 5° smaller than planned. Because of the lack of broadcast ephemerides and single-frequency transmission, the IOV-4 satellite currently has the status "not available." As of June 2016, the other three IOV satellites are declared "available," namely, providing broadcast ephemerides and healthy signals that can be used in real-time navigation. Although the two FOC satellites FOC-1 and FOC-2 are not listed in the constellation status and in the eccentric orbit, they are generally transmitting broadcast ephemerides and navigation signals. The other 20 FOC satellites were successively declared "available" over the past 4 years [4].

Figure 1 shows a 24-h ground track of quad-constellation GNSS satellites on 18 March 2019. The quad-constellation mixed precise satellite orbit file is used to derive the satellite coordinates for all GNSS satellites. Different satellites are identified by different colors. The coordinate transformation from reference frame "IGS08" (three-dimensional positions x, y, and z), namely, the realization of international terrestrial reference frame 2008, to geodetic coordinate system (geodetic

effects of ionospheric delay based on the dispersive nature of the ionosphere. The modernized GNSS satellites have the capability of transmitting multi-frequency signals. Table 1 details the navigation signals of the four GNSS constellations [6]. k denotes the frequency factor of GLONASS satellites. The GPS Block IIF satellites can provide signals on L1, L2, and L5 frequencies, while the other GPS satellites are still transmitting L1/L2 signals. The last seven satellites of GLONASS-M series and all GLONASS-K satellites operate with three frequency bands, namely G1, G2, and G3, while the other GLONASS satellites can only offer G1 and G2 signals. All Galileo satellites are able to broadcast E1, E5A, E5B, E5 (A + B), and E6 signals. The BDS-2 satellites are capable of providing B1, B2, and B3 triple-frequency signals. In addition to the B1 and B3 signals, the BDS-3 satellites can transmit four new navigation

Frequency (MHz) GPS GLONASS Galileo BDS-2 BDS-3

1575.42 L1 E1 B1C 1561.098 B1 B1

1268.52 B3 B3

1207.14 E5B B2 B2b

1191.795 E5(A + B) B2a + B2b 1176.45 L5 E5A B2a

3. Kinematic single point positioning with quad-constellations

With the use of single-frequency code observations and broadcast ephemeris, the SPP technology can provide meter-level positioning accuracy. Many researchers have focused on error mitigations to improve the SPP performance. The emerging multi-GNSS integration opens new prospects. In this section, the quad-constellation integrated SPP (QISPP) model with GPS, GLONASS, BDS, and Galileo measurements is developed, and its performance in the kinematic mode is evaluated.

Alignment of the coordinate and time references of the four GNSS systems is a key issue for the QISPP. With respect to the coordinate references, the coordinate systems of Galileo, BDS, GLONASS, and GPS satellites adopt the broadcast orbits of GTRF, CGCS2000, PZ90.11, and WGS-84, respectively. Although different coordinate references are employed for the four GNSS systems, the differences among them are only at a level of several centimeters [7, 8]. In view that the code-based positioning solutions using broadcast ephemeris can only achieve an accuracy at a meter level, such a small difference is negligible. In other words, the four GNSS

signals, namely, B1C, B2a, B2b, and B2a + B2b [6].

1602 + k 9/16 G1

DOI: http://dx.doi.org/10.5772/intechopen.86368

1246 + k 7/16 G2

1202.025 G3

1227.60 L2

Navigation signals of quad-constellations.

Table 1.

1278.75 E6

Kinematic Absolute Positioning with Quad-Constellation GNSS

3.1 QISPP model

91

Figure 1. Ground tracks of quad-constellation navigation systems on 18 March 2019.

latitude B, longitude L, and height H) is first performed. Then, the tracks of GNSS satellites are projected to the grounds in the map based on their geodetic coordinates of latitudes and longitudes. The GNSS satellites can provide better PNT services for the stations near their ground tracks. All MEO satellites offer complete global coverage, including Galileo satellites, BDS MEO satellites, GLONASS satellites, and GPS satellites. The ground tracks of the four types of MEO satellites are confined from 57.2°S to 57.2°N latitude, from 56.1°S to 56.1°N latitude, from 65.6°S to 65.6°N latitude, and from 56.8°S to 56.8°N latitude, respectively. Due to the higher inclination angle of GLONASS satellites, they have the most extensive latitude coverage. More satellites at high-latitude areas will be visible by this extensive latitude coverage. The ground tracks of two Galileo FOC satellites in eccentric orbit (FOC-1 and FOC-2), which are represented by the thick lines, show notable asymmetric shape. For the two Galileo satellites, the orbital inclination is 5° lower than nominal. As a result, the peak values of the latitude coverage are reduced. Nevertheless, the scientific, geodetic, and surveying applications can still use the FOC-1 and FOC-2 satellites. The ground tracks of the BDS IGSO satellites are restricted from about 76.2°E to 138.0°E longitude and 57.5°S to 57.5°N latitude. Two figure-ofeight loops can be used to describe the IGSO satellite tracks. The two loops have an average longitude difference of around 30°, so as to effectively cover the western and eastern parts of China as well as the adjoining regions. The intersection points of the two loops are at the longitude of approximately 118°E and 95°E, respectively. The availability of satellites with high elevation angles can be improved by the employment of the inclined geosynchronous orbit. Thus, for users in densely populated areas, the "urban canyon" problems can be alleviated. As shown in Figure 1, the BDS GEO satellites have a movement within a range of 1.8°S–1.8°N latitude, but they are fixed in longitude. The south-north movement can be attributed to the non-zero inclination. To ensure enough visible satellites for users in Asian-Pacific regions, five GEO satellites are distributed in the Indian and Pacific oceans over the equator as supplements for the IGSO satellites [5].

#### 2.2 Navigation signals

The earlier navigation satellites provide signals on two frequencies so that the users can form dual-frequency observation combination to remove first-order

Kinematic Absolute Positioning with Quad-Constellation GNSS DOI: http://dx.doi.org/10.5772/intechopen.86368


#### Table 1.

latitude B, longitude L, and height H) is first performed. Then, the tracks of GNSS satellites are projected to the grounds in the map based on their geodetic coordinates of latitudes and longitudes. The GNSS satellites can provide better PNT services for the stations near their ground tracks. All MEO satellites offer complete global coverage, including Galileo satellites, BDS MEO satellites, GLONASS satellites, and GPS satellites. The ground tracks of the four types of MEO satellites are confined from 57.2°S to 57.2°N latitude, from 56.1°S to 56.1°N latitude, from 65.6°S to 65.6°N latitude, and from 56.8°S to 56.8°N latitude, respectively. Due to the higher inclination angle of GLONASS satellites, they have the most extensive latitude coverage. More satellites at high-latitude areas will be visible by this extensive latitude coverage. The ground tracks of two Galileo FOC satellites in eccentric orbit (FOC-1 and FOC-2), which are represented by the thick lines, show notable asymmetric shape. For the two Galileo satellites, the orbital inclination is 5° lower than nominal. As a result, the peak values of the latitude coverage are reduced. Nevertheless, the scientific, geodetic, and surveying applications can still use the FOC-1 and FOC-2 satellites. The ground tracks of the BDS IGSO satellites are restricted from about 76.2°E to 138.0°E longitude and 57.5°S to 57.5°N latitude. Two figure-ofeight loops can be used to describe the IGSO satellite tracks. The two loops have an average longitude difference of around 30°, so as to effectively cover the western and eastern parts of China as well as the adjoining regions. The intersection points of the two loops are at the longitude of approximately 118°E and 95°E, respectively. The availability of satellites with high elevation angles can be improved by the employment of the inclined geosynchronous orbit. Thus, for users in densely populated areas, the "urban canyon" problems can be alleviated. As shown in Figure 1, the BDS GEO satellites have a movement within a range of 1.8°S–1.8°N latitude, but they are fixed in longitude. The south-north movement can be attributed to the non-zero inclination. To ensure enough visible satellites for users in Asian-Pacific regions, five GEO satellites are distributed in the Indian and Pacific oceans over the

Ground tracks of quad-constellation navigation systems on 18 March 2019.

Kinematics - Analysis and Applications

The earlier navigation satellites provide signals on two frequencies so that the users can form dual-frequency observation combination to remove first-order

equator as supplements for the IGSO satellites [5].

2.2 Navigation signals

90

Figure 1.

Navigation signals of quad-constellations.

effects of ionospheric delay based on the dispersive nature of the ionosphere. The modernized GNSS satellites have the capability of transmitting multi-frequency signals. Table 1 details the navigation signals of the four GNSS constellations [6]. k denotes the frequency factor of GLONASS satellites. The GPS Block IIF satellites can provide signals on L1, L2, and L5 frequencies, while the other GPS satellites are still transmitting L1/L2 signals. The last seven satellites of GLONASS-M series and all GLONASS-K satellites operate with three frequency bands, namely G1, G2, and G3, while the other GLONASS satellites can only offer G1 and G2 signals. All Galileo satellites are able to broadcast E1, E5A, E5B, E5 (A + B), and E6 signals. The BDS-2 satellites are capable of providing B1, B2, and B3 triple-frequency signals. In addition to the B1 and B3 signals, the BDS-3 satellites can transmit four new navigation signals, namely, B1C, B2a, B2b, and B2a + B2b [6].

### 3. Kinematic single point positioning with quad-constellations

With the use of single-frequency code observations and broadcast ephemeris, the SPP technology can provide meter-level positioning accuracy. Many researchers have focused on error mitigations to improve the SPP performance. The emerging multi-GNSS integration opens new prospects. In this section, the quad-constellation integrated SPP (QISPP) model with GPS, GLONASS, BDS, and Galileo measurements is developed, and its performance in the kinematic mode is evaluated.

#### 3.1 QISPP model

Alignment of the coordinate and time references of the four GNSS systems is a key issue for the QISPP. With respect to the coordinate references, the coordinate systems of Galileo, BDS, GLONASS, and GPS satellites adopt the broadcast orbits of GTRF, CGCS2000, PZ90.11, and WGS-84, respectively. Although different coordinate references are employed for the four GNSS systems, the differences among them are only at a level of several centimeters [7, 8]. In view that the code-based positioning solutions using broadcast ephemeris can only achieve an accuracy at a meter level, such a small difference is negligible. In other words, the four GNSS

systems can directly use their satellite coordinates without coordinate transformations in the QISPP. On the other hand, it is not the case for the time scales employed by the four GNSS systems. The GPS Master Control Station establishes the GPS time, which refers to the US Naval Observatory (USNO) Coordinated Universal Time (UTC) with a small difference of <1 μs. Besides, the UTC (USNO) is periodically corrected with integer leap seconds, and thus GPS time differs from it [1]. An atomic time scale UTC (Soviet Union, SU), which is maintained by Russia with an integer difference of 3 h and a fractional difference of <1 ms, is adopted by the GLONASS system [2]. Therefore, in addition to a tiny fractional difference, the GPS time differs from the GLONASS time by leap seconds. The BDS time system (BDT) was synchronized with UTC within 100 ns at 00:00:00 on 1 January 2006, and there exists a constant offset of 14 s between the GPS time and BDT [3]. Apart from a difference of 10 nanoseconds, the Galileo System Time (GST) is nearly identical to the GPS time [4]. The differences among the time references of the four GNSS systems will significantly affect the positioning solutions. Thus, unlike the coordinate reference frames, the inconsistent time scales cannot be ignored and must be properly handled in the QISPP.

using them. Hence, in this study, we only use the single-frequency code observations of each system on the L1/G1/B1/E1 frequencies. The broadcast clock offsets and satellite orbits are referred to the ionosphere-free code combination on two frequencies. Consequently, the hardware delay biases in a form of ionosphere-free combination are contained in the satellite clock offsets derived from the broadcast ephemeris. When the ionosphere-free combined code observables are used, the hardware delay biases can be canceled out for dual-frequency users. But for singlefrequency users, the hardware delay biases must be corrected. Fortunately, the broadcast navigation messages on a satellite-by-satellite basis have provided the time group delays, which can be employed to carry out the hardware delay bias

corrections in the single-frequency pseudorange-based positioning.

AMP <sup>¼</sup> <sup>∑</sup>

8 ><

>:

PER <sup>¼</sup> <sup>∑</sup>

<sup>N</sup>Eð Þ¼ <sup>h</sup> <sup>4</sup>NmE

<sup>1</sup> <sup>þ</sup> exp <sup>h</sup>�hmF<sup>1</sup>

<sup>N</sup>F2ð Þ¼ <sup>h</sup> <sup>4</sup>NmF<sup>2</sup>

<sup>N</sup>F1ð Þ¼ <sup>h</sup> <sup>4</sup>NmF<sup>1</sup>

<sup>1</sup> <sup>þ</sup> exp <sup>h</sup>�hmE

<sup>B</sup><sup>1</sup> <sup>ξ</sup>ð Þ <sup>h</sup> � � � � <sup>2</sup> � exp

<sup>1</sup> <sup>þ</sup> exp <sup>h</sup>�hmF<sup>2</sup>

N hð Þ¼ <sup>4</sup>NmF<sup>2</sup>

B2 � � � � <sup>2</sup> � exp

ð Þ 1 þ exp ð Þz

For more details and specific meaning of the above parameters, refer to Nava

As to the stochastic model for the QISPP, we can use the following covariance

corrections of single-frequency users can be described as

Kinematic Absolute Positioning with Quad-Constellation GNSS

DOI: http://dx.doi.org/10.5772/intechopen.86368

8 ><

>:

details, refer to Klobuchar [12].

equations provided below:

et al. [13].

93

matrix of observations:

dion <sup>¼</sup> <sup>F</sup> � <sup>5</sup>:<sup>0</sup> � <sup>10</sup>�<sup>9</sup> � AMP � <sup>1</sup> � <sup>x</sup><sup>2</sup>

Following Hoque et al. [15], the Klobuchar model for ionospheric error

2 þ x4

<sup>m</sup> AMP≥0

<sup>m</sup> PER≥72, 000

Nbotð Þ¼ h NEð Þþ h NF1ð Þþ h NF2ð Þ h (9)

h � hmE BE <sup>ξ</sup>ð Þ <sup>h</sup>

h � hmF2 B2

h � hmF1 <sup>B</sup><sup>1</sup> <sup>ξ</sup>ð Þ <sup>h</sup>

� � (10)

� � (11)

<sup>2</sup> exp ð Þz (13)

� � (12)

0 AMP<0

72, 000 PER<72, 000

where α and β are coefficients included as part of satellite message, ϕ<sup>m</sup> is the geomagnetic latitude, t is the local time, and E is the elevation angle. For more

The second version of the NeQuick model can be depicted by several main

BE <sup>ξ</sup>ð Þ <sup>h</sup> � � � � <sup>2</sup> � exp

PER (6)

(5)

(7)

(8)

<sup>24</sup> � � � � j j <sup>x</sup> <sup>&</sup>lt;1:<sup>57</sup>

<sup>F</sup> � <sup>5</sup>:<sup>0</sup> � <sup>10</sup>�<sup>9</sup> � � j j <sup>x</sup> <sup>≥</sup>1:<sup>57</sup>

<sup>x</sup> <sup>¼</sup> <sup>2</sup>πð Þ <sup>t</sup> � <sup>50400</sup>

3 n¼0 αnϕ<sup>n</sup>

8 ><

>:

3 n¼0 βnϕ<sup>n</sup>

Although there is only a physical clock in the multi-GNSS receiver, receiver clock parameters with respect to their respective time scales have to be estimated for each satellite system, since different time scales are adopted by the four GNSS systems. Alternatively, instead of adding a receiver clock parameter, a system time difference parameter with respect to a reference time scale can also be introduced [9]. We can directly estimate the GPS receiver clock offset as an unknown parameter, and the receiver clock offsets of the other satellite systems are regarded as a sum of the system time difference parameter and GPS receiver clock, provided that the GPS time scale is chosen as the reference. Following Pan et al. [10], the QISPP observation model reads

$$P^{\mathcal{g}} = \rho^{\mathcal{g}} + cdt^{\mathcal{g}} - cdT^{\mathcal{g}} + d\_{orb}^{\mathcal{g}} + d\_{trvp}^{\mathcal{g}} + d\_{ion}^{\mathcal{g}} + \epsilon\_P^{\mathcal{g}} \tag{1}$$

$$P' = \rho^r + cdt^g + cdt^{r,g}\_{\text{sys}} - cdT^r + d^r\_{orb} + d^r\_{trop} + d^r\_{ion} + e^r\_p \tag{2}$$

$$P^b = \rho^b + cdt^g + cdt\_{sys}^{b,g} - cdT^b + d\_{orb}^b + d\_{trop}^b + d\_{ion}^b + e\_P^b \tag{3}$$

$$P^{\epsilon} = \rho^{\epsilon} + cdt^{\S} + cdt\_{\text{sys}}^{\epsilon, \S} - cdT^{\epsilon} + d\_{orb}^{\epsilon} + d\_{trvp}^{\epsilon} + d\_{ion}^{\epsilon} + \epsilon\_P^{\epsilon} \tag{4}$$

where the superscripts e, b, r, and g refer to Galileo, BDS, GLONASS, and GPS satellites, respectively. P is the measured pseudorange; ρ is the geometric range; c is the speed of light in vacuum; dt<sup>g</sup> is the GPS receiver clock offset; dtr, <sup>g</sup> sys , dtb, <sup>g</sup> sys and dte, <sup>g</sup> sys are the GPS-GLONASS, GPS-BDS, and GPS-Galileo system time differences, respectively; dT is the satellite clock offset; dorb is the satellite orbit error; dtrop is the tropospheric delay; dion is the ionospheric delay; and ɛ is the measurement noise including multipath.

The broadcast ephemeris data is used to compute the clock offset and satellite position, as given in Eqs. (1)–(4). The Saastamoinen model is used to correct the tropospheric delay errors [11]. For the GPS, GLONASS, and BDS systems, the Klobuchar model is used to correct the ionospheric delay errors [12]. The second version of the NeQuick model is employed to perform the ionospheric error correction for the Galileo observations [13]. Regarding the Galileo ionospheric error corrections, the NeQuick model is better suited than the Klobuchar model [14]. Therefore, in the QISPP model, the unknown parameters to be estimated include three system time differences, one GPS receiver clock offset, and three receiver coordinates. Due to the low cost of single-frequency receivers, most SPP users are Kinematic Absolute Positioning with Quad-Constellation GNSS DOI: http://dx.doi.org/10.5772/intechopen.86368

systems can directly use their satellite coordinates without coordinate transformations in the QISPP. On the other hand, it is not the case for the time scales employed by the four GNSS systems. The GPS Master Control Station establishes the GPS time, which refers to the US Naval Observatory (USNO) Coordinated Universal Time (UTC) with a small difference of <1 μs. Besides, the UTC (USNO) is periodically corrected with integer leap seconds, and thus GPS time differs from it [1]. An atomic time scale UTC (Soviet Union, SU), which is maintained by Russia with an integer difference of 3 h and a fractional difference of <1 ms, is adopted by the GLONASS system [2]. Therefore, in addition to a tiny fractional difference, the GPS time differs from the GLONASS time by leap seconds. The BDS time system (BDT) was synchronized with UTC within 100 ns at 00:00:00 on 1 January 2006, and there exists a constant offset of 14 s between the GPS time and BDT [3]. Apart from a difference of 10 nanoseconds, the Galileo System Time (GST) is nearly identical to the GPS time [4]. The differences among the time references of the four GNSS systems will significantly affect the positioning solutions. Thus, unlike the coordinate reference frames, the inconsistent time scales cannot be ignored and must be

Although there is only a physical clock in the multi-GNSS receiver, receiver clock parameters with respect to their respective time scales have to be estimated for each satellite system, since different time scales are adopted by the four GNSS systems. Alternatively, instead of adding a receiver clock parameter, a system time difference parameter with respect to a reference time scale can also be introduced [9]. We can directly estimate the GPS receiver clock offset as an unknown parameter, and the receiver clock offsets of the other satellite systems are regarded as a sum of the system time difference parameter and GPS receiver clock, provided that the GPS time scale is chosen as the reference. Following Pan et al. [10], the QISPP

sys � cdT<sup>r</sup> <sup>þ</sup> dr

sys � cdT<sup>b</sup> <sup>þ</sup> <sup>d</sup><sup>b</sup>

sys � cdT <sup>e</sup> <sup>þ</sup> <sup>d</sup><sup>e</sup>

where the superscripts e, b, r, and g refer to Galileo, BDS, GLONASS, and GPS satellites, respectively. P is the measured pseudorange; ρ is the geometric range; c is

The broadcast ephemeris data is used to compute the clock offset and satellite position, as given in Eqs. (1)–(4). The Saastamoinen model is used to correct the tropospheric delay errors [11]. For the GPS, GLONASS, and BDS systems, the Klobuchar model is used to correct the ionospheric delay errors [12]. The second version of the NeQuick model is employed to perform the ionospheric error correction for the Galileo observations [13]. Regarding the Galileo ionospheric error corrections, the NeQuick model is better suited than the Klobuchar model [14]. Therefore, in the QISPP model, the unknown parameters to be estimated include three system time differences, one GPS receiver clock offset, and three receiver coordinates. Due to the low cost of single-frequency receivers, most SPP users are

orb <sup>þ</sup> <sup>d</sup> <sup>g</sup>

orb <sup>þ</sup> <sup>d</sup><sup>r</sup>

orb <sup>þ</sup> db

orb <sup>þ</sup> de

trop <sup>þ</sup> <sup>d</sup> <sup>g</sup>

trop <sup>þ</sup> dr

trop <sup>þ</sup> <sup>d</sup><sup>b</sup>

trop <sup>þ</sup> <sup>d</sup><sup>e</sup>

ion <sup>þ</sup> <sup>ε</sup> <sup>g</sup>

ion <sup>þ</sup> <sup>ε</sup><sup>r</sup>

ion <sup>þ</sup> <sup>ε</sup><sup>b</sup>

ion <sup>þ</sup> <sup>ε</sup><sup>e</sup>

<sup>P</sup> (1)

<sup>P</sup> (2)

<sup>P</sup> (3)

<sup>P</sup> (4)

sys and dte, <sup>g</sup> sys

sys , dtb, <sup>g</sup>

<sup>P</sup> <sup>g</sup> <sup>¼</sup> <sup>ρ</sup> <sup>g</sup> <sup>þ</sup> cdt <sup>g</sup> � cdT <sup>g</sup> <sup>þ</sup> <sup>d</sup> <sup>g</sup>

the speed of light in vacuum; dt<sup>g</sup> is the GPS receiver clock offset; dtr, <sup>g</sup>

are the GPS-GLONASS, GPS-BDS, and GPS-Galileo system time differences, respectively; dT is the satellite clock offset; dorb is the satellite orbit error; dtrop is the tropospheric delay; dion is the ionospheric delay; and ɛ is the measurement noise

Pr <sup>¼</sup> <sup>ρ</sup><sup>r</sup> <sup>þ</sup> cdt <sup>g</sup> <sup>þ</sup> cdtr, <sup>g</sup>

<sup>P</sup><sup>b</sup> <sup>¼</sup> <sup>ρ</sup><sup>b</sup> <sup>þ</sup> cdt<sup>g</sup> <sup>þ</sup> cdtb, <sup>g</sup>

<sup>P</sup> <sup>e</sup> <sup>¼</sup> <sup>ρ</sup><sup>e</sup> <sup>þ</sup> cdt <sup>g</sup> <sup>þ</sup> cdte, <sup>g</sup>

properly handled in the QISPP.

Kinematics - Analysis and Applications

observation model reads

including multipath.

92

using them. Hence, in this study, we only use the single-frequency code observations of each system on the L1/G1/B1/E1 frequencies. The broadcast clock offsets and satellite orbits are referred to the ionosphere-free code combination on two frequencies. Consequently, the hardware delay biases in a form of ionosphere-free combination are contained in the satellite clock offsets derived from the broadcast ephemeris. When the ionosphere-free combined code observables are used, the hardware delay biases can be canceled out for dual-frequency users. But for singlefrequency users, the hardware delay biases must be corrected. Fortunately, the broadcast navigation messages on a satellite-by-satellite basis have provided the time group delays, which can be employed to carry out the hardware delay bias corrections in the single-frequency pseudorange-based positioning.

Following Hoque et al. [15], the Klobuchar model for ionospheric error corrections of single-frequency users can be described as

$$d\_{\rm ion} = \begin{cases} F \times \left[ 5.0 \times 10^{-9} \times \text{AMP} \times \left( 1 - \frac{\pi^2}{2} + \frac{\pi^4}{24} \right) \right] & |\mathbf{x}| < 1.57\\ F \times \left( 5.0 \times 10^{-9} \right) & |\mathbf{x}| \ge 1.57 \end{cases} \tag{5}$$

$$\infty = \frac{2\pi(t - 50400)}{PER} \tag{6}$$

$$AMP = \begin{cases} \sum\_{n=0}^{3} \alpha\_n \phi\_m^n & AMP \ge 0 \\ 0 & AMP < 0 \end{cases} \tag{7}$$

$$PER = \begin{cases} \sum\_{n=0}^{3} \beta\_n \phi\_m^n & \text{PER} \ge 72,000\\ 72,000 & \text{PER} < 72,000 \end{cases} \tag{8}$$

where α and β are coefficients included as part of satellite message, ϕ<sup>m</sup> is the geomagnetic latitude, t is the local time, and E is the elevation angle. For more details, refer to Klobuchar [12].

The second version of the NeQuick model can be depicted by several main equations provided below:

$$N\_{\rm bot}(h) = N\_{\rm E}(h) + N\_{\rm F1}(h) + N\_{\rm F2}(h) \tag{9}$$

$$N\_{\rm E}(h) = \frac{4NmE}{\left(1 + \exp\left(\frac{h - hmE}{BE}\xi(h)\right)\right)^2} \times \exp\left(\frac{h - hmE}{BE}\xi(h)\right) \tag{10}$$

$$N\_{\rm F1}(h) = \frac{4NmF1}{\left(1 + \exp\left(\frac{h - hmF1}{B1}\xi(h)\right)\right)^2} \times \exp\left(\frac{h - hmF1}{B1}\xi(h)\right) \tag{11}$$

$$N\_{\rm F2}(h) = \frac{4NmF2}{\left(1 + \exp\left(\frac{h - hmF2}{B2}\right)\right)^2} \times \exp\left(\frac{h - hmF2}{B2}\right) \tag{12}$$

$$N(h) = \frac{4NmF2}{\left(1 + \exp\left(z\right)\right)^2} \exp\left(z\right) \tag{13}$$

For more details and specific meaning of the above parameters, refer to Nava et al. [13].

As to the stochastic model for the QISPP, we can use the following covariance matrix of observations:

$$Cov = \begin{bmatrix} Q\_g & 0 & 0 & 0 \\ 0 & Q\_r & 0 & 0 \\ 0 & 0 & Q\_b & 0 \\ 0 & 0 & 0 & Q\_\epsilon \end{bmatrix} \tag{14}$$
 
$$Q = \begin{bmatrix} \sigma\_1^2 & 0 & \cdots & 0 \\ 0 & \sigma\_2^2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \sigma\_n^2 \end{bmatrix} \tag{15}$$

where the subscript n denotes the number of satellites for each satellite system and σ<sup>2</sup> denotes the variance of code observations, which can be written as

$$
\sigma^2 = \sigma\_0^2 / \left(\sin E\right)^2 \tag{16}
$$

why the positioning accuracies can be improved in the multi-constellation

Kinematic positioning errors of four different combination cases for the SPP processing at NNOR.

Kinematic Absolute Positioning with Quad-Constellation GNSS

DOI: http://dx.doi.org/10.5772/intechopen.86368

and 2.96 m in the east, north, and up directions, respectively. In the triple-

and 2% to 0.66, 1.78, and 2.70 m in the three directions, respectively [10].

4. Kinematic precise point positioning with quad-constellations

The PPP technique adopts an absolute positioning approach to achieve centimeter-level positioning accuracy using code and carrier phase observations as well as precise satellite orbit and clock offset corrections [18]. Both industrial

constellation SPP, an accuracy improvement of 9, 6, and 7% over the GPS/GLONASS case to 0.68, 1.82, and 2.75 m in the three directions is achieved, respectively. After a further integration with Galileo, the positioning accuracy is only improved by 3, 2,

The datasets from 47 International GNSS Service (IGS) stations are processed for further analysis, and the average root mean square (RMS) statistics of epoch-wise SPP positioning errors, number of satellites, and PDOP are obtained. The average satellite numbers of the above four processing cases are 8.6, 15.6, 21.1, and 23.4, respectively, while the corresponding average PDOPs are 2.0, 1.4, 1.2, and 1.1, respectively. The combination of GPS and GLONASS improves the positioning accuracy over the GPS-only case by 7, 5, and 5% from 0.81, 2.05, and 3.13 m to 0.75, 1.94,

integrated cases.

Figure 2.

95

where σ<sup>0</sup> denotes the standard deviation (STD) of code observations, which differs among different satellite systems, and E denotes the satellite elevation angle.

The code observation precision is set to 0.3 m for GPS satellites [16]. Following Cai et al. [17], an initial weight ratio of 1:1 is appropriate for BDS and GPS code observations. Thus, the code observation precision is also set to 0.3 m for BDS satellites. The code observation precision is set to 0.6 m for GLONASS satellites due to their twice lower code chipping rate than the GPS code observations. The code observations of Galileo satellites are down-weighted by a factor of four, considering that the broadcast ephemeris has relatively lower accuracies [8]. That is, the precision of the Galileo code observations is also set to 0.6 m.

#### 3.2 Performance analysis of QISPP solutions

Four different constellation combinations are employed for the purpose of comparison, namely, GPS/GLONASS/BDS/Galileo, GPS/GLONASS/BDS, GPS/ GLONASS, and GPS-only. In the data processing, the receiver coordinates as well as other unknown parameters are estimated epoch-by-epoch without imposing any constraints between the epochs in order to analyze the single-epoch SPP performance. For brevity, in the following figures and tables, "GLO" and "GAL" are used to represent GLONASS and Galileo systems, respectively.

Figure 2 shows the epoch-wise SPP positioning errors for the four different combination cases at station NNOR on 8 April 2015. In the east, north, and up directions, the variations of positioning errors are consistent for the four cases, but the series of position errors show less fluctuation for the triple- and quadconstellation cases. It is seen that the GPS-only SPP achieves larger positioning errors than the GPS/GLONASS case. The positioning errors of GPS/GLONASS/BDS SPP at almost all epochs are further reduced by combining with BDS. The further introduction of Galileo observations does not exhibit significant change, since the blue lines are almost completely covered by the orange ones.

Figure 3 presents the position dilution of precision (PDOP) and the number of visible satellites for the four cases. It is clear that the number of visible satellites is obviously increased by the multi-constellation combination and the PDOP value is simultaneously decreased. The quad-constellation integrated case increases the average number of visible satellites from 6.9 to 27.8 in comparison to the GPS-only SPP, and thus the average PDOP values are significantly decreased from 2.3 to 1.1. The increased number of available satellites and decreased PDOP values explain

Cov ¼

Kinematics - Analysis and Applications

Q ¼

sion of the Galileo code observations is also set to 0.6 m.

to represent GLONASS and Galileo systems, respectively.

blue lines are almost completely covered by the orange ones.

94

3.2 Performance analysis of QISPP solutions

σ2

and σ<sup>2</sup> denotes the variance of code observations, which can be written as

<sup>σ</sup><sup>2</sup> <sup>¼</sup> <sup>σ</sup><sup>2</sup>

0 σ<sup>2</sup>

Qg 000 0 Qr 0 0 0 0 Qb 0 000 Qe

<sup>1</sup> 0 ⋯ 0

where the subscript n denotes the number of satellites for each satellite system

where σ<sup>0</sup> denotes the standard deviation (STD) of code observations, which differs among different satellite systems, and E denotes the satellite elevation angle. The code observation precision is set to 0.3 m for GPS satellites [16]. Following Cai et al. [17], an initial weight ratio of 1:1 is appropriate for BDS and GPS code observations. Thus, the code observation precision is also set to 0.3 m for BDS satellites. The code observation precision is set to 0.6 m for GLONASS satellites due to their twice lower code chipping rate than the GPS code observations. The code observations of Galileo satellites are down-weighted by a factor of four, considering that the broadcast ephemeris has relatively lower accuracies [8]. That is, the preci-

Four different constellation combinations are employed for the purpose of com-

GLONASS, and GPS-only. In the data processing, the receiver coordinates as well as other unknown parameters are estimated epoch-by-epoch without imposing any constraints between the epochs in order to analyze the single-epoch SPP performance. For brevity, in the following figures and tables, "GLO" and "GAL" are used

Figure 2 shows the epoch-wise SPP positioning errors for the four different combination cases at station NNOR on 8 April 2015. In the east, north, and up directions, the variations of positioning errors are consistent for the four cases, but

Figure 3 presents the position dilution of precision (PDOP) and the number of visible satellites for the four cases. It is clear that the number of visible satellites is obviously increased by the multi-constellation combination and the PDOP value is simultaneously decreased. The quad-constellation integrated case increases the average number of visible satellites from 6.9 to 27.8 in comparison to the GPS-only SPP, and thus the average PDOP values are significantly decreased from 2.3 to 1.1. The increased number of available satellites and decreased PDOP values explain

parison, namely, GPS/GLONASS/BDS/Galileo, GPS/GLONASS/BDS, GPS/

the series of position errors show less fluctuation for the triple- and quadconstellation cases. It is seen that the GPS-only SPP achieves larger positioning errors than the GPS/GLONASS case. The positioning errors of GPS/GLONASS/BDS SPP at almost all epochs are further reduced by combining with BDS. The further introduction of Galileo observations does not exhibit significant change, since the

<sup>2</sup> ⋯ 0 ⋮ ⋮⋱⋮ 0 0 ⋯ σ<sup>2</sup>

n

<sup>0</sup>=ð Þ sin <sup>E</sup> <sup>2</sup> (16)

(14)

(15)

Figure 2. Kinematic positioning errors of four different combination cases for the SPP processing at NNOR.

why the positioning accuracies can be improved in the multi-constellation integrated cases.

The datasets from 47 International GNSS Service (IGS) stations are processed for further analysis, and the average root mean square (RMS) statistics of epoch-wise SPP positioning errors, number of satellites, and PDOP are obtained. The average satellite numbers of the above four processing cases are 8.6, 15.6, 21.1, and 23.4, respectively, while the corresponding average PDOPs are 2.0, 1.4, 1.2, and 1.1, respectively. The combination of GPS and GLONASS improves the positioning accuracy over the GPS-only case by 7, 5, and 5% from 0.81, 2.05, and 3.13 m to 0.75, 1.94, and 2.96 m in the east, north, and up directions, respectively. In the tripleconstellation SPP, an accuracy improvement of 9, 6, and 7% over the GPS/GLONASS case to 0.68, 1.82, and 2.75 m in the three directions is achieved, respectively. After a further integration with Galileo, the positioning accuracy is only improved by 3, 2, and 2% to 0.66, 1.78, and 2.70 m in the three directions, respectively [10].

### 4. Kinematic precise point positioning with quad-constellations

The PPP technique adopts an absolute positioning approach to achieve centimeter-level positioning accuracy using code and carrier phase observations as well as precise satellite orbit and clock offset corrections [18]. Both industrial

To remove the first-order effects of ionospheric delays, the PPP normally utilizes

2 <sup>2</sup> � P<sup>2</sup> = f

2 <sup>2</sup> � Φ<sup>2</sup> = f

where ΦIF and PIF denote the IF combined carrier phase and code observables,

In order to investigate emerging new satellite systems such as Galileo and BDS, the Multi-GNSS Experiment (MGEX) has been established by the IGS [19]. The correction of satellite clock and orbit errors for QIPPP is conducted with the use of precise satellite clock and orbit products provided by the MGEX. The first-order effects of ionospheric delays are removed using the IF combined observables, as shown in Eqs. (18) and (19). The tropospheric delays can be divided into a wet part and a dry part [20]. The wet part is estimated from the measurements, while the Hopfield tropospheric model is employed to correct the dry part. The projection from slant delays to zenith delays adopts the Niell mapping functions [21]. As some literatures such as Kouba and Héroux [18] have well described other error mitigations, they are not provided here. In PPP, the code-specific hardware delays at the receiver and the receiver clock offsets are usually estimated as a lumped term, as they are linearly correlated with each other. For different navigation systems, both the frequency and signal structures differ. Consequently, within a multi-GNSS receiver, the receiver-dependent code hardware delays are different for the four navigation systems. For the purpose of solving this issue, we should design a receiver clock offset parameter for each satellite system. Alternatively, the differences between receiver clock estimates of different satellite systems can be compensated by introducing an inter-system bias (ISB). The QIPPP observation model can be written as follows, provided that the GPS system is chosen as the

respectively, and f denotes the carrier phase frequency, which differs among

2 <sup>1</sup> � f 2 2

> 2 <sup>1</sup> � f 2 2

(18)

(19)

the ionosphere-free (IF) combined observables, that is:

Kinematic Absolute Positioning with Quad-Constellation GNSS

DOI: http://dx.doi.org/10.5772/intechopen.86368

PIF ¼ f

ΦIF ¼ f

P g

Φ g

Pb

Pr

P e

(ZWD), and M is the tropospheric mapping function.

Φb

Φr

Φe

IF <sup>¼</sup> <sup>ρ</sup> <sup>g</sup> <sup>þ</sup> cdt <sup>þ</sup> <sup>M</sup> <sup>g</sup>

IF <sup>¼</sup> <sup>ρ</sup><sup>b</sup> <sup>þ</sup> cdt <sup>þ</sup> ISBb, <sup>g</sup> <sup>þ</sup> <sup>M</sup>bdzwd <sup>þ</sup> <sup>ε</sup><sup>b</sup>

IF <sup>¼</sup> <sup>ρ</sup><sup>b</sup> <sup>þ</sup> cdt <sup>þ</sup> ISBb, <sup>g</sup> <sup>þ</sup> <sup>M</sup>bdzwd <sup>þ</sup> <sup>N</sup><sup>b</sup>

IF <sup>¼</sup> <sup>ρ</sup><sup>r</sup> <sup>þ</sup> cdt <sup>þ</sup> ISBr, <sup>g</sup> <sup>þ</sup> Mr

IF <sup>¼</sup> <sup>ρ</sup><sup>e</sup> <sup>þ</sup> cdt <sup>þ</sup> ISBe, <sup>g</sup> <sup>þ</sup> <sup>M</sup><sup>e</sup>

where ISBb,g, ISBr,g, and ISBe,g are the GPS-BDS, GPS-GLONASS, and GPS-Galileo inter-system biases, respectively, dzwd is the tropospheric zenith wet delay

For the QIPPP processing, we employ a Kalman filter approach. Actually, the geometric range ρ in Eqs. (20)–(27) is a function of receiver coordinates and satellite coordinates. After the precise ephemeris data is used to determine the satellite coordinates and the ρ is linearized, the unknown parameters include phase ambiguity parameters equal to the number of the observed GNSS satellites, one

IF <sup>¼</sup> <sup>ρ</sup> <sup>g</sup> <sup>þ</sup> cdt <sup>þ</sup> <sup>M</sup><sup>g</sup>

IF <sup>¼</sup> <sup>ρ</sup><sup>r</sup> <sup>þ</sup> cdt <sup>þ</sup> ISBr, <sup>g</sup> <sup>þ</sup> <sup>M</sup><sup>r</sup>

IF <sup>¼</sup> <sup>ρ</sup><sup>e</sup> <sup>þ</sup> cdt <sup>þ</sup> ISBe, <sup>g</sup> <sup>þ</sup> <sup>M</sup><sup>e</sup>

dzwd <sup>þ</sup> <sup>ε</sup> <sup>g</sup>

IF <sup>þ</sup> <sup>ε</sup> <sup>g</sup>

dzwd <sup>þ</sup> <sup>ε</sup><sup>r</sup>

dzwd <sup>þ</sup> <sup>ε</sup><sup>e</sup>

dzwd <sup>þ</sup> <sup>N</sup><sup>r</sup>

dzwd <sup>þ</sup> <sup>N</sup><sup>e</sup>

IF <sup>þ</sup> <sup>ε</sup><sup>b</sup>

IF <sup>þ</sup> <sup>ε</sup><sup>r</sup>

IF <sup>þ</sup> <sup>ε</sup><sup>e</sup>

dzwd <sup>þ</sup> <sup>N</sup> <sup>g</sup>

PIF (20)

<sup>Φ</sup>IF (21)

PIF (22)

PIF (24)

PIF (26)

<sup>Φ</sup>IF (23)

<sup>Φ</sup>IF (25)

<sup>Φ</sup>IF (27)

GLONASS satellites.

reference [22]:

97

2 <sup>1</sup> � P<sup>1</sup> � f

2 <sup>1</sup> � Φ<sup>1</sup> � f

Figure 3. Number of satellites and PDOP values for the SPP processing at NNOR.

applications and scientific research widely use GPS-based PPP. But the position solutions of GPS-based PPP require a long time to converge. Due to high measurement redundancy, significantly reduced convergence time and improved positioning accuracy can be expected by using multi-constellation GNSS PPP. The positioning model and processing method of quad-constellation integrated PPP (QIPPP) with GPS, GLONASS, BDS, and Galileo measurements are developed, and then the improvement in precision and convergence time from quad-constellations is verified by comparing the solutions of different constellation combinations.

#### 4.1 QIPPP model

In view that the code observation equation is detailed in Eqs. (1)–(4), only the carrier phase observations on ith (i = 1, 2) frequency are formulated, that is:

$$
\Phi\_i = \rho + cdt - cdT + d\_{orb} + d\_{trop} - d\_{ion/L\_i} + \lambda\_i N\_i + \varepsilon\_{\Phi\_i} \tag{17}
$$

where Φ is the measured carrier phase, dion/Li is the ionospheric delay on ith frequency, λ is the wavelength, and N is the phase ambiguity term grouped with hardware delays. It should be noted that the wavelength λ is different for different GLONASS satellites because GLONASS employs the frequency-division multiple access (FDMA) technique.

To remove the first-order effects of ionospheric delays, the PPP normally utilizes the ionosphere-free (IF) combined observables, that is:

$$P\_{IF} = \left(f\_1^2 \cdot P\_1 - f\_2^2 \cdot P\_2\right) / \left(f\_1^2 - f\_2^2\right) \tag{18}$$

$$\Phi\_{IF} = \left(f\_1^2 \cdot \Phi\_1 - f\_2^2 \cdot \Phi\_2\right) / \left(f\_1^2 - f\_2^2\right) \tag{19}$$

where ΦIF and PIF denote the IF combined carrier phase and code observables, respectively, and f denotes the carrier phase frequency, which differs among GLONASS satellites.

In order to investigate emerging new satellite systems such as Galileo and BDS, the Multi-GNSS Experiment (MGEX) has been established by the IGS [19]. The correction of satellite clock and orbit errors for QIPPP is conducted with the use of precise satellite clock and orbit products provided by the MGEX. The first-order effects of ionospheric delays are removed using the IF combined observables, as shown in Eqs. (18) and (19). The tropospheric delays can be divided into a wet part and a dry part [20]. The wet part is estimated from the measurements, while the Hopfield tropospheric model is employed to correct the dry part. The projection from slant delays to zenith delays adopts the Niell mapping functions [21]. As some literatures such as Kouba and Héroux [18] have well described other error mitigations, they are not provided here. In PPP, the code-specific hardware delays at the receiver and the receiver clock offsets are usually estimated as a lumped term, as they are linearly correlated with each other. For different navigation systems, both the frequency and signal structures differ. Consequently, within a multi-GNSS receiver, the receiver-dependent code hardware delays are different for the four navigation systems. For the purpose of solving this issue, we should design a receiver clock offset parameter for each satellite system. Alternatively, the differences between receiver clock estimates of different satellite systems can be compensated by introducing an inter-system bias (ISB). The QIPPP observation model can be written as follows, provided that the GPS system is chosen as the reference [22]:

$$P\_{\rm IF}^{\rm g} = \rho^{\rm g} + cdt + \mathcal{M}^{\rm g} d\_{xud} + \varepsilon\_{\rm P\_{\rm IF}}^{\rm g} \tag{20}$$

$$
\Phi\_{\rm IF}^{\rm g} = \rho^{\rm g} + cdt + \mathcal{M}^{\rm g} d\_{xwd} + \mathcal{N}\_{\rm IF}^{\rm g} + \varepsilon\_{\Phi\_{\rm IF}}^{\rm g} \tag{21}
$$

$$P\_{IF}^b = \rho^b + cdt + \text{ISB}\_{b,g} + \mathcal{M}^b d\_{xwd} + \varepsilon\_{P\_{IF}}^b \tag{22}$$

$$
\Phi\_{\rm IF}^{b} = \rho^{b} + cdt + \text{ISB}\_{b,\rm g} + \mathcal{M}^{b}d\_{xwd} + \mathcal{N}\_{\rm IF}^{b} + \varepsilon\_{\Phi\_{\rm IF}}^{b} \tag{23}
$$

$$P\_{\rm IF}^{r} = \rho^{r} + cdt + \text{ISB}\_{r,\rm g} + \mathcal{M}^{r}d\_{xud} + \varepsilon\_{\rm P\_{\rm IF}}^{r} \tag{24}$$

$$
\Phi\_{\rm IF}^{r} = \rho^{r} + cdt + \text{ISB}\_{r,\rm g} + \mathcal{M}^{r}d\_{xvd} + \mathcal{N}\_{\rm IF}^{r} + \varepsilon\_{\Phi\_{\rm IF}}^{r} \tag{25}
$$

$$P\_{\rm IF}^{\epsilon} = \rho^{\epsilon} + cdt + \text{ISB}\_{\epsilon, \text{g}} + \mathcal{M}^{\epsilon} d\_{xud} + \varepsilon\_{P\_{\rm IF}}^{\epsilon} \tag{26}$$

$$
\Phi\_{\rm IF}^{\epsilon} = \rho^{\epsilon} + cdt + \text{ISB}\_{\epsilon, \text{g}} + \mathcal{M}^{\epsilon} d\_{xwd} + \mathcal{N}\_{\rm IF}^{\epsilon} + \varepsilon\_{\Phi\_{\rm IF}}^{\epsilon} \tag{27}
$$

where ISBb,g, ISBr,g, and ISBe,g are the GPS-BDS, GPS-GLONASS, and GPS-Galileo inter-system biases, respectively, dzwd is the tropospheric zenith wet delay (ZWD), and M is the tropospheric mapping function.

For the QIPPP processing, we employ a Kalman filter approach. Actually, the geometric range ρ in Eqs. (20)–(27) is a function of receiver coordinates and satellite coordinates. After the precise ephemeris data is used to determine the satellite coordinates and the ρ is linearized, the unknown parameters include phase ambiguity parameters equal to the number of the observed GNSS satellites, one

applications and scientific research widely use GPS-based PPP. But the position solutions of GPS-based PPP require a long time to converge. Due to high measurement redundancy, significantly reduced convergence time and improved position-

In view that the code observation equation is detailed in Eqs. (1)–(4), only the

where Φ is the measured carrier phase, dion/Li is the ionospheric delay on ith frequency, λ is the wavelength, and N is the phase ambiguity term grouped with hardware delays. It should be noted that the wavelength λ is different for different GLONASS satellites because GLONASS employs the frequency-division multiple

Φ<sup>i</sup> ¼ ρ þ cdt � cdT þ dorb þ dtrop � dion=Li þ λiNi þ ε<sup>Φ</sup><sup>i</sup> (17)

carrier phase observations on ith (i = 1, 2) frequency are formulated, that is:

ing accuracy can be expected by using multi-constellation GNSS PPP. The positioning model and processing method of quad-constellation integrated PPP (QIPPP) with GPS, GLONASS, BDS, and Galileo measurements are developed, and then the improvement in precision and convergence time from quad-constellations is verified by comparing the solutions of different constellation combinations.

Number of satellites and PDOP values for the SPP processing at NNOR.

Kinematics - Analysis and Applications

4.1 QIPPP model

Figure 3.

access (FDMA) technique.

96

ZWD, three ISB, one receiver clock offset, and three receiver coordinates. We should also provide appropriate dynamic models for the state vector and proper stochastic models for the measurements in the Kalman filter. The IF combined phase and code observables are obtained using the raw phase and code measurements on two different frequencies in a form of linear combination. It is assumed that the measurements on different frequencies are independent. Based on the law of random error propagation [23], we can derive the initial variances of the IF combined observables. The satellite elevation angles can be further introduced to acquire the actual variances [24]. As for the dynamic models, the phase ambiguity parameters can be modeled as constants, whereas the ZWD, ISB, receiver clock offset, and kinematic receiver coordinates may be modeled as a random walk (RW) process [25–27].

#### 4.2 Performance analysis of QIPPP solutions

The datasets from stations SEYG, JFNG, and MAR7 on 1 March 2019 are used for numerical analysis in this section. The three stations are located at low-, middle-, and high-latitude regions, respectively, and all of them are able to provide multiconstellation observations. The in-house MIPS-PPP software capable of processing quad-system observation data as well as single-system measurements of Galileo, BDS, GLONASS, and GPS, which is developed at Central South University, China, is employed for the quad-constellation GNSS PPP processing. Regarding the specific PPP position determination, the Galileo E1/E5A, BDS B1/B2, GLONASS G1/G2, and GPS L1/L2 dual-frequency observations are used. The cutoff satellite angle is set to 10°, while the observations are recorded at a sampling rate of 30 s. The IGS Analysis Center, German Research Centre for Geosciences (GFZ), has been generating and releasing the mixed multi-GNSS final precise satellite clock offset and orbit products with a sampling interval of 30 s and 5 min, respectively, and the mitigation of the satellite clock and orbit errors is carried out using them in this study. For the ISB, receiver clock offset, and ZWD parameters, the spectral density values are set to 10<sup>7</sup> , 10<sup>5</sup> , and 10<sup>9</sup> m<sup>2</sup> /s, respectively [28]. The initial STD values for GPS and GLONASS code observations are set to 0.3 and 0.6 m, while for phase observations they are both set to 2 mm [28]. Since there is a relatively lower accuracy for the satellite orbits and clocks of Galileo and BDS [29, 30], their observations are downweighted with a factor of four. That is, for both Galileo and BDS, the code and phase observation accuracies are set to 0.6 m and 4 mm, respectively. Kinematic processing is made on an epoch-by-epoch basis using the static data. No constraints between epochs are imposed so as to simulate kinematic situations. In the Kalman filtering, the coordinates of the dynamic receiver are modeled as a RW process, and the spectral density is set as 10<sup>2</sup> m<sup>2</sup> /s.

Figure 4 shows the positioning errors of GPS, GPS/BDS, GPS/BDS/GLONASS, and four-system PPP in three directions of east, north, and up in the simulated kinematic test. It can be seen from Figure 4 that, compared with the PPP solutions of GPS-only system, the error curve of the PPP of multi-constellation combinations converges to the stable value faster in the east, north, and up directions. For all processing schemes, the positioning errors in the vertical direction are larger than those of horizontal directions.

improved by 11, 55, and 18%, respectively. Compared with the GPS/BDS PPP, the positioning accuracy of the combination of three systems is significantly improved. After further adding Galileo observations, the three-dimensional (3D) accuracy of PPP of the four-system combination is slightly improved by 1.2 cm. Table 3 shows the convergence time in three directions. Taking the JFNG station as an example, the PPP solutions of the four-system combination requires 22.5, 28.5, and 79.0 min

PPP kinematic positioning errors at stations SEYG, JFNG, and MAR7 for four different processing cases.

Kinematic Absolute Positioning with Quad-Constellation GNSS

DOI: http://dx.doi.org/10.5772/intechopen.86368

JFNG East 0.160 0.142 0.052 0.050

MAR7 East 0.066 0.054 0.040 0.032

Stations Directions GPS GPS/BDS GPS/BDS/GLO GPS/BDS/GLO/GAL SEYG East 0.090 0.061 0.058 0.043

> North 0.062 0.057 0.028 0.023 Up 0.158 0.126 0.082 0.071 3D 0.192 0.151 0.104 0.086

> North 0.097 0.044 0.020 0.020 Up 0.192 0.158 0.119 0.106 3D 0.268 0.217 0.131 0.119

> North 0.045 0.035 0.028 0.021 Up 0.081 0.068 0.034 0.030 3D 0.114 0.094 0.060 0.049

to converge to the accuracy level of 1 dm. Compared with the single- and

Figure 4.

Table 2.

99

Convergence accuracy of kinematic PPP (unit: m).

In order to assess the kinematic positioning accuracy, Table 2 provides the RMS statistical values using the position errors in the last 1 h in which the position solutions in all three components have reached stable values. The results show that, taking JFNG station as an example, the positioning accuracy of the GPS-only PPP in three directions is 0.160, 0.097, and 0.192 m, respectively. After the combination of GPS and BDS, compared with the single GPS system, the positioning accuracy is

Kinematic Absolute Positioning with Quad-Constellation GNSS DOI: http://dx.doi.org/10.5772/intechopen.86368

Figure 4.

ZWD, three ISB, one receiver clock offset, and three receiver coordinates. We should also provide appropriate dynamic models for the state vector and proper stochastic models for the measurements in the Kalman filter. The IF combined phase and code observables are obtained using the raw phase and code measurements on two different frequencies in a form of linear combination. It is assumed that the measurements on different frequencies are independent. Based on the law of random error propagation [23], we can derive the initial variances of the IF combined observables. The satellite elevation angles can be further introduced to acquire the actual variances [24]. As for the dynamic models, the phase ambiguity parameters can be modeled as constants, whereas the ZWD, ISB, receiver clock offset, and kinematic receiver coordinates may be modeled as a random walk (RW)

The datasets from stations SEYG, JFNG, and MAR7 on 1 March 2019 are used for numerical analysis in this section. The three stations are located at low-, middle-, and high-latitude regions, respectively, and all of them are able to provide multiconstellation observations. The in-house MIPS-PPP software capable of processing quad-system observation data as well as single-system measurements of Galileo, BDS, GLONASS, and GPS, which is developed at Central South University, China, is employed for the quad-constellation GNSS PPP processing. Regarding the specific PPP position determination, the Galileo E1/E5A, BDS B1/B2, GLONASS G1/G2, and GPS L1/L2 dual-frequency observations are used. The cutoff satellite angle is set to 10°, while the observations are recorded at a sampling rate of 30 s. The IGS Analysis Center, German Research Centre for Geosciences (GFZ), has been generating and releasing the mixed multi-GNSS final precise satellite clock offset and orbit products with a sampling interval of 30 s and 5 min, respectively, and the mitigation of the satellite clock and orbit errors is carried out using them in this study. For the ISB, receiver clock offset, and ZWD parameters, the spectral density values are set

GLONASS code observations are set to 0.3 and 0.6 m, while for phase observations they are both set to 2 mm [28]. Since there is a relatively lower accuracy for the satellite orbits and clocks of Galileo and BDS [29, 30], their observations are downweighted with a factor of four. That is, for both Galileo and BDS, the code and phase

processing is made on an epoch-by-epoch basis using the static data. No constraints between epochs are imposed so as to simulate kinematic situations. In the Kalman filtering, the coordinates of the dynamic receiver are modeled as a RW process, and

Figure 4 shows the positioning errors of GPS, GPS/BDS, GPS/BDS/GLONASS, and four-system PPP in three directions of east, north, and up in the simulated kinematic test. It can be seen from Figure 4 that, compared with the PPP solutions of GPS-only system, the error curve of the PPP of multi-constellation combinations converges to the stable value faster in the east, north, and up directions. For all processing schemes, the positioning errors in the vertical direction are larger than

In order to assess the kinematic positioning accuracy, Table 2 provides the RMS

statistical values using the position errors in the last 1 h in which the position solutions in all three components have reached stable values. The results show that, taking JFNG station as an example, the positioning accuracy of the GPS-only PPP in three directions is 0.160, 0.097, and 0.192 m, respectively. After the combination of GPS and BDS, compared with the single GPS system, the positioning accuracy is

observation accuracies are set to 0.6 m and 4 mm, respectively. Kinematic

/s.

/s, respectively [28]. The initial STD values for GPS and

process [25–27].

to 10<sup>7</sup>

98

, 10<sup>5</sup>

, and 10<sup>9</sup> m<sup>2</sup>

the spectral density is set as 10<sup>2</sup> m<sup>2</sup>

those of horizontal directions.

4.2 Performance analysis of QIPPP solutions

Kinematics - Analysis and Applications

PPP kinematic positioning errors at stations SEYG, JFNG, and MAR7 for four different processing cases.


#### Table 2.

Convergence accuracy of kinematic PPP (unit: m).

improved by 11, 55, and 18%, respectively. Compared with the GPS/BDS PPP, the positioning accuracy of the combination of three systems is significantly improved. After further adding Galileo observations, the three-dimensional (3D) accuracy of PPP of the four-system combination is slightly improved by 1.2 cm. Table 3 shows the convergence time in three directions. Taking the JFNG station as an example, the PPP solutions of the four-system combination requires 22.5, 28.5, and 79.0 min to converge to the accuracy level of 1 dm. Compared with the single- and


List of abbreviations

BDS-2 regional BDS BDS-3 global BDS

BDT BDS time system

BDS BeiDou Navigation Satellite System BDS-1 BeiDou navigation demonstration system

Kinematic Absolute Positioning with Quad-Constellation GNSS

FDMA frequency-division multiple access

GFZ German Research Centre for Geosciences GLONASS Global Navigation Satellite System GNSS Global Navigation Satellite System

BDS-3S BDS-3 demonstration system

FOC full operational capability GEO geostationary earth orbit

DOI: http://dx.doi.org/10.5772/intechopen.86368

GST Galileo System Time IF ionosphere-free

IOV in-orbit validation ISB inter-system bias MEO medium earth orbit MGEX multi-GNSS experiment PDOP position dilution of precision PNT positioning, navigation, and timing

IGS International GNSS Service IGSO inclined geosynchronous orbit

PPP precise point positioning

SPP single point positioning STD standard deviation SU Soviet Union

RMS root mean square RW random walk

ZWD zenith wet delay 3D three-dimensional

Author details

101

QIPPP quad-constellation integrated PPP QISPP quad-constellation integrated SPP

USNO United States Naval Observatory UTC Coordinated Universal Time

Lin Pan, Changsheng Cai\*, Jianjun Zhu and Xianqiang Cui

\*Address all correspondence to: cscai@hotmail.com

provided the original work is properly cited.

School of Geosciences and Info-Physics, Central South University, Changsha, China

© 2019 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium,

#### Table 3.

Convergence time of kinematic PPP (unit: min).

dual-system cases, the convergence time in the horizontal directions of the tripleand quad-system integrated PPP is dramatically shortened. In view of the fact that the positioning accuracy of kinematic PPP in the vertical direction is worse than that of the horizontal directions, the convergence standard may be too strict for the vertical direction, which even leads to the failure of effective convergence in a short period of time sometimes.

### 5. Conclusions

The GNSS-based absolute positioning technologies can provide reliable kinematic position services anywhere, in all weather, and anytime using a single receiver. With single-frequency code measurements and broadcast satellite ephemeris, the SPP technology can provide meter-level positioning accuracy. With dualfrequency code and carrier phase measurements as well as precise satellite orbit and clock products, the PPP technology can offer centimeter-level positioning accuracy. In recent years, the satellite systems have been booming. In view that both SPP and PPP belong to the satellite-based kinematic absolute positioning technologies, the multi-constellation combination provides new prospects for their performance improvement, due to more visible satellites, increased measurement redundancy, and enhanced satellite sky distribution. The quad-constellation integrated SPP and PPP models with GPS, GLONASS, BDS, and Galileo measurements are developed, respectively. The results indicate significantly improved positioning performance of the multi-GNSS integration, which will further promote the applications of SPP and PPP technologies.

#### Acknowledgements

The contribution of data and products from IGS is appreciated.

#### Conflict of interest

The authors declare no conflict of interest.

Kinematic Absolute Positioning with Quad-Constellation GNSS DOI: http://dx.doi.org/10.5772/intechopen.86368

### List of abbreviations


## Author details

dual-system cases, the convergence time in the horizontal directions of the tripleand quad-system integrated PPP is dramatically shortened. In view of the fact that the positioning accuracy of kinematic PPP in the vertical direction is worse than that of the horizontal directions, the convergence standard may be too strict for the vertical direction, which even leads to the failure of effective convergence in a short

Stations Directions GPS GPS/BDS GPS/BDS/GLO GPS/BDS/GLO/GAL

North 37.5 33.0 15.0 14.5 Up 90.0 88.0 77.0 76.5

North 86.0 79.5 28.5 28.5 Up 88.0 82.5 79.0 79.0

North 43.0 41.0 23.5 23.5 Up 86.0 84.5 61.0 60.0

SEYG East 70.5 64.0 51.5 51.5

JFNG East 94.0 85.0 23.0 22.5

MAR7 East 79.0 78.0 68.5 68.5

The GNSS-based absolute positioning technologies can provide reliable kine-

matic position services anywhere, in all weather, and anytime using a single receiver. With single-frequency code measurements and broadcast satellite ephemeris, the SPP technology can provide meter-level positioning accuracy. With dualfrequency code and carrier phase measurements as well as precise satellite orbit and clock products, the PPP technology can offer centimeter-level positioning accuracy. In recent years, the satellite systems have been booming. In view that both SPP and PPP belong to the satellite-based kinematic absolute positioning technologies, the multi-constellation combination provides new prospects for their performance improvement, due to more visible satellites, increased measurement redundancy, and enhanced satellite sky distribution. The quad-constellation integrated SPP and PPP models with GPS, GLONASS, BDS, and Galileo measurements are developed, respectively. The results indicate significantly improved positioning performance of the multi-GNSS integration, which will further promote the applications of SPP and

The contribution of data and products from IGS is appreciated.

The authors declare no conflict of interest.

period of time sometimes.

Convergence time of kinematic PPP (unit: min).

Kinematics - Analysis and Applications

5. Conclusions

Table 3.

PPP technologies.

Acknowledgements

Conflict of interest

100

Lin Pan, Changsheng Cai\*, Jianjun Zhu and Xianqiang Cui School of Geosciences and Info-Physics, Central South University, Changsha, China

\*Address all correspondence to: cscai@hotmail.com

© 2019 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

## References

[1] GPS Directorate. Navstar GPS space segment navigation user segment interfaces, Interface specification (IS-GPS-200). Revision G. Washington: Global Positioning System Directorate; 2012

[2] RISDE. Global Navigation Satellite System GLONASS Interface Control Document. Version 5.1. Moscow: Russian Institute of Space Device Engineering; 2008

[3] CSNO. BeiDou Navigation Satellite System Signal in Space Interface Control Document (Open Service Signal). Version 2.0. Beijing: China Satellite Navigation Office; 2013

[4] EU. European GNSS (Galileo) open service signal in space interface control document (OS-SIS-ICD). Issue 1.1. European Union; 2010

[5] Pan L, Zhang X, Li X, Li X, Lu C, Liu J, et al. Satellite availability and point positioning accuracy evaluation on a global scale for integration of GPS, GLONASS, BeiDou and Galileo. Advances in Space Research. 2019; 63(9):2696-2710. DOI: 10.1016/j.asr. 2017.07.029

[6] RTCM-SC104. The Receiver Independent Exchange Format (RINEX). Version 3.04. International GNSS Service (IGS), RINEX Working Group and Radio Technical Commission for Maritime Services Special Committee 104 (RTCM-SC104). 2018

[7] Torre AD, Caporali A. An analysis of intersystem biases for multi-GNSS positioning. GPS Solutions. 2015;19(2): 297-307

[8] Montenbruck O, Steigenberger P, Hauschild A. Broadcast versus precise ephemerides: A multi-GNSS perspective. GPS Solutions. 2015;19(2): 321-333

[9] Cai C, Gao Y. A combined GPS/ GLONASS navigation algorithm for use with limited satellite visibility. Journal of Navigation. 2009;62(4): 671-685

[17] Cai C, Pan L, Gao Y. A precise weighting approach with application to

positioning. Journal of Navigation. 2014;

DOI: http://dx.doi.org/10.5772/intechopen.86368

Kinematic Absolute Positioning with Quad-Constellation GNSS

Positioning System: Theory and Applications. Progress in Astronautics and Aeronautics. Virginia: American

Institute of Astronautics and Aeronautics; 1996. pp. 409-433

10.1007/BF00873700

s10291-012-0273-9

[28] Cai C, Gao Y. Modeling and

assessment of combined GPS/GLONASS precise point positioning. GPS Solutions. 2013;17(2):223-236. DOI: 10.1007/

[29] Zhao Q, Guo J, Li M, Qu L, Hu Z, Shi C, et al. Initial results of precise orbit and clock determination for COMPASS navigation satellite system. Journal of Geodesy. 2013;87(5):475-486. DOI: 10.1007/s00190-013-0622-7

[30] Steigenberger P, Hugentobler U, Loyer S, Perosanz F, Prange L, Dach R, et al. Galileo orbit and clock quality of the IGS multi-GNSS experiment. Advances in Space Research. 2015;55(1): 269-281. DOI: 10.1016/j.asr.2014.06.030

[26] Brown RG, Hwang PY. Introduction to Random Signals and Applied Kalman Filtering. 3rd ed. New York: Wiley; 1997

[27] Dodson AH, Shardlow PJ, Hubbard LCM, Elgered G, Jarlemark POJ. Wet tropospheric effects on precise relative GPS height determination. Journal of Geodesy. 1996;70(4):188-202. DOI:

[18] Kouba J, Héroux P. Precise point positioning using IGS orbit and clock products. GPS Solutions. 2001;5(2): 12-28. DOI: 10.1007/PL00012883

[19] Rizos C, Montenbruck O, Weber R, Weber G, Neilan R, Hugentobler U. The IGS MGEX experiment as a milestone for a comprehensive multi-GNSS service. In: Proceedings of the ION 2013 Pacific PNT Meeting (ION-PNT-2013); 23-25 April 2013; Honolulu, Hawaii,

[20] Davis JL, Herring TA, Shapiro II, Rogers AEE, Elgered G. Geodesy by radio interferometry: Effects of atmospheric modeling errors on estimates of baseline length. Radio Science. 1985;20(6):1593-1607. DOI:

[21] Niell AE. Global mapping functions for the atmosphere delay at radio wavelengths. Journal of Geophysical Research. 1996;101(B2):3227-3246. DOI:

[22] Cai C, Gao Y, Pan L, Zhu J. Precise

constellations: GPS, BeiDou, GLONASS

[23] Xu G. GPS: Theory, Algorithms and Applications. 2nd ed. Berlin: Springer;

[24] Gerdan GP. A comparison of four methods of weighting double difference pseudorange measurements. The Australian Surveyor. 1995;40(4):60-66

[25] Axelrad P, Brown RG. GPS navigation algorithms. In: Parkinson BW, Spilker JJ, editors. Global

point positioning with quad-

10.1016/j.asr. 2015.04.001

2007

103

and Galileo. Advances in Space Research. 2015;56(1):133-143. DOI:

10.1029/RS020i006p01593

10.1029/95JB03048

combined L1/B1 GPS/BeiDou

67(5):911-925

USA. pp. 289-295

[10] Pan L, Cai C, Santerre R, Zhang X. Performance evaluation of singlefrequency point positioning with GPS, GLONASS, BeiDou and Galileo. Survey Review. 2017;49(354):197-205. DOI: 10.1080/00396265.2016.1151628

[11] Saastamoinen J. Contribution to the theory of atmospheric refraction. Bulletin Géodésique. 1973;107(1):13-34

[12] Klobuchar J. Ionospheric time-delay algorithms for single-frequency GPS users. IEEE Transactions on Aerospace and Electronic Systems. 1987;AES-23 (3):325-331

[13] Nava B, Coïsson P, Radicella SM. A new version of the NeQuick ionosphere electron density model. Journal of Atmospheric and Solar-Terrestrial Physics. 2008;70(15):1856-1862

[14] Oladipo OA, Schüler T. GNSS single frequency ionospheric range delay corrections: NeQuick data ingestion technique. Advances in Space Research. 2012;50(9):1204-1212

[15] Hoque MM, Jakowski N, Berdermann J. An ionosphere broadcast model for next generation GNSS. In: Proceedings of the 28th International Technical Meeting of the ION Satellite Division, ION GNSS+ 2015; 14-18 September 2015; Tampa, Florida, USA. pp. 3755-3765

[16] Cai C, Gao Y, Pan L, Dai W. An analysis on combined GPS/COMPASS data quality and its effect on single point positioning accuracy under different observing conditions. Advances in Space Research. 2014;54(5):818-829

Kinematic Absolute Positioning with Quad-Constellation GNSS DOI: http://dx.doi.org/10.5772/intechopen.86368

[17] Cai C, Pan L, Gao Y. A precise weighting approach with application to combined L1/B1 GPS/BeiDou positioning. Journal of Navigation. 2014; 67(5):911-925

References

Engineering; 2008

Navigation Office; 2013

European Union; 2010

2017.07.029

297-307

321-333

102

2012

[1] GPS Directorate. Navstar GPS space segment navigation user segment interfaces, Interface specification (IS-GPS-200). Revision G. Washington: Global Positioning System Directorate;

Kinematics - Analysis and Applications

[9] Cai C, Gao Y. A combined GPS/ GLONASS navigation algorithm for use with limited satellite visibility. Journal of Navigation. 2009;62(4):

[10] Pan L, Cai C, Santerre R, Zhang X. Performance evaluation of singlefrequency point positioning with GPS, GLONASS, BeiDou and Galileo. Survey Review. 2017;49(354):197-205. DOI: 10.1080/00396265.2016.1151628

[11] Saastamoinen J. Contribution to the theory of atmospheric refraction. Bulletin Géodésique. 1973;107(1):13-34

[12] Klobuchar J. Ionospheric time-delay algorithms for single-frequency GPS users. IEEE Transactions on Aerospace and Electronic Systems. 1987;AES-23

[13] Nava B, Coïsson P, Radicella SM. A new version of the NeQuick ionosphere electron density model. Journal of Atmospheric and Solar-Terrestrial Physics. 2008;70(15):1856-1862

[14] Oladipo OA, Schüler T. GNSS single frequency ionospheric range delay corrections: NeQuick data ingestion technique. Advances in Space Research.

Berdermann J. An ionosphere broadcast model for next generation GNSS. In: Proceedings of the 28th International Technical Meeting of the ION Satellite Division, ION GNSS+ 2015; 14-18 September 2015; Tampa, Florida, USA.

[16] Cai C, Gao Y, Pan L, Dai W. An analysis on combined GPS/COMPASS data quality and its effect on single point positioning accuracy under different observing conditions. Advances in Space

Research. 2014;54(5):818-829

2012;50(9):1204-1212

pp. 3755-3765

[15] Hoque MM, Jakowski N,

671-685

(3):325-331

[2] RISDE. Global Navigation Satellite System GLONASS Interface Control Document. Version 5.1. Moscow: Russian Institute of Space Device

[3] CSNO. BeiDou Navigation Satellite System Signal in Space Interface Control Document (Open Service Signal). Version 2.0. Beijing: China Satellite

[4] EU. European GNSS (Galileo) open service signal in space interface control document (OS-SIS-ICD). Issue 1.1.

[5] Pan L, Zhang X, Li X, Li X, Lu C, Liu J, et al. Satellite availability and point positioning accuracy evaluation on a global scale for integration of GPS, GLONASS, BeiDou and Galileo. Advances in Space Research. 2019; 63(9):2696-2710. DOI: 10.1016/j.asr.

[6] RTCM-SC104. The Receiver Independent Exchange Format (RINEX). Version 3.04. International GNSS Service (IGS), RINEX Working Group and Radio Technical Commission

for Maritime Services Special

Committee 104 (RTCM-SC104). 2018

[7] Torre AD, Caporali A. An analysis of intersystem biases for multi-GNSS positioning. GPS Solutions. 2015;19(2):

[8] Montenbruck O, Steigenberger P, Hauschild A. Broadcast versus precise

perspective. GPS Solutions. 2015;19(2):

ephemerides: A multi-GNSS

[18] Kouba J, Héroux P. Precise point positioning using IGS orbit and clock products. GPS Solutions. 2001;5(2): 12-28. DOI: 10.1007/PL00012883

[19] Rizos C, Montenbruck O, Weber R, Weber G, Neilan R, Hugentobler U. The IGS MGEX experiment as a milestone for a comprehensive multi-GNSS service. In: Proceedings of the ION 2013 Pacific PNT Meeting (ION-PNT-2013); 23-25 April 2013; Honolulu, Hawaii, USA. pp. 289-295

[20] Davis JL, Herring TA, Shapiro II, Rogers AEE, Elgered G. Geodesy by radio interferometry: Effects of atmospheric modeling errors on estimates of baseline length. Radio Science. 1985;20(6):1593-1607. DOI: 10.1029/RS020i006p01593

[21] Niell AE. Global mapping functions for the atmosphere delay at radio wavelengths. Journal of Geophysical Research. 1996;101(B2):3227-3246. DOI: 10.1029/95JB03048

[22] Cai C, Gao Y, Pan L, Zhu J. Precise point positioning with quadconstellations: GPS, BeiDou, GLONASS and Galileo. Advances in Space Research. 2015;56(1):133-143. DOI: 10.1016/j.asr. 2015.04.001

[23] Xu G. GPS: Theory, Algorithms and Applications. 2nd ed. Berlin: Springer; 2007

[24] Gerdan GP. A comparison of four methods of weighting double difference pseudorange measurements. The Australian Surveyor. 1995;40(4):60-66

[25] Axelrad P, Brown RG. GPS navigation algorithms. In: Parkinson BW, Spilker JJ, editors. Global

Positioning System: Theory and Applications. Progress in Astronautics and Aeronautics. Virginia: American Institute of Astronautics and Aeronautics; 1996. pp. 409-433

[26] Brown RG, Hwang PY. Introduction to Random Signals and Applied Kalman Filtering. 3rd ed. New York: Wiley; 1997

[27] Dodson AH, Shardlow PJ, Hubbard LCM, Elgered G, Jarlemark POJ. Wet tropospheric effects on precise relative GPS height determination. Journal of Geodesy. 1996;70(4):188-202. DOI: 10.1007/BF00873700

[28] Cai C, Gao Y. Modeling and assessment of combined GPS/GLONASS precise point positioning. GPS Solutions. 2013;17(2):223-236. DOI: 10.1007/ s10291-012-0273-9

[29] Zhao Q, Guo J, Li M, Qu L, Hu Z, Shi C, et al. Initial results of precise orbit and clock determination for COMPASS navigation satellite system. Journal of Geodesy. 2013;87(5):475-486. DOI: 10.1007/s00190-013-0622-7

[30] Steigenberger P, Hugentobler U, Loyer S, Perosanz F, Prange L, Dach R, et al. Galileo orbit and clock quality of the IGS multi-GNSS experiment. Advances in Space Research. 2015;55(1): 269-281. DOI: 10.1016/j.asr.2014.06.030

Chapter 6

Abstract

and inverse kinematics.

attitude, orbit

1. Introduction

105

Kinematics for Spacecraft-Type

The scope of this chapter is the study of the forward and inverse kinematics for a space robot. The main focus is to compute the position and orientation of manipulators' end-effectors relative to their platform. Such platform plays the role of workstations referred in the literature approaching ground manipulators. In this study, the method is to write the manipulator kinematics' equations as functions of the joint variables by following the Denavit-Hartenberg convention. The homogeneous transform technique is used to study the kinematics. The set of coordinate frames defined in this chapter follows the convention for frames that appears in the literature for ground robot manipulators. The kinematics related to the spacecraft attitude is added in the formulation because the manipulator studied in this chapter is type spacecraft. The objective is to provide an overview and clear understanding of the kinematics' equations for spacecraft-type manipulators. To be consistent with orbital dynamics area, the inertial, orbital, and body-fixed coordinate frames are included in this kinematics study. The forward and inverse kinematics formulations are derived. The MATLAB®/Simulink tools are presented for the computer simulations of the forward

Keywords: space robot manipulators, forward kinematics, inverse kinematics,

Kinematics is mainly concerned with the geometry of motion. The subject of kinematics is to some extent mathematical and does not consider any forces associated with the motion. A comprehensive study of kinematics is a prerequisite to the successful formulation of the equations of motion and the dynamic analysis. Consistent with the geometry of motion, kinematics involves the definition of systems of reference, methods of establishing relationships between frames, and algebraic manipulations of matrices and vectors. Kinematics is primarily concerned with describing the orientation of a body with respect to a known reference or coordinate frame. When dealing with robotic manipulators, we have to include also translational kinematics. In this chapter, the calculation goes beyond the algebraic manipulations of rotation matrices connecting frames. This includes the derivation of the equations for the rotational and the translational kinematics that characterize robot manipulators' kinematics. Space robotic manipulators are designed to implement orbital operations. There is a considerable difference regarding Earth-based manipulators, as they increase the type and number of reference systems for proper

Ijar Milagre da Fonseca, Maurício Nacib Pontuschka

Robotic Manipulators

and Glaydson Luiz Bertoze Lima

### Chapter 6

## Kinematics for Spacecraft-Type Robotic Manipulators

Ijar Milagre da Fonseca, Maurício Nacib Pontuschka and Glaydson Luiz Bertoze Lima

### Abstract

The scope of this chapter is the study of the forward and inverse kinematics for a space robot. The main focus is to compute the position and orientation of manipulators' end-effectors relative to their platform. Such platform plays the role of workstations referred in the literature approaching ground manipulators. In this study, the method is to write the manipulator kinematics' equations as functions of the joint variables by following the Denavit-Hartenberg convention. The homogeneous transform technique is used to study the kinematics. The set of coordinate frames defined in this chapter follows the convention for frames that appears in the literature for ground robot manipulators. The kinematics related to the spacecraft attitude is added in the formulation because the manipulator studied in this chapter is type spacecraft. The objective is to provide an overview and clear understanding of the kinematics' equations for spacecraft-type manipulators. To be consistent with orbital dynamics area, the inertial, orbital, and body-fixed coordinate frames are included in this kinematics study. The forward and inverse kinematics formulations are derived. The MATLAB®/Simulink tools are presented for the computer simulations of the forward and inverse kinematics.

Keywords: space robot manipulators, forward kinematics, inverse kinematics, attitude, orbit

#### 1. Introduction

Kinematics is mainly concerned with the geometry of motion. The subject of kinematics is to some extent mathematical and does not consider any forces associated with the motion. A comprehensive study of kinematics is a prerequisite to the successful formulation of the equations of motion and the dynamic analysis. Consistent with the geometry of motion, kinematics involves the definition of systems of reference, methods of establishing relationships between frames, and algebraic manipulations of matrices and vectors. Kinematics is primarily concerned with describing the orientation of a body with respect to a known reference or coordinate frame. When dealing with robotic manipulators, we have to include also translational kinematics. In this chapter, the calculation goes beyond the algebraic manipulations of rotation matrices connecting frames. This includes the derivation of the equations for the rotational and the translational kinematics that characterize robot manipulators' kinematics. Space robotic manipulators are designed to implement orbital operations. There is a considerable difference regarding Earth-based manipulators, as they increase the type and number of reference systems for proper

problem formulation [1–4]. The main differences with regard to Earth-based robotic manipulators are the inclusion of an Earth-centered inertial coordinate frame (ECI), besides the local vertical, local horizontal (LVLH) reference frame, and the spacecraft fixed reference frame. These frames are used to describe the kinematics' differential equations related to orbit and attitude motions. The other reference frames are typical of those defined for ground manipulators. For space manipulators, we need to consider the relationships between all these reference systems.

### 2. Forward kinematics of manipulators and D-H convention homogeneous transforms

Forward kinematics is the problem of computing the position and orientation of the robot tool with respect to a fixed reference frame defined on the base of the space robot. Such manipulators can be seen as a set of links connected in a chain of joints. The joints are points of connection between a pair of links. The joints of a manipulator may be revolute (rotatory), prismatic (sliding), or a combination of both. According to that link to joint configuration, any robot manipulator kinematics can be described in terms of four quantities associated with each link. Two of these quantities describe the link itself. The other two describe the link's connection to the neighboring link. For revolute joints, the angle θ<sup>i</sup> is the joint variable, and the other three quantities are fixed link parameters. These parameters are illustrated in Figure 1 as the twist angle α<sup>i</sup>�<sup>1</sup>,the distances ai�1, and di. For prismatic joints, the joint variable is di, while the angle θ<sup>i</sup> is zero. In this case, α<sup>i</sup>�<sup>1</sup>, and ai�<sup>1</sup> are fixed link parameters (see Figure 1). The description of manipulators by these quantities is referred as the Denavit-Hartenberg convention, denoted here as D-H [5, 6]. The understanding of those four quantities requires the frame convention shown in Figure 2. The frames are denoted by the symbol {}, for example {i-1} denotes a frame whose origin is in joint i-1. A summary of the link parameters considering frames convention shown in Figure 2 can be written as:

The D-H convention describes the manipulator rotational plus translational kinematics by using homogeneous transform. This transform is a 4 � 4 matrix in

position of the end-effector is given with respect to the origin of the frame <sup>1</sup>�<sup>1</sup>      <sup>i</sup>R:

The 3 � <sup>3</sup> <sup>1</sup>�<sup>1</sup>      <sup>i</sup><sup>R</sup> block matrix is related to the rotation between frames f g <sup>i</sup> � <sup>1</sup> and f g<sup>i</sup> by the joint angle <sup>θ</sup>i. The 3 � 1 vector <sup>i</sup>�<sup>1</sup>po is related to the end-effector position with respect o the frame f g i � 1 . The advantage of this approach is the use of matrix

Such matrix notation facilitates the mathematical formulation in a compact and comprehensive format. Also, the D-H formulation is appropriate to implement

Let us discuss now the manipulation of the frames as defined in Figure 2 and their parameters. The D-H convention considers that all the joint angle rotations take place along z-axis. The rule is the same for the joint variables di (prismatic joints). In the case of prismatic joints, the joint angle is zero. This convention regarding joint variables along z-axis means that when z-axis is not parallel to the next one, a twist angle rotation α must be done so as to generate a new z-axis parallel to the next one. Angle α can be positive or negative. Another convention applies to the distance ai. Once we have the new z-axis parallel to the next one, the distance ai is determined by drawing a mutual perpendicular line between both zaxes. Such line will determine the distance between the joints and is also the place of a new x-axis. A y-axis completes the right-handed frame (Figure 2). If the line containing ai crosses the z-axis offset of the joint, it defines a distance di as shown in Figure 2. Finally, we rotate the z-axis by the joint variable θ<sup>i</sup> and a new frame is

On the base of the space robot attachment, the Ox0y0z<sup>0</sup> frame or 0f g is defined. Note that the frame is f g i � 1 , and for i = 1, this frame is 0f g. In the literature of robot manipulators, 0f g frame does not move with respect to the space robot platform. Accordingly, the joint associated to {0} is also named joint 0 from which the link 0

i po 1 � �

¼

<sup>p</sup>þ<sup>i</sup>�<sup>1</sup>po and 1 = 1. The subscript <sup>o</sup> is just to say that the

<sup>1</sup>�<sup>1</sup> <sup>i</sup>R <sup>i</sup>�<sup>1</sup>po 0 1 " # <sup>i</sup>

p 1 � � (1)

cos <sup>θ</sup><sup>i</sup> � sin <sup>θ</sup><sup>i</sup> <sup>0</sup> <sup>i</sup>�<sup>1</sup>px sin θ<sup>i</sup> cos θ<sup>i</sup> 0 <sup>i</sup>�<sup>1</sup>py 0 01 <sup>i</sup>�<sup>1</sup>pz 0 0 01

notation including the position in addition to the orientation.

problems in the MATLAB® computational environment.

the form:

107

Figure 2.

<sup>i</sup>�<sup>1</sup>p 1 � �

¼

The result is <sup>i</sup>�<sup>1</sup>

Frames, links, and D-H quantities.

    <sup>i</sup> <sup>p</sup> <sup>¼</sup> <sup>i</sup>�<sup>1</sup>     <sup>i</sup> <sup>R</sup><sup>i</sup>

Kinematics for Spacecraft-Type Robotic Manipulators DOI: http://dx.doi.org/10.5772/intechopen.85636

generated, taking into account the offset di.


Figure 1. D-H convention showing the notation, joint variable, and parameters.

Figure 2. Frames, links, and D-H quantities.

problem formulation [1–4]. The main differences with regard to Earth-based robotic manipulators are the inclusion of an Earth-centered inertial coordinate frame (ECI), besides the local vertical, local horizontal (LVLH) reference frame, and the spacecraft fixed reference frame. These frames are used to describe the kinematics' differential equations related to orbit and attitude motions. The other reference frames are typical of those defined for ground manipulators. For space manipulators, we need to con-

Forward kinematics is the problem of computing the position and orientation of the robot tool with respect to a fixed reference frame defined on the base of the space robot. Such manipulators can be seen as a set of links connected in a chain of joints. The joints are points of connection between a pair of links. The joints of a manipulator may be revolute (rotatory), prismatic (sliding), or a combination of both. According to that link to joint configuration, any robot manipulator kinematics can be described in terms of four quantities associated with each link. Two of these quantities describe the link itself. The other two describe the link's connection to the neighboring link. For revolute joints, the angle θ<sup>i</sup> is the joint variable, and the other three quantities are fixed link parameters. These parameters are illustrated in Figure 1 as the twist angle α<sup>i</sup>�<sup>1</sup>,the distances ai�1, and di. For prismatic joints, the joint variable is di, while the angle θ<sup>i</sup> is zero. In this case, α<sup>i</sup>�<sup>1</sup>, and ai�<sup>1</sup> are fixed link parameters (see Figure 1). The description of manipulators by these quantities is referred as the Denavit-Hartenberg convention, denoted here as D-H [5, 6]. The understanding of those four quantities requires the frame convention shown in Figure 2. The frames are denoted by the symbol {}, for example {i-1} denotes a frame whose origin is in joint i-1. A summary of the link parameters considering

sider the relationships between all these reference systems.

frames convention shown in Figure 2 can be written as:

• ai is the distance measured along x-axis from zi to ziþ<sup>1</sup>

• α<sup>i</sup> is the angle rotation measured about xi from zi to ziþ<sup>1</sup>

• di is the distance measured along zi from xi�<sup>1</sup> to xi

• θ<sup>i</sup> is a joint angle measured about zi from xi�<sup>1</sup> to xi.

D-H convention showing the notation, joint variable, and parameters.

Figure 1.

106

homogeneous transforms

Kinematics - Analysis and Applications

2. Forward kinematics of manipulators and D-H convention—

The D-H convention describes the manipulator rotational plus translational kinematics by using homogeneous transform. This transform is a 4 � 4 matrix in the form:

$$
\begin{Bmatrix} i^{-1}\mathbf{p} \\ \mathbf{1} \end{Bmatrix} = \begin{bmatrix} \cos\theta\_i & -\sin\theta\_i & \mathbf{0} & ^{i-1}\mathbf{p}\_x \\ \sin\theta\_i & \cos\theta\_i & \mathbf{0} & ^{i-1}\mathbf{p}\_y \\ \mathbf{0} & \mathbf{0} & \mathbf{1} & ^{i-1}\mathbf{p}\_x \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{1} \end{bmatrix} \begin{Bmatrix} i\mathbf{p}\_o \\ \mathbf{1} \end{Bmatrix} = \begin{bmatrix} \mathbf{1}^{-1}\mathbf{R} & i^{-1}\mathbf{p}\_o \\ \mathbf{0} & \mathbf{1} \end{bmatrix} \begin{Bmatrix} i\mathbf{p} \\ \mathbf{1} \end{Bmatrix} \tag{1}
$$

The result is <sup>i</sup>�<sup>1</sup>     <sup>i</sup> <sup>p</sup> <sup>¼</sup> <sup>i</sup>�<sup>1</sup>     <sup>i</sup> <sup>R</sup><sup>i</sup> <sup>p</sup>þ<sup>i</sup>�<sup>1</sup>po and 1 = 1. The subscript <sup>o</sup> is just to say that the position of the end-effector is given with respect to the origin of the frame <sup>1</sup>�<sup>1</sup>      <sup>i</sup>R:

The 3 � <sup>3</sup> <sup>1</sup>�<sup>1</sup>      <sup>i</sup><sup>R</sup> block matrix is related to the rotation between frames f g <sup>i</sup> � <sup>1</sup> and f g<sup>i</sup> by the joint angle <sup>θ</sup>i. The 3 � 1 vector <sup>i</sup>�<sup>1</sup>po is related to the end-effector position with respect o the frame f g i � 1 . The advantage of this approach is the use of matrix notation including the position in addition to the orientation.

Such matrix notation facilitates the mathematical formulation in a compact and comprehensive format. Also, the D-H formulation is appropriate to implement problems in the MATLAB® computational environment.

Let us discuss now the manipulation of the frames as defined in Figure 2 and their parameters. The D-H convention considers that all the joint angle rotations take place along z-axis. The rule is the same for the joint variables di (prismatic joints). In the case of prismatic joints, the joint angle is zero. This convention regarding joint variables along z-axis means that when z-axis is not parallel to the next one, a twist angle rotation α must be done so as to generate a new z-axis parallel to the next one. Angle α can be positive or negative. Another convention applies to the distance ai. Once we have the new z-axis parallel to the next one, the distance ai is determined by drawing a mutual perpendicular line between both zaxes. Such line will determine the distance between the joints and is also the place of a new x-axis. A y-axis completes the right-handed frame (Figure 2). If the line containing ai crosses the z-axis offset of the joint, it defines a distance di as shown in Figure 2. Finally, we rotate the z-axis by the joint variable θ<sup>i</sup> and a new frame is generated, taking into account the offset di.

On the base of the space robot attachment, the Ox0y0z<sup>0</sup> frame or 0f g is defined. Note that the frame is f g i � 1 , and for i = 1, this frame is 0f g. In the literature of robot manipulators, 0f g frame does not move with respect to the space robot platform. Accordingly, the joint associated to {0} is also named joint 0 from which the link 0

connects to the joint 1. It is common to define joint 0 coincident with joint 1. In this case, a<sup>0</sup> is zero. For ground robots, it is common to define Ox0y0z<sup>0</sup> as an inertial frame. However, such definition is not consistent with the concept of inertial frames, mainly when we have space robot orbiting the Earth. For a space robotic manipulator, the Ox0y0z<sup>0</sup> frame is defined with its origin somewhere on the platform and locates the placement of joint zero. Figures 2 and 3 illustrate the frames, joint axes, joint variables, and the other fixed quantities according to the D-H convention.

i�1

i�1      i T ¼

i�1     <sup>i</sup> <sup>T</sup> <sup>¼</sup>

109

<sup>P</sup> <sup>¼</sup> <sup>i</sup>�<sup>1</sup>    <sup>R</sup> <sup>T</sup><sup>R</sup>

where <sup>i</sup>�<sup>1</sup>

    <sup>i</sup> <sup>T</sup> <sup>¼</sup> <sup>i</sup>�<sup>1</sup>     RT<sup>R</sup>

In Eq. (2), the first matrix,

¼

<sup>Q</sup> <sup>T</sup><sup>Q</sup> <sup>P</sup> T<sup>P</sup> <sup>i</sup> T    <sup>i</sup>

10 0 0 0 cα<sup>i</sup>�<sup>1</sup> �sα<sup>i</sup>�<sup>1</sup> 0 0 sα<sup>i</sup>�<sup>1</sup> cα<sup>i</sup>�<sup>1</sup> 0 00 0 1

10 0 ai�<sup>1</sup> 0 cα<sup>i</sup>�<sup>1</sup> �sα<sup>i</sup>�<sup>1</sup> 0 0 sα<sup>i</sup>�<sup>1</sup> cα<sup>i</sup>�<sup>1</sup> 0 00 0 1

<sup>P</sup> <sup>¼</sup> <sup>i</sup>�<sup>1</sup> <sup>i</sup> <sup>T</sup>    <sup>i</sup>

Kinematics for Spacecraft-Type Robotic Manipulators DOI: http://dx.doi.org/10.5772/intechopen.85636

> cθ<sup>i</sup> �sθ<sup>i</sup> 0 ai�<sup>1</sup> sθicα<sup>i</sup>�<sup>1</sup> cθicα<sup>i</sup>�<sup>1</sup> �sα<sup>i</sup>�<sup>1</sup> �sα<sup>i</sup>�<sup>1</sup>di sθisα<sup>i</sup>�<sup>1</sup> cθisα<sup>i</sup>�<sup>1</sup> cα<sup>i</sup>�<sup>1</sup> cα<sup>i</sup>�<sup>1</sup>di 0 00 1

> > <sup>Q</sup> <sup>T</sup><sup>Q</sup> <sup>P</sup> T<sup>P</sup>

from the previous to the next coordinate frame.

c<sup>123</sup> ¼ c ∑

3 i θi

both the previous and the next z-axis. The rotation matrix,

� � and <sup>s</sup><sup>123</sup> <sup>¼</sup> <sup>s</sup> <sup>∑</sup>

10 0 0 0 cα<sup>i</sup>�<sup>1</sup> �sα<sup>i</sup>�<sup>1</sup> 0 0 sα<sup>i</sup>�<sup>1</sup> cα<sup>i</sup>�<sup>1</sup> 0 00 0 1

is given by a rotation of α<sup>i</sup>�<sup>1</sup> about xi�1. Note that no translation appears in the matrix. That matrix makes the zi�<sup>1</sup> parallel to the next, zi-axis. The next matrix,

refers to the translation ai�1. This translation is along the xi�1-axis, defining the distance ai�1. As zi�<sup>1</sup> is parallel to zi, the distance ai�<sup>1</sup> is a mutual perpendicular to

100 ai�<sup>1</sup> 010 0 001 0 000 1

P

> cθ<sup>i</sup> �sθ<sup>i</sup> 0 0 sθ<sup>i</sup> cθ<sup>i</sup> 0 0 0 01 di 0 0 01

<sup>i</sup> T and the subscript to superscript means transform

� � for <sup>i</sup> <sup>¼</sup> <sup>1</sup>, <sup>2</sup>, <sup>3</sup> (7)

5, (8)

, (9)

si ¼ sθ<sup>i</sup> ¼ sin θ<sup>i</sup> (3) ci ¼ cθ<sup>i</sup> ¼ cos θ<sup>i</sup> (4)

c12 ¼ cð Þ¼ θ<sup>1</sup> þ θ<sup>2</sup> c1c2 � s1s2 ¼ cθ1cθ<sup>2</sup> � sθ1sθ<sup>2</sup> (5) s12 ¼ sð Þ¼ θ<sup>1</sup> þ θ<sup>2</sup> c1s2 þ s1c2 ¼ cθ1sθ<sup>2</sup> � sθ1cθ<sup>2</sup> (6)

3 i θi

cθ<sup>i</sup> �sθ<sup>i</sup> 0 0 sθ<sup>i</sup> cθ<sup>i</sup> 0 0 0 0 10 0 0 01

(2)

These abbreviations will be used throughout this chapter.

To include the translation in the same matrix notation, we use the homogeneous transform as shown in Eq. (1).

Before going to the other frames definition for the spacecraft, let us implement the matrix algebra, considering Figure 2. To clarify the procedure to obtain the general expression for the rotational-translational kinematics, let us redefine some frames in Figure 2.


Consider now Figure 2 that shows the various frames in addition to the D-H parameters.

To write the kinematics' equations in the frame f g i � 1 , we need to compute the rotation matrices using the homogeneous transform as below

Figure 3. Space robot on-orbit configuration and frames.

Kinematics for Spacecraft-Type Robotic Manipulators DOI: http://dx.doi.org/10.5772/intechopen.85636

connects to the joint 1. It is common to define joint 0 coincident with joint 1. In this case, a<sup>0</sup> is zero. For ground robots, it is common to define Ox0y0z<sup>0</sup> as an inertial frame. However, such definition is not consistent with the concept of inertial frames, mainly when we have space robot orbiting the Earth. For a space robotic manipulator, the Ox0y0z<sup>0</sup> frame is defined with its origin somewhere on the platform and locates the placement of joint zero. Figures 2 and 3 illustrate the frames, joint axes, joint variables, and the other fixed quantities according to the D-H convention.

To include the translation in the same matrix notation, we use the homogeneous

Before going to the other frames definition for the spacecraft, let us implement the matrix algebra, considering Figure 2. To clarify the procedure to obtain the general expression for the rotational-translational kinematics, let us redefine some

<sup>i</sup>�1, obtained be a rotation α<sup>i</sup>�<sup>1</sup> about xi�1, be

<sup>i</sup>�<sup>1</sup> translated by ai�<sup>1</sup> as f g Q ;

<sup>i</sup>�<sup>1</sup>, y<sup>0</sup> i, z<sup>0</sup>

• Consider frame f gi obtained by a translation di.

rotation matrices using the homogeneous transform as below

<sup>i</sup>�<sup>1</sup>, y<sup>0</sup>

<sup>i</sup>�1, z<sup>0</sup>

• Consider the frame obtained by a rotation θ<sup>i</sup> and parallel to frame f gi as f gP ;

Consider now Figure 2 that shows the various frames in addition to the D-H

To write the kinematics' equations in the frame f g i � 1 , we need to compute the

transform as shown in Eq. (1).

Kinematics - Analysis and Applications

• Let the frame formed by x<sup>0</sup>

• Consider frame formed by x<sup>0</sup>

frames in Figure 2.

parameters.

Figure 3.

108

Space robot on-orbit configuration and frames.

defined as f gR ;

i�1 <sup>P</sup> <sup>¼</sup> <sup>i</sup>�<sup>1</sup>    <sup>R</sup> <sup>T</sup><sup>R</sup> <sup>Q</sup> <sup>T</sup><sup>Q</sup> <sup>P</sup> T<sup>P</sup> <sup>i</sup> T    <sup>i</sup> <sup>P</sup> <sup>¼</sup> <sup>i</sup>�<sup>1</sup> <sup>i</sup> <sup>T</sup>    <sup>i</sup> P i�1      i T ¼ 10 0 0 0 cα<sup>i</sup>�<sup>1</sup> �sα<sup>i</sup>�<sup>1</sup> 0 0 sα<sup>i</sup>�<sup>1</sup> cα<sup>i</sup>�<sup>1</sup> 0 00 0 1 2 6 6 6 6 6 4 3 7 7 7 7 7 5 100ai�<sup>1</sup> 010 0 001 0 000 1 2 6 6 6 6 6 4 3 7 7 7 7 7 5 cθ<sup>i</sup> �sθ<sup>i</sup> 0 0 sθ<sup>i</sup> cθ<sup>i</sup> 0 0 0 0 10 0 0 01 2 6 6 6 6 6 4 3 7 7 7 7 7 5 1000 0100 001di 000 1 2 6 6 6 6 6 4 3 7 7 7 7 7 5 i�1     <sup>i</sup> <sup>T</sup> <sup>¼</sup> 10 0 ai�<sup>1</sup> 0 cα<sup>i</sup>�<sup>1</sup> �sα<sup>i</sup>�<sup>1</sup> 0 0 sα<sup>i</sup>�<sup>1</sup> cα<sup>i</sup>�<sup>1</sup> 0 00 0 1 2 6 6 6 6 6 4 3 7 7 7 7 7 5 cθ<sup>i</sup> �sθ<sup>i</sup> 0 0 sθ<sup>i</sup> cθ<sup>i</sup> 0 0 0 01 di 0 0 01 2 6 6 6 6 6 4 3 7 7 7 7 7 5 ¼ cθ<sup>i</sup> �sθ<sup>i</sup> 0 ai�<sup>1</sup> sθicα<sup>i</sup>�<sup>1</sup> cθicα<sup>i</sup>�<sup>1</sup> �sα<sup>i</sup>�<sup>1</sup> �sα<sup>i</sup>�<sup>1</sup>di sθisα<sup>i</sup>�<sup>1</sup> cθisα<sup>i</sup>�<sup>1</sup> cα<sup>i</sup>�<sup>1</sup> cα<sup>i</sup>�<sup>1</sup>di 0 00 1 2 6 6 6 6 6 4 3 7 7 7 7 7 5 (2)

where <sup>i</sup>�<sup>1</sup>     <sup>i</sup> <sup>T</sup> <sup>¼</sup> <sup>i</sup>�<sup>1</sup>     RT<sup>R</sup> <sup>Q</sup> <sup>T</sup><sup>Q</sup> <sup>P</sup> T<sup>P</sup> <sup>i</sup> T and the subscript to superscript means transform from the previous to the next coordinate frame.

These abbreviations will be used throughout this chapter.

$$\mathbf{s}\_{\mathbf{i}} = \mathbf{s}\theta\_{\mathbf{i}} = \sin\theta\_{\mathbf{i}} \tag{3}$$

$$\mathbf{c}\_{\mathbf{i}} = \mathbf{c}\theta\_{\mathbf{i}} = \cos\theta\_{\mathbf{i}}\tag{4}$$

$$\mathbf{c}\_{12} = \mathbf{c}(\theta\_1 + \theta\_2) = \mathbf{c}\_1 \mathbf{c}\_2 - \mathbf{s}\_1 \mathbf{s}\_2 = \mathbf{c}\theta\_1 \mathbf{c}\theta\_2 - \mathbf{s}\theta\_1 \mathbf{s}\theta\_2 \tag{5}$$

$$\mathbf{s}\_{12} = \mathbf{s}(\theta\_1 + \theta\_2) = \mathbf{c}\_1 \mathbf{s}\_2 + \mathbf{s}\_1 \mathbf{c}\_2 = \mathbf{c}\theta\_1 \mathbf{s}\theta\_2 - \mathbf{s}\theta\_1 \mathbf{c}\theta\_2 \tag{6}$$

$$\mathfrak{c}\_{123} = \mathfrak{c}\left(\sum\_{i}^{3} \theta\_{i}\right) \text{ and } \mathfrak{s}\_{123} = \mathfrak{s}\left(\sum\_{i}^{3} \theta\_{i}\right) \text{ for } i = 1, 2, 3 \tag{7}$$

In Eq. (2), the first matrix,

$$
\begin{bmatrix}
\mathbf{1} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\
\mathbf{0} & ca\_{i-1} & -sa\_{i-1} & \mathbf{0} \\
\mathbf{0} & sa\_{i-1} & ca\_{i-1} & \mathbf{0} \\
\mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{1}
\end{bmatrix},
\tag{8}
$$

is given by a rotation of α<sup>i</sup>�<sup>1</sup> about xi�1. Note that no translation appears in the matrix. That matrix makes the zi�<sup>1</sup> parallel to the next, zi-axis. The next matrix,

$$
\begin{bmatrix}
\mathbf{1} & \mathbf{0} & \mathbf{0} & a\_{i-1} \\
\mathbf{0} & \mathbf{1} & \mathbf{0} & \mathbf{0} \\
\mathbf{0} & \mathbf{0} & \mathbf{1} & \mathbf{0} \\
\mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{1}
\end{bmatrix},\tag{9}
$$

refers to the translation ai�1. This translation is along the xi�1-axis, defining the distance ai�1. As zi�<sup>1</sup> is parallel to zi, the distance ai�<sup>1</sup> is a mutual perpendicular to both the previous and the next z-axis. The rotation matrix,

$$
\begin{bmatrix} c\theta\_i & -s\theta\_i & \mathbf{0} & \mathbf{0} \\ s\theta\_i & c\theta\_i & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{1} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{1} \end{bmatrix},\tag{10}
$$

orbital velocity only. The practice is to use the Euler angles to describe the angular orientation between both frames (the attitude of the spacecraft with respect the

zl frame). According to the Euler Theorem, we can do this through three sequences of rotations from one frame to another frame. Figure 3 shows both frames in a nominal configuration in which the attitude angles are zero and both frames are coincident. As the spacecraft moves, the fixed frames also move. Figure 4 illustrates the rotation between those frames. So, we can derive the

Let us consider that the spacecraft fixed frame is misaligned from the Cxlyl

and Cxpypzp are misaligned. Figure 4b illustrates the sequence of rotations from

The sequence of three rotations is φ⟵θ⟵ψ [7]. The arrows indicate the direction of the sequence as: first, a rotation about zl by the angle ψ generating an

generating the auxiliary frame Cx}py}pz}p. Finally, the third rotation is done about x} to coincide with Cxpypzp. There are 12 sets of Euler angle sequences of rotations

φ⟵θ⟵ψ is one the most used transformations for satellites' rotational kinematics

Cxpypzp frames. The notation of the sequence used here is the same as φ�!θ�!ψ. However, the first notation emphasizes the resulting structure of matrix multipli-

> cψ sψ 0 �sψ cψ 0 0 01

8 >><

>>:

xl yl zl

9 >>=

>>;

to describe the rotational kinematics. However, the sequence of rotations

description. This transform provides the relationship between both Cxlyl

cation when implementing the transformation from Cxlyl

x0 p y0 p z0 p 9 >>=

>>; ¼

8 >><

>>:

(a) Satellite on-orbit configuration with Cxpypzp misaligned from Cxlyl

<sup>p</sup>. Then, a second rotation is done about y<sup>0</sup>

the angles as shown in Figure 4. Figure 4a illustrates the space robot in-orbit configuration with attitude different from zero. Note that in this figure, both Cxlyl

zl by

<sup>l</sup> by the angle θ

zl and

(12)

zl to Cxpypzp frames as we

zl and (b) the Euler sequence of

zl

kinematic expressions by using the spacecraft fixed frame, xpypzp.

Cxlyl

Cxlyl

zl to Cxpypzp.

are going to show next. First rotation, about zl

Figure 4.

111

rotations φ�!θ�!ψ:

py0 pz0

Kinematics for Spacecraft-Type Robotic Manipulators DOI: http://dx.doi.org/10.5772/intechopen.85636

auxiliary frame Cx<sup>0</sup>

results of rotation about zi by the joint angle θ<sup>i</sup> . This rotation generates a frame parallel to f gi at a distance di (offset) from the joint i. The 4th matrix,

$$
\begin{bmatrix}
\mathbf{1} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\
\mathbf{0} & \mathbf{1} & \mathbf{0} & \mathbf{0} \\
\mathbf{0} & \mathbf{0} & \mathbf{1} & d\_i \\
\mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{1}
\end{bmatrix},\tag{11}
$$

translates the frame by di completing the transform from frame f gi to f g i � 1 . Sometimes, it is more convenient to write the kinematic differential equations in the frame f gi . In this case, we have to transpose the rotation matrices so as to have the equations written in frame f gi but in coordinates of f g i � 1 .

For space robots, we must compute the kinematics not only of the manipulator but also for the whole spacecraft (rotational and translational kinematics) [3, 4]. To do this, it is necessary to define three more frames:


Figure 3 illustrates these new frames. The Cxlyl zl is fixed in the orbit. So, it follows with the orbit synchronized with the orbital velocity and with zl pointing toward the local vertical (Nadir direction). In the literature, this frame is also named RPY (roll, pitch, and yaw). Roll is the rotation angle about the orbital velocity direction (x-direction), pitch is the rotation angle about the negative orbit normal (y-axis), and yaw is the rotation angle about the local vertical (z-axis). The roll, pitch, and yaw angles yield the space robot attitude with respect to the Cxlyl zl frame. Both the Cxlyl zl and the Cxpypzp frames coincide with each other in a nominal attitude specification, that is, the space robot is pointing toward the local vertical, and the spacecraft attitude is zero. As the Cxpypzp frame is fixed in the spacecraft, the attitude is given by the angles between Cxlyl zl and Cxpypzp. Other conventions may be used for the definition of these frames. The reference systems ⨁XYZ, Cxlyl zl, and Cxpypzp are included in the formulation to allow for describing the rotational and translational kinematics in frame Cxpypzp or in the frame Cxlyl zl. In general, the equations of dynamics are written in the body fixed frame, Cxpypzp. We can write the kinematics and/or the kinematics' differential equations in any of those three frames by using transformation matrices. The {ECI} allows writing the translational kinematic expressions in terms of the orbit parameters [7]. If the coupling between the orbit and attitude motion is negligible (in general, it is), we can write the rotational kinematics in the Cxpypzp or Cxlyl zl frames, including the

#### Kinematics for Spacecraft-Type Robotic Manipulators DOI: http://dx.doi.org/10.5772/intechopen.85636

cθ<sup>i</sup> �sθ<sup>i</sup> 0 0 sθ<sup>i</sup> cθ<sup>i</sup> 0 0 0 0 10 0 0 01

results of rotation about zi by the joint angle θ<sup>i</sup> . This rotation generates a frame

1000 0100 001 di 000 1

translates the frame by di completing the transform from frame f gi to f g i � 1 . Sometimes, it is more convenient to write the kinematic differential equations in the frame f gi . In this case, we have to transpose the rotation matrices so as to have

For space robots, we must compute the kinematics not only of the manipulator but also for the whole spacecraft (rotational and translational kinematics) [3, 4]. To

• the Earth-centered inertial (ECI) frame ⨁XYZ, defined at the center of the

• the local vertical, local horizontal frame (LVLH) denoted as Cxlyl

center of mass of the space robot. The subscript l is a short for LVLH;

follows with the orbit synchronized with the orbital velocity and with zl pointing toward the local vertical (Nadir direction). In the literature, this frame is also named RPY (roll, pitch, and yaw). Roll is the rotation angle about the orbital velocity direction (x-direction), pitch is the rotation angle about the negative orbit normal (y-axis), and yaw is the rotation angle about the local vertical (z-axis). The roll, pitch, and yaw angles yield the space robot attitude with respect to the Cxlyl

nominal attitude specification, that is, the space robot is pointing toward the local vertical, and the spacecraft attitude is zero. As the Cxpypzp frame is fixed in the

conventions may be used for the definition of these frames. The reference systems

the rotational and translational kinematics in frame Cxpypzp or in the frame Cxlyl

In general, the equations of dynamics are written in the body fixed frame, Cxpypzp. We can write the kinematics and/or the kinematics' differential equations in any of those three frames by using transformation matrices. The {ECI} allows writing the translational kinematic expressions in terms of the orbit parameters [7]. If the coupling between the orbit and attitude motion is negligible (in general, it is), we

• and a spacecraft fixed system of reference, defined in the center of mass of the

zl and the Cxpypzp frames coincide with each other in a

zl, and Cxpypzp are included in the formulation to allow for describing

, (10)

5, (11)

zl is fixed in the orbit. So, it

zl and Cxpypzp. Other

zl frames, including the

zl. C is the

zl

zl.

Kinematics - Analysis and Applications

parallel to f gi at a distance di (offset) from the joint i. The 4th matrix,

the equations written in frame f gi but in coordinates of f g i � 1 .

do this, it is necessary to define three more frames:

space robots [7] and denoted by Cxpypzp.

frame. Both the Cxlyl

⨁XYZ, Cxlyl

110

Figure 3 illustrates these new frames. The Cxlyl

spacecraft, the attitude is given by the angles between Cxlyl

can write the rotational kinematics in the Cxpypzp or Cxlyl

Earth (the symbol ⨁ means Earth)

orbital velocity only. The practice is to use the Euler angles to describe the angular orientation between both frames (the attitude of the spacecraft with respect the Cxlyl zl frame). According to the Euler Theorem, we can do this through three sequences of rotations from one frame to another frame. Figure 3 shows both frames in a nominal configuration in which the attitude angles are zero and both frames are coincident. As the spacecraft moves, the fixed frames also move. Figure 4 illustrates the rotation between those frames. So, we can derive the kinematic expressions by using the spacecraft fixed frame, xpypzp.

Let us consider that the spacecraft fixed frame is misaligned from the Cxlyl zl by the angles as shown in Figure 4. Figure 4a illustrates the space robot in-orbit configuration with attitude different from zero. Note that in this figure, both Cxlyl zl and Cxpypzp are misaligned. Figure 4b illustrates the sequence of rotations from Cxlyl zl to Cxpypzp.

The sequence of three rotations is φ⟵θ⟵ψ [7]. The arrows indicate the direction of the sequence as: first, a rotation about zl by the angle ψ generating an auxiliary frame Cx<sup>0</sup> py0 pz0 <sup>p</sup>. Then, a second rotation is done about y<sup>0</sup> <sup>l</sup> by the angle θ generating the auxiliary frame Cx}py}pz}p. Finally, the third rotation is done about x} to coincide with Cxpypzp. There are 12 sets of Euler angle sequences of rotations to describe the rotational kinematics. However, the sequence of rotations φ⟵θ⟵ψ is one the most used transformations for satellites' rotational kinematics description. This transform provides the relationship between both Cxlyl zl and Cxpypzp frames. The notation of the sequence used here is the same as φ�!θ�!ψ. However, the first notation emphasizes the resulting structure of matrix multiplication when implementing the transformation from Cxlyl zl to Cxpypzp frames as we are going to show next.

First rotation, about zl

$$\begin{Bmatrix} \mathbf{x}'\_p \\ \mathbf{y}'\_p \\ \mathbf{z}'\_p \end{Bmatrix} = \begin{bmatrix} c\boldsymbol{\eta} & s\boldsymbol{\eta} & \mathbf{0} \\ -s\boldsymbol{\eta} & c\boldsymbol{\eta} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{1} \end{bmatrix} \begin{Bmatrix} \boldsymbol{x}\_l \\ \boldsymbol{\eta}\_l \\ \mathbf{z}\_l \end{Bmatrix} \tag{12}$$

#### Figure 4.

(a) Satellite on-orbit configuration with Cxpypzp misaligned from Cxlyl zl and (b) the Euler sequence of rotations φ�!θ�!ψ:

Second rotation, about y<sup>0</sup> <sup>p</sup> � y<sup>00</sup> p

$$
\begin{Bmatrix} \mathbf{x}^{\prime\prime}\_{\ p} \\ \mathbf{y}^{\prime\prime}\_{\ p} \\ \mathbf{z}^{\prime\prime}\_{\ p} \end{Bmatrix} = \begin{bmatrix} c\theta & \mathbf{0} & -s\theta \\ \mathbf{0} & \mathbf{1} & \mathbf{0} \\ s\theta & \mathbf{0} & c\theta \end{bmatrix} \begin{Bmatrix} \mathbf{x}^{\prime}\_{p} \\ \mathbf{y}^{\prime}\_{p} \\ \mathbf{z}^{\prime}\_{p} \end{Bmatrix} \tag{13}$$

p <sup>0</sup>T ¼

i�1      i 2 4

�100 0 01 0 10

T is given by Eq. (2)

robot is stabilized with Cxlyl

problem.

Table 1.

113

D-H parameters for Example 1.

3

Kinematics for Spacecraft-Type Robotic Manipulators DOI: http://dx.doi.org/10.5772/intechopen.85636

story and may be subject of research on another occasion.

0 <sup>1</sup>T ¼

1 <sup>2</sup>T ¼

5—this transform relates OxOyOzO with Cxpypzp (see Figure 3)

i ¼ 1, 2, …, n

The vector f g¼ rc rcx rcy rcz � �<sup>T</sup> represents the translation from the center of mass (CM) to the robot base or the OxOyOzO frame. This is a constant vector in this formulation. However, when the robot manipulator is in orbit operation, the links move, changing the location of the system CM, in other words, the CM location becomes a variable of the problem. By the time the space shuttle was in operation, this type of problem appeared. When it grasped significant massive target, the system CM motion affected the efficiency of the robot arm to put the target where it was planned. Another problem experienced by the space shuttle manipulator was the interaction between the attitude control system and the robot arm structural flexibility. The attitude control was in action to keep the shuttle stabilized while the manipulator was commanded to grasp a target. The structural flexibility of the long arms was excited by the attitude actuators and entered in elastic vibration mode causing problems to accomplish a safety grasping of the target. But this is another

Example 1. Consider the robot manipulator shown in Figure 3 where the space

written in the OxOyOzO frame. First, we built the D-H as shown in Table 1 for the

By inspection in Figure 3, we can see that all the z-axes are parallel and the joints

cθ<sup>1</sup> �sθ<sup>1</sup> 0 1 sθ<sup>1</sup> cθ<sup>1</sup> 0 0 0 0 10 0 0 01

cθ<sup>2</sup> �sθ<sup>2</sup> 0 l<sup>1</sup> sθ<sup>2</sup> cθ<sup>2</sup> 0 0 0 0 10 0 0 01

are revolute. So, the parameters α and d are zeros. Following the D-H table:

zl and Cypzp frames. Find the kinematics relationship

(20)

(21)

Third rotation, about x}<sup>p</sup> � xp

$$
\begin{Bmatrix} \boldsymbol{x}\_p \\ \boldsymbol{y}\_p \\ \boldsymbol{z}\_p \end{Bmatrix} = \begin{bmatrix} \mathbf{1} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & c\boldsymbol{\rho} & s\boldsymbol{\rho} \\ \mathbf{0} & -s\boldsymbol{\rho} & c\boldsymbol{\rho} \end{bmatrix} \begin{Bmatrix} \boldsymbol{x}\_p^\prime \\ \boldsymbol{y}^\prime \end{Bmatrix} \tag{14}$$

By substituting Eq. (13) into Eq. (14) and substituting Eq. (12) into the result, we have:

$$
\begin{Bmatrix} \mathbf{x}\_p \\ \mathbf{y}\_p \\ \mathbf{z}\_p \end{Bmatrix} = \begin{bmatrix} \mathbf{1} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & c\boldsymbol{\rho} & s\boldsymbol{\rho} \\ \mathbf{0} & -s\boldsymbol{\rho} & c\boldsymbol{\rho} \end{bmatrix} \begin{bmatrix} c\boldsymbol{\theta} & \mathbf{0} & -s\boldsymbol{\theta} \\ \mathbf{0} & \mathbf{1} & \mathbf{0} \\ s\boldsymbol{\theta} & \mathbf{0} & c\boldsymbol{\theta} \end{bmatrix} \begin{bmatrix} c\boldsymbol{\nu} & s\boldsymbol{\nu} & \mathbf{0} \\ -s\boldsymbol{\nu} & c\boldsymbol{\nu} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{1} \end{bmatrix} \begin{Bmatrix} \mathbf{x}\_l \\ \mathbf{y}\_l \\ \mathbf{z}\_l \end{Bmatrix} \tag{15}
$$

So, the structure is really like φ⟵θ⟵ψ.

The final matrix allows writing the equations in the spacecraft fixed frame Cxpypz<sup>p</sup> in coordinates of Cxlyl zl as:

$$\begin{Bmatrix} \mathbf{x}\_p \\ \mathbf{y}\_p \\ \mathbf{z}\_p \end{Bmatrix} = \begin{bmatrix} c\theta c\psi & c\theta s\psi & -s\theta \\ s\rho s\theta c\psi - c\rho s\psi & s\rho s\theta s\psi + c\rho c\psi & s\rho c\theta \\ c\rho s\theta c\psi + s\rho s\psi & c\rho s\theta s\psi - s\rho c\psi & c\rho c\theta \end{bmatrix} \begin{Bmatrix} \mathbf{x}\_l \\ \mathbf{y}\_l \\ \mathbf{z}\_l \end{Bmatrix} = \prescript{P}{}{\mathbf{R}} \tag{16}$$

By transposing this matrix, we can write the kinematics' equation in the LVLH frame:

$$\begin{Bmatrix} \mathbf{x}\_{\parallel} \\ \mathbf{y}\_{1} \\ \mathbf{z}\_{\parallel} \end{Bmatrix} = \begin{bmatrix} \mathbf{c}\theta\mathbf{c}\boldsymbol{\eta} & \mathbf{c}\theta\mathbf{s}\boldsymbol{\eta} & -\mathbf{c}\theta \\\\ \mathbf{s}\mathbf{q}\mathbf{s}\theta\mathbf{c}\boldsymbol{\eta} - \mathbf{c}\mathbf{q}\mathbf{s}\boldsymbol{\eta} & \mathbf{s}\mathbf{q}\mathbf{s}\theta\mathbf{s}\boldsymbol{\eta} + \mathbf{c}\mathbf{q}\mathbf{c}\boldsymbol{\eta} & \mathbf{s}\mathbf{q}\mathbf{c}\boldsymbol{\theta} \\\\ \mathbf{c}\mathbf{q}\mathbf{s}\theta\mathbf{c}\boldsymbol{\eta} + \mathbf{s}\mathbf{q}\mathbf{s}\boldsymbol{\eta} & \mathbf{c}\mathbf{q}\mathbf{s}\theta\mathbf{s}\boldsymbol{\eta} - \mathbf{s}\mathbf{q}\mathbf{c}\boldsymbol{\eta} & \mathbf{c}\mathbf{q}\mathbf{c}\boldsymbol{\theta} \end{bmatrix}^{\mathrm{T}} \begin{Bmatrix} \mathbf{x}\_{\mathbf{p}} \\ \mathbf{y}\_{\mathbf{p}} \\ \mathbf{z}\_{\mathbf{p}} \end{Bmatrix} = {}\_{\mathrm{p}}\mathbf{R}\{\mathbf{p}\} \tag{17}$$

The relationships from nf g to 0f g can be written as:

$$T\_n^0 T = \prescript{0}{1}{T\_2^1} T\_3^2 T \dots \prescript{n-1}{n}{T} \tag{18}$$

To write the kinematics of the spacecraft-like robot manipulator in the LVLH system of reference, we have to concatenate all the matrices. By using Eqs. (17) and (18) and Eq. (2), we can establish the relationships between i the coordinate frames as:

$$\mathbf{T}\_i^l T = \prescript{l}{}{T}\_0^p T^{i-1} \prescript{}{i}{T} \tag{19}$$

where

$$\begin{aligned} \prescript{I}{}{\prescript{I}{}}{T} = \begin{bmatrix} \prescript{I}{}{R}\mathbf{3}\times\mathbf{1} & \mathbf{0}\_{3\times 1} \\ \mathbf{0}\_{1\times 3} & \mathbf{1} \end{bmatrix} = \begin{bmatrix} c\theta c\psi & s\rho s\theta c\psi - c\rho s\psi & c\rho s\theta c\psi + s\rho s\psi & r\_{cx} \\ c\theta s\psi & s\rho s\theta s\psi + c\rho c\psi & s\rho c\theta & r\_{cy} \\ -s\theta & s\rho c\theta & c\rho c\theta & r\_{cx} \\ 0 & 0 & 0 & 1 \end{bmatrix}; \end{aligned}$$

Kinematics for Spacecraft-Type Robotic Manipulators DOI: http://dx.doi.org/10.5772/intechopen.85636

Second rotation, about y<sup>0</sup>

Kinematics - Analysis and Applications

Third rotation, about x}<sup>p</sup> � xp

have:

xp yp zp

> xp yp zp

9 >= >; ¼

8 ><

>:

xl yl zl

where

l <sup>p</sup>T ¼

112

l

<sup>p</sup>R<sup>3</sup>�<sup>1</sup> 03�<sup>1</sup> 01�<sup>3</sup> 1

" #

9 >= >; ¼ 2 6 4

frame:

8 ><

>:

9 >= >; ¼ 2 6 4

Cxpypz<sup>p</sup> in coordinates of Cxlyl

2 6 4

8 ><

>:

<sup>p</sup> � y<sup>00</sup> p

> 9 >= >; ¼

2 6 4

2 6 4 cθ 0 �sθ 01 0 sθ 0 cθ

10 0 0 cφ sφ 0 �sφ cφ

By substituting Eq. (13) into Eq. (14) and substituting Eq. (12) into the result, we

3 7 5

2 6 4

cθ 0 �sθ 01 0 sθ 0 cθ

The final matrix allows writing the equations in the spacecraft fixed frame

cθcψ cθsψ �sθ sφsθcψ � cφsψ sφsθsψ þ cφcψ sφcθ cφsθcψ þ sφsψ cφsθsψ � sφcψ cφcθ

cθcψ cθsψ �sθ sφsθcψ � cφsψ sφsθsψ þ cφcψ sφcθ cφsθcψ þ sφsψ cφsθsψ � sφcψ cφcθ

The relationships from nf g to 0f g can be written as:

¼

0 nT <sup>¼</sup> <sup>0</sup> 1T<sup>1</sup> 2T<sup>2</sup>

> l i <sup>T</sup> <sup>¼</sup> <sup>l</sup> PT<sup>p</sup> <sup>0</sup>T<sup>i</sup>�<sup>1</sup>      <sup>i</sup>

By transposing this matrix, we can write the kinematics' equation in the LVLH

To write the kinematics of the spacecraft-like robot manipulator in the LVLH system of reference, we have to concatenate all the matrices. By using Eqs. (17) and (18) and Eq. (2), we can establish the relationships between i the coordinate frames as:

3 7 5

3 7 5

8 ><

>:

8 ><

>:

x0 p y0 p z0 p

x}<sup>p</sup> y}<sup>p</sup> z}<sup>p</sup>

9 >=

>;

9 >=

>;

cψ sψ 0 �sψ cψ 0 0 01

> 3 7 5

3 7 5

cθcψ sφsθcψ � cφsψ cφsθcψ þ sφsψ rcx cθsψ sφsθsψ þ cφcψ sφcθ rcy �sθ sφcθ cφcθ rcz 0 0 01

8 ><

>:

<sup>T</sup> xp yp zp

8 ><

>:

xl yl zl

9 >= >; ¼ p

9 >= >; ¼ l

<sup>3</sup>T…<sup>n</sup>�<sup>1</sup>      nT (18)

T (19)

3 7 5

8 ><

>:

xl yl zl 9 >=

>;

<sup>l</sup> R (16)

<sup>p</sup>R pf g (17)

(13)

(14)

(15)

x00 p y00 p z00 p

xp yp zp 9 >= >; ¼

> 3 7 5

2 6 4

zl as:

8 ><

>:

8 ><

>:

10 0 0 cφ sφ 0 �sφ cφ

So, the structure is really like φ⟵θ⟵ψ.

p <sup>0</sup>T ¼ �100 0 01 0 10 2 4 3 5—this transform relates OxOyOzO with Cxpypzp (see Figure 3) i�1      i T is given by Eq. (2)

$$\mathbf{i} = \mathbf{1}, \mathbf{2}, \dots, \mathbf{n}$$

The vector f g¼ rc rcx rcy rcz � �<sup>T</sup> represents the translation from the center of mass (CM) to the robot base or the OxOyOzO frame. This is a constant vector in this formulation. However, when the robot manipulator is in orbit operation, the links move, changing the location of the system CM, in other words, the CM location becomes a variable of the problem. By the time the space shuttle was in operation, this type of problem appeared. When it grasped significant massive target, the system CM motion affected the efficiency of the robot arm to put the target where it was planned. Another problem experienced by the space shuttle manipulator was the interaction between the attitude control system and the robot arm structural flexibility. The attitude control was in action to keep the shuttle stabilized while the manipulator was commanded to grasp a target. The structural flexibility of the long arms was excited by the attitude actuators and entered in elastic vibration mode causing problems to accomplish a safety grasping of the target. But this is another story and may be subject of research on another occasion.

Example 1. Consider the robot manipulator shown in Figure 3 where the space robot is stabilized with Cxlyl zl and Cypzp frames. Find the kinematics relationship written in the OxOyOzO frame. First, we built the D-H as shown in Table 1 for the problem.

By inspection in Figure 3, we can see that all the z-axes are parallel and the joints are revolute. So, the parameters α and d are zeros. Following the D-H table:

$$\begin{aligned} \,^0\_1T &= \begin{bmatrix} c\theta\_1 & -s\theta\_1 & 0 & 1\\ s\theta\_1 & c\theta\_1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{bmatrix} \\ \,^1\_2T &= \begin{bmatrix} c\theta\_2 & -s\theta\_2 & 0 & l\_1\\ s\theta\_2 & c\theta\_2 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{bmatrix} \end{aligned} \tag{21}$$


Table 1. D-H parameters for Example 1.

$$\begin{aligned} \;^2\_3T = \begin{bmatrix} c\theta\_3 & -s\theta\_3 & 0 & l\_2 \\ s\theta\_3 & c\theta\_3 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \end{aligned} \tag{22}$$

new z-axis for the joint variable d2. The sequence of rotations is shown in

 

> 

  cθ<sup>1</sup> �sθ<sup>1</sup> 0 sθ<sup>1</sup> cθ<sup>1</sup> 0 0 0 a1

10 0 cα<sup>1</sup> �sα<sup>1</sup> cα<sup>1</sup> cα<sup>1</sup>

cθ<sup>3</sup> �sθ<sup>3</sup> 0 sθ<sup>3</sup> cθ<sup>3</sup> 0 0 01

To account for the translation d2, we can use the homogeneous transformation as:

cθ<sup>1</sup> �sθ<sup>1</sup> 0 0 sθ<sup>1</sup> cθ<sup>1</sup> 0 0 0 0 10 0 0 01

10 0 0 cα<sup>1</sup> �sα<sup>1</sup> 0 sα<sup>1</sup> cα<sup>1</sup> 0 00 0 1

cθ<sup>3</sup> �sθ<sup>3</sup> 0 0 sθ<sup>3</sup> cθ<sup>3</sup> 0 0 0 01 l<sup>2</sup> 0 0 01

Note that this notation of homogeneous matrix allows including the two translations of this problem. Also, note that the offset translation is along z.

Þ yields:

 

 

 

(27)

(28)

(29)

<sup>α</sup><sup>1</sup> <sup>¼</sup> <sup>90</sup><sup>o</sup> ð Þ (32)

(30)

(31)

(33)

 T ¼

(a) RPR manipulator showing the joints, (b) D-H parameters, and (c) frames.

Kinematics for Spacecraft-Type Robotic Manipulators DOI: http://dx.doi.org/10.5772/intechopen.85636

> T ¼

 T ¼

The transform matrix (α<sup>2</sup> <sup>¼</sup> <sup>90</sup><sup>o</sup>

100 0 d<sup>2</sup> 000 0

Eqs. (18)–(20).

Figure 5.

so

$$\begin{aligned} \;^0\_3T = \;^0\_1T^1\_2T^2\_3T = \begin{bmatrix} c\theta\_1 & -s\theta\_1 & 0 & 1\\ s\theta\_1 & c\theta\_1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} c\theta\_2 & -s\theta\_2 & 0 & l\_1\\ s\theta\_2 & c\theta\_2 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} c\theta\_3 & -s\theta\_3 & 0 & l\_2\\ s\theta\_3 & c\theta\_3 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{bmatrix} \tag{23} \end{aligned} \tag{23}$$

Compacting the result through the definition of c = cos, s = sin, and i ¼ θi, i ¼ 1, 2, 3 we have

$$\begin{aligned} \;\_3F\_3 = \begin{bmatrix} c\_3 c\_{12} - s\_3 s\_{12} & -c\_3 s\_{12} - s\_3 c\_{12} & 0 & l\_1 c\_1 + l\_2 c\_{12} \\ c\_3 s\_{12} + s\_3 c\_{12} & c\_3 c\_{12} - s\_3 s\_{12} & 0 & l\_1 s\_1 + l\_2 s\_{12} \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} = \begin{bmatrix} c\_{123} & -s\_{123} & 0 & l\_1 c\_1 + l\_2 c\_{12} \\ s\_{123} & c\_{123} & 0 & l\_1 s\_1 + l\_2 s\_{12} \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \end{aligned} \tag{24}$$

where

$$c\_{123} = \cos\left(\sum\_{i=1}^{3} \theta\_i\right) = c\_{3}c\_{12} - s\_{3}s\_{12} \tag{25}$$

$$s\_{123} = \sin\left(\sum\_{i=1}^{3} \theta\_i\right) = c\_{3}s\_{12} + s\_{3}c\_{12} \tag{26}$$
 
$$\mathbf{i} = \mathbf{1, 2, 3}$$

Note that we considered the spacecraft stable in attitude with the Cxpypzp aligned to the Cxlyl zl frame. But be aware that if the manipulator moves in attitude, the kinematics' differential equations of the spacecraft shall account for the manipulator links of kinematics' differential equations. In particular, if the manipulator is to grasp an on-orbit target, the attitude of the manipulator-like spacecraft must be synchronized to that of the target so that the relative attitude of both spacecraft and target is zero. As the dynamic analysis is out of the scope of this chapter, we will assume that the robot arm kinematics' equations differ from ground robots only by Cxlyl zl and xpypzp frames. For orbital robots these frames are used to describe the kinematics differential equations related to orbit and attitude motions.

Example 2. Consider Figure 5 which illustrates a revolute-prismatic-revolute (RPR) manipulator. In this case, we have a rotation ð Þ θ<sup>1</sup> , a translational variable ð Þ d<sup>2</sup> , and another rotation ð Þ θ<sup>2</sup> as illustrated in Table 2. The convention is that those degrees-of-freedom occur along z-axis.

For this case, we have one rotation of θ1, translation of d2, and a rotation of θ3. Following the D-H convention above, we should make a few comments. First, we have a rotation on z-axis. Then, we use the twist angle rotation to make a

Kinematics for Spacecraft-Type Robotic Manipulators DOI: http://dx.doi.org/10.5772/intechopen.85636

 T ¼

cθ<sup>1</sup> �sθ<sup>1</sup> 0 1 sθ<sup>1</sup> cθ<sup>1</sup> 0 0 0 0 10 0 0 01

c3c<sup>12</sup> � s3s<sup>12</sup> �c3s<sup>12</sup> � s3c<sup>12</sup> 0 l1c<sup>1</sup> þ l2c<sup>12</sup> c3s<sup>12</sup> þ s3c<sup>12</sup> c3c<sup>12</sup> � s3s<sup>12</sup> 0 l1s<sup>1</sup> þ l2s<sup>12</sup> 0 0 10 0 0 01

c<sup>123</sup> ¼ cos ∑

s<sup>123</sup> ¼ sin ∑

 i¼1 θi � �

 i¼1 θi � �

i ¼ 1, 2, 3

the kinematics' differential equations of the spacecraft shall account for the manipulator links of kinematics' differential equations. In particular, if the manipulator is to grasp an on-orbit target, the attitude of the manipulator-like spacecraft must be synchronized to that of the target so that the relative attitude of both spacecraft and target is zero. As the dynamic analysis is out of the scope of this chapter, we will assume that the robot arm kinematics' equations differ from ground robots only by

zl and xpypzp frames. For orbital robots these frames are used to describe the

Example 2. Consider Figure 5 which illustrates a revolute-prismatic-revolute (RPR) manipulator. In this case, we have a rotation ð Þ θ<sup>1</sup> , a translational variable ð Þ d<sup>2</sup> , and another rotation ð Þ θ<sup>2</sup> as illustrated in Table 2. The convention is that

For this case, we have one rotation of θ1, translation of d2, and a rotation of θ3. Following the D-H convention above, we should make a few comments. First,

we have a rotation on z-axis. Then, we use the twist angle rotation to make a

kinematics differential equations related to orbit and attitude motions.

those degrees-of-freedom occur along z-axis.

zl frame. But be aware that if the manipulator moves in attitude,

Note that we considered the spacecraft stable in attitude with the Cxpypzp

Kinematics - Analysis and Applications

i ¼ θi, i ¼ 1, 2, 3 we have

where

aligned to the Cxlyl

Cxlyl

so

 <sup>T</sup> <sup>¼</sup> <sup>0</sup> T<sup>1</sup> T<sup>2</sup> T ¼

 T ¼

Compacting the result through the definition of c = cos, s = sin, and

cθ<sup>3</sup> �sθ<sup>3</sup> 0 l<sup>2</sup> sθ<sup>3</sup> cθ<sup>3</sup> 0 0 0 0 10 0 0 01

cθ<sup>2</sup> �sθ<sup>2</sup> 0 l<sup>1</sup> sθ<sup>2</sup> cθ<sup>2</sup> 0 0 0 0 10 0 0 01

(22)

(23)

(24)

cθ<sup>3</sup> �sθ<sup>3</sup> 0 l<sup>2</sup> sθ<sup>3</sup> cθ<sup>3</sup> 0 0 0 0 10 0 0 01

c<sup>123</sup> �s<sup>123</sup> 0 l1c<sup>1</sup> þ l2c<sup>12</sup> s<sup>123</sup> c<sup>123</sup> 0 l1s<sup>1</sup> þ l2s<sup>12</sup> 0 01 0 0 00 1

c3c<sup>12</sup> � s3s<sup>12</sup> (25)

c3s<sup>12</sup> þ s3c<sup>12</sup> (26)

Figure 5. (a) RPR manipulator showing the joints, (b) D-H parameters, and (c) frames.

new z-axis for the joint variable d2. The sequence of rotations is shown in Eqs. (18)–(20).

$$\begin{aligned} \, \, \_1^0 T &= \begin{bmatrix} c\theta\_1 & -s\theta\_1 & 0\\ s\theta\_1 & c\theta\_1 & 0\\ 0 & 0 & a\mathbf{1} \end{bmatrix} \\ \, \_2^1 T &= \begin{bmatrix} \mathbf{1} & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & ca\_1 & -s\alpha\_1\\ \mathbf{0} & ca\_1 & ca\_1 \end{bmatrix} \\ \, \_3^2 T &= \begin{bmatrix} c\theta\_3 & -s\theta\_3 & \mathbf{0} \\ s\theta\_3 & c\theta\_3 & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{1} \end{bmatrix} \end{aligned} \tag{28}$$

To account for the translation d2, we can use the homogeneous transformation as:

$$\begin{bmatrix} c\theta\_1 & -s\theta\_1 & 0 & 0\\ s\theta\_1 & c\theta\_1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{bmatrix} \tag{30}$$

$$\begin{bmatrix} 1 & 0 & 0 & 0\\ 0 & c\alpha\_1 & -s\alpha\_1 & 0\\ 0 & s\alpha\_1 & c\alpha\_1 & 0\\ 0 & 0 & 0 & 1 \end{bmatrix} \tag{31}$$

$$\begin{bmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & d\_2\\ 0 & 0 & 0 & 0 \end{bmatrix} (a\_1 = 90^\circ) \tag{32}$$

$$\begin{bmatrix} c\theta\_3 & -s\theta\_3 & 0 & 0\\ s\theta\_3 & c\theta\_3 & 0 & 0\\ 0 & 0 & 1 & l\_2\\ 0 & 0 & 0 & 1 \end{bmatrix} \tag{33}$$

Note that this notation of homogeneous matrix allows including the two translations of this problem. Also, note that the offset translation is along z. The transform matrix (α<sup>2</sup> <sup>¼</sup> <sup>90</sup><sup>o</sup> Þ yields:

Kinematics - Analysis and Applications

$$
\begin{aligned}
\begin{bmatrix}1 & 0 & 0 & 0 \\ 0 & ca\_1 & -sa\_1 & 0 \\ 0 & sa\_1 & ca\_1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix}1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & d\_2 \\ 0 & 0 & 0 & 0 \end{bmatrix} &= \begin{bmatrix}1 & 0 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix}1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & d\_2 \\ 0 & 0 & 0 & 0 \end{bmatrix} \\ &= \begin{bmatrix}1 & 0 & 0 & 0 \\ 0 & 0 & -1 & -d\_2 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \end{aligned} \tag{34}
$$

the space robotic system. We know the location of that system of reference. Also, we

previous section, we considered the relationship between frames to write the equations in the {0}. In the inverse problem, we want to find joint variables given the specified position and orientation of the tool relative to {0} or {p}. We assume that the robot attitude and orbit control subsystem keeps the space robot stabilized while it works. This is to simplify the formulation. This assumption is reasonable since this is exactly what must be done when the space robot is in operation. We have already presented and discussed the several frames defined to study the robot kinematics. We

There are some aspects of inverse kinematics that make it more complex than forward kinematics [5, 6]. We will briefly report them, but we will not go into too many details. In this chapter, our focus is on the solution methods for inverse

One of the main aspects of the inverse kinematics is solvability because we have to face the problem of existence of solutions and of multiple solutions. Also, this problem brings to light the concept of manipulator workspace which is the volume of space that the end-effector of the manipulator can reach. For space robot, this prob-

In solving inverse problems, the first thing we must consider is that there are no general algorithms [5, 8, 9] that may be employed to solve manipulator kinematics. We state that a manipulator is solvable when the joint variables can be determined

The solutions may be classified as close solution and numerical solutions. Numerical solutions are in general slower than close solutions and would require another chapter for discussion. We will concentrate in the close solutions methods. These methods can be subdivided into algebraic and geometric solutions. This classification is somehow hazy since each solution method uses both algebraic and geometric kinematic formulations. In presenting and discussing the methods, we point out that according to the concept of solution we stated here, the revolute and prismatic manipulators having six degrees of freedom (DOF) in series chain are solvable.

Consider the robot manipulator mounted on the spacecraft platform and the onorbit configurations in the nominal attitude as shown in Figure 3. Let us apply the

by an algorithm for a given position and orientation of the tool frame, {T}.

lem is critical because if the object to be manipulated is beyond space robot workspace, the space mission may not accomplish goals like grasping a target, docking to another spacecraft, put a piece of structure in right location during assembling space structures, and so on. For a solution to exist, the point specified to be reached by the end-effector must be inside the manipulator workspace. Another problem related to the workspace of robots is to plan a trajectory [9] preventing collisions with objects, another robot, or astronaut inside the workspace. There are two useful definitions related to workspace. One is the dextrous workspace [5] and the other is the reachable workspace. The first refers to the space volume that the end-effector reaches with all orientation. The second is that volume of space that the robot can reach in at least one orientation. If the position and orientation specified for the robot task is inside the workspace, then at least one solution exists. Workspace also depends on the tool-frame transformation, since it is frequently the tool-tip that is discussed when we discuss reachable points in space. Details about the solvability

will refer to just one more, the frame defined on the tool, {T}.

problem for inverse kinematics may be seen in [5].

3.1 The inverse kinematics solution methods

3.2 Algebraic solution

zl frame because we know the space robot orbit. In the

know the location of the Cxlyl

Kinematics for Spacecraft-Type Robotic Manipulators DOI: http://dx.doi.org/10.5772/intechopen.85636

kinematics problem.

This result is a transformation, <sup>1</sup> T <sup>2</sup> , that is, the transform from frame {2} to frame {1}. Combining all the transformations, we can write from frame {3} to frame {0} as:

$$\begin{aligned} \;^0\_3T &= {}^0\_1T\_2 {}^2\_2T\_3^2T \\ &= \begin{bmatrix} c\theta\_1 & -s\theta\_2 & 0 & 0\\ s\theta\_2 & c\theta\_1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 & 0\\ 0 & 0 & -1 & -d\_2\\ 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} c\theta\_3 & -s\theta\_3 & 0 & 0\\ s\theta\_3 & c\theta\_3 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{bmatrix} \\ &= \begin{bmatrix} c\theta\_1 c\theta\_3 & -c\theta\_2 s\theta\_3 & s\theta\_3 & s\theta\_1 (l\_2 + d\_2)\\ s\theta\_1 c\theta\_3 & -s\theta\_2 s\theta\_3 & -c\theta\_1 & -c\theta\_1 (l\_2 + d\_2)\\ s\theta\_3 & c\theta\_3 & 0 & 0\\ 0 & 0 & 0 & 1 \end{bmatrix} \end{aligned} \tag{35}$$

We can also find the transpose of this transformation to have the transformation:

$$\begin{aligned} \;^3\_0T &= \,^3\_2T^2\_1T^1\_0T \\ \begin{bmatrix} \;^3\_0T & \;^3\_0\theta\_3 & \;0 & 0\\ -s\theta\_3 & c\theta\_3 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 & 0\\ 0 & 0 & 1 & -d\_2\\ 0 & -1 & 0 & 0\\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} c\theta\_1 & s\theta\_1 & 0 & 0\\ -s\theta\_1 & c\theta\_1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{bmatrix} \\ &= \begin{bmatrix} c\theta\_1c\theta\_3 & s\theta\_1c\theta\_3 & s\theta\_3 & 0\\ -c\theta\_1s\theta\_3 & -s\theta\_2s\theta\_3 & c\theta\_3 & 0\\ s\theta\_1 & -c\theta\_1 & 0 & -(l\_2+d\_2)\\ 0 & 0 & 0 & 1 \end{bmatrix} \end{aligned} \tag{36}$$

#### 3. Inverse kinematics

Inverse kinematics consists of finding the joint variables given the robot manipulator end-effector position and orientation with respect to the user workstation. The workstation means a known location and reference system. It could be the platform of the robotic system. In this chapter, the platform reference system is the frame Cxpypzp represented by {P}. In general, such frame is defined in the center of mass of

#### Kinematics for Spacecraft-Type Robotic Manipulators DOI: http://dx.doi.org/10.5772/intechopen.85636

10 0 0 cα<sup>1</sup> �sα<sup>1</sup> 0 sα<sup>1</sup> cα<sup>1</sup> 0 00 0 1

Kinematics - Analysis and Applications

This result is a transformation, <sup>1</sup>

T <sup>2</sup>

cθ<sup>1</sup> �sθ<sup>2</sup> 0 0 sθ<sup>2</sup> cθ<sup>1</sup> 0 0 0 0 10 0 0 01

cθ<sup>3</sup> sθ<sup>3</sup> 0 0 �sθ<sup>3</sup> cθ<sup>3</sup> 0 0 0 0 10 0 0 01

cθ1cθ<sup>3</sup> sθ1cθ<sup>3</sup> sθ<sup>3</sup> 0 �cθ1sθ<sup>3</sup> �sθ1sθ<sup>3</sup> cθ<sup>3</sup> 0

sθ<sup>1</sup> �cθ<sup>1</sup> 0 �ð Þ l<sup>2</sup> þ d<sup>2</sup> 0 00 1

{1}. Combining all the transformations, we can write from frame {3} to frame {0} as:

10 0 0 0 0 �1 �d<sup>2</sup> 01 0 0 00 0 1

cθ1cθ<sup>3</sup> �cθ1sθ<sup>3</sup> sθ<sup>3</sup> sθ1ð Þ l<sup>2</sup> þ d<sup>2</sup> sθ1cθ<sup>3</sup> �sθ1sθ<sup>3</sup> �cθ<sup>1</sup> �cθ1ð Þ l<sup>2</sup> þ d<sup>2</sup> sθ<sup>3</sup> cθ<sup>3</sup> 0 0 0 00 1

We can also find the transpose of this transformation to have the transformation:

Inverse kinematics consists of finding the joint variables given the robot manipulator end-effector position and orientation with respect to the user workstation. The workstation means a known location and reference system. It could be the platform of the robotic system. In this chapter, the platform reference system is the frame Cxpypzp represented by {P}. In general, such frame is defined in the center of mass of

100 0 �d<sup>2</sup> �10 0 000 1

10 0 0 0 0 �1 �d<sup>2</sup> 01 0 0 00 0 1

, that is, the transform from frame {2} to frame

cθ<sup>3</sup> �sθ<sup>3</sup> 0 0 sθ<sup>3</sup> cθ<sup>3</sup> 0 0 0 0 10 0 0 01

cθ<sup>1</sup> sθ<sup>1</sup> 0 0 �sθ<sup>1</sup> cθ<sup>1</sup> 0 0 0 0 10 0 0 01

(34)

(35)

(36)

 <sup>T</sup> <sup>¼</sup><sup>0</sup> T<sup>1</sup> T<sup>2</sup> T ¼

 <sup>T</sup> <sup>¼</sup> <sup>3</sup> T<sup>2</sup> T<sup>1</sup> T

3. Inverse kinematics

the space robotic system. We know the location of that system of reference. Also, we know the location of the Cxlyl zl frame because we know the space robot orbit. In the previous section, we considered the relationship between frames to write the equations in the {0}. In the inverse problem, we want to find joint variables given the specified position and orientation of the tool relative to {0} or {p}. We assume that the robot attitude and orbit control subsystem keeps the space robot stabilized while it works. This is to simplify the formulation. This assumption is reasonable since this is exactly what must be done when the space robot is in operation. We have already presented and discussed the several frames defined to study the robot kinematics. We will refer to just one more, the frame defined on the tool, {T}.

There are some aspects of inverse kinematics that make it more complex than forward kinematics [5, 6]. We will briefly report them, but we will not go into too many details. In this chapter, our focus is on the solution methods for inverse kinematics problem.

One of the main aspects of the inverse kinematics is solvability because we have to face the problem of existence of solutions and of multiple solutions. Also, this problem brings to light the concept of manipulator workspace which is the volume of space that the end-effector of the manipulator can reach. For space robot, this problem is critical because if the object to be manipulated is beyond space robot workspace, the space mission may not accomplish goals like grasping a target, docking to another spacecraft, put a piece of structure in right location during assembling space structures, and so on. For a solution to exist, the point specified to be reached by the end-effector must be inside the manipulator workspace. Another problem related to the workspace of robots is to plan a trajectory [9] preventing collisions with objects, another robot, or astronaut inside the workspace. There are two useful definitions related to workspace. One is the dextrous workspace [5] and the other is the reachable workspace. The first refers to the space volume that the end-effector reaches with all orientation. The second is that volume of space that the robot can reach in at least one orientation. If the position and orientation specified for the robot task is inside the workspace, then at least one solution exists. Workspace also depends on the tool-frame transformation, since it is frequently the tool-tip that is discussed when we discuss reachable points in space. Details about the solvability problem for inverse kinematics may be seen in [5].

#### 3.1 The inverse kinematics solution methods

In solving inverse problems, the first thing we must consider is that there are no general algorithms [5, 8, 9] that may be employed to solve manipulator kinematics. We state that a manipulator is solvable when the joint variables can be determined by an algorithm for a given position and orientation of the tool frame, {T}.

The solutions may be classified as close solution and numerical solutions. Numerical solutions are in general slower than close solutions and would require another chapter for discussion. We will concentrate in the close solutions methods. These methods can be subdivided into algebraic and geometric solutions. This classification is somehow hazy since each solution method uses both algebraic and geometric kinematic formulations. In presenting and discussing the methods, we point out that according to the concept of solution we stated here, the revolute and prismatic manipulators having six degrees of freedom (DOF) in series chain are solvable.

#### 3.2 Algebraic solution

Consider the robot manipulator mounted on the spacecraft platform and the onorbit configurations in the nominal attitude as shown in Figure 3. Let us apply the

algebraic solution. Note that the frame defined in the last link is the frame for n = 3. So, the transform from frame {0} to frame {N} is given by Eq. (24) repeated here for didactic reasons.

$$\begin{aligned} \, \_N^O T = \begin{bmatrix} c\_{123} & -s\_{123} & 0 & l\_1 c\_1 + l\_2 c\_{12} \\ s\_{123} & c\_{123} & 0 & l\_1 s\_1 + l\_2 s\_{12} \\ \mathbf{0} & \mathbf{0} & \mathbf{1} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{1} \end{bmatrix} \end{aligned} \tag{37}$$

Let us assume that the orientation and position specified in inside the

<sup>s</sup><sup>2</sup> ¼ � ffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 � c<sup>2</sup>

<sup>θ</sup><sup>2</sup> <sup>¼</sup> Atan2 <sup>s</sup><sup>2</sup>

Note that there are signs to be considered in Eq. (45). The signs are related to

Next, we make a change of variable by writing k<sup>1</sup> and k<sup>2</sup> in the form:

q

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k1 <sup>2</sup> <sup>þ</sup> <sup>k</sup><sup>2</sup> 2

r ¼ þ

Using these new variables, we can write Eqs. (51) and (52) as:

<sup>γ</sup> <sup>þ</sup> <sup>θ</sup><sup>1</sup> <sup>¼</sup> Atan2 <sup>y</sup>

Then, using the definition of γ, we obtain the expression for θ<sup>1</sup> as:

r ; x r

x

y

Once we have obtained θ2, we can obtain θ<sup>1</sup> by manipulating Eqs. (41) and (42).

c2

c2

<sup>2</sup> p (45)

� � (46)

� �, but uses the signs of both <sup>s</sup><sup>2</sup> and <sup>c</sup><sup>2</sup> to identify

k<sup>1</sup> ¼ l<sup>1</sup> þ l2c<sup>2</sup> (47) k<sup>2</sup> ¼ s<sup>2</sup> (48)

c<sup>12</sup> ¼ c1c<sup>2</sup> � s1s<sup>2</sup> (49) s<sup>12</sup> ¼ c1s<sup>2</sup> þ s1c<sup>2</sup> (50)

x ¼ k1c<sup>1</sup> � k2s<sup>1</sup> (51) x ¼ k1s<sup>1</sup> þ k2c<sup>1</sup> (52)

γ ¼ Atan2ð Þ k1; k<sup>2</sup> (54)

<sup>r</sup> <sup>¼</sup> <sup>c</sup>γcθ<sup>1</sup> � <sup>s</sup>γcsθ<sup>1</sup> <sup>¼</sup> <sup>c</sup>ð Þ <sup>γ</sup> <sup>þ</sup> <sup>θ</sup><sup>1</sup> (55)

<sup>r</sup> <sup>¼</sup> <sup>c</sup>γsθ<sup>1</sup> <sup>þ</sup> <sup>s</sup>γθ<sup>1</sup> <sup>¼</sup> <sup>s</sup>ð Þ <sup>γ</sup> <sup>þ</sup> <sup>θ</sup><sup>1</sup> (56)

θ<sup>1</sup> ¼ Atan2ð Þ� y; x Atan2ð Þ k1; k<sup>2</sup> (58)

� � <sup>¼</sup> Atan2ð Þ <sup>y</sup>; <sup>x</sup> (57)

(53)

workspace. In this case, we can compute s<sup>2</sup> as:

Kinematics for Spacecraft-Type Robotic Manipulators DOI: http://dx.doi.org/10.5772/intechopen.85636

The joint angle can then be computed by:

� � computes tan �<sup>1</sup> <sup>s</sup><sup>2</sup>

To solve the problem, we should define:

and use the trigonometric identities

to rewrite Eqs. (41) and (42) as:

the quadrant in which the resulting angle lies.

Atan2 <sup>s</sup><sup>2</sup> c2

multiple solutions.

and

from which

119

Now consider that for inverse kinematics, the specified position of the endeffector is given:

$$\begin{aligned} \, \_N^O T = \begin{bmatrix} c\rho & -s\rho & 0 & \mathbf{x} \\ s\rho & c\rho & \mathbf{0} & \mathbf{y} \\ \mathbf{0} & \mathbf{0} & \mathbf{1} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{1} \end{bmatrix} \end{aligned} \tag{38}$$

Note that this specified orientation and position are the angle φ and the Cartesian position (x, y, 0) of the tool. As Eq. (37) must be equal to Eq. (38), the problem can be solved by comparing the elements of the matrices. The joint angles and the parameters li, i ¼ 1, 2 are unknown. The specified orientation given by φ and the components of the vector f g xyz <sup>T</sup>, in the planar case, x and y are known and represent the tool position with respect to {0}. By equating a<sup>11</sup> element of both matrices we have:

$$c\rho = c\_{123} = \cos\left(\sum\_{i=1}^{3} \theta\_i\right) \tag{39}$$

$$s\rho = s\_{123} = \cos\left(\sum\_{i=1}^{3} \theta\_i\right) \tag{40}$$

$$\mathbf{x} = l\_1 \mathbf{c}\_1 + l\_2 \mathbf{c}\_{12} \tag{41}$$

$$y = l\_1 s\_1 + l\_2 s\_{12} \tag{42}$$

where

$$\mathbf{c}\_{12} = \mathbf{c}\_1 \mathbf{c}\_2 - \mathbf{s}\_1 \mathbf{s}\_2$$

$$\mathbf{s}\_{12} = \mathbf{c}\_1 \mathbf{s}\_2 + \mathbf{s}\_1 \mathbf{c}\_2$$

The algebraic solution requires the solution of Eq. (39) to Eq. (42). We can compute the square of Eqs. (41) and (42) to obtain:

$$\mathbf{x}^2 + \mathbf{y}^2 = l\_1^2 + l\_2^2 + 2l\_1l\_2c\_2 \tag{43}$$

Solving Eq. (43) by c2, we compute c<sup>2</sup> as:

$$c\_2 = \frac{x^2 + y^2 - \left(l\_1^2 + l\_2^2\right)}{2l\_1l\_2} \tag{44}$$

Important information can be extracted from Eq. (44). We know that the cosine function must be between �1 and + 1. By checking such constraint, we obtain information whether or not the solution exists. In case a solution does not exist, the meaning is that the point is not inside the robot manipulator workspace.

algebraic solution. Note that the frame defined in the last link is the frame for n = 3. So, the transform from frame {0} to frame {N} is given by Eq. (24) repeated here

Now consider that for inverse kinematics, the specified position of the end-

c<sup>123</sup> �s<sup>123</sup> 0 l1c<sup>1</sup> þ l2c<sup>12</sup> s<sup>123</sup> c<sup>123</sup> 0 l1s<sup>1</sup> þ l2s<sup>12</sup> 0 01 0 0 00 1

> cφ �sφ 0 x sφ cφ 0 y 0 0 10 0 0 01

Note that this specified orientation and position are the angle φ and the Cartesian position (x, y, 0) of the tool. As Eq. (37) must be equal to Eq. (38), the problem can be solved by comparing the elements of the matrices. The joint angles and the parameters li, i ¼ 1, 2 are unknown. The specified orientation given by φ and the

represent the tool position with respect to {0}. By equating a<sup>11</sup> element of both

cφ ¼ c<sup>123</sup> ¼ cos ∑

sφ ¼ s<sup>123</sup> ¼ cos ∑

c12 ¼ c1c2 � s1s2 s12 ¼ c1s2 þ s1c2

The algebraic solution requires the solution of Eq. (39) to Eq. (42). We can

<sup>2</sup> <sup>þ</sup> <sup>l</sup><sup>2</sup>

2l1l<sup>2</sup>

Important information can be extracted from Eq. (44). We know that the cosine

<sup>2</sup> <sup>þ</sup> <sup>l</sup><sup>2</sup> <sup>2</sup> � �

<sup>x</sup><sup>2</sup> <sup>þ</sup> <sup>y</sup><sup>2</sup> <sup>¼</sup> <sup>l</sup><sup>1</sup>

<sup>c</sup><sup>2</sup> <sup>¼</sup> <sup>x</sup><sup>2</sup> <sup>þ</sup> <sup>y</sup><sup>2</sup> � <sup>l</sup><sup>1</sup>

function must be between �1 and + 1. By checking such constraint, we obtain information whether or not the solution exists. In case a solution does not exist, the

meaning is that the point is not inside the robot manipulator workspace.

compute the square of Eqs. (41) and (42) to obtain:

Solving Eq. (43) by c2, we compute c<sup>2</sup> as:

3 i¼1 θi � �

3 i¼1 θi � �

<sup>T</sup>, in the planar case, x and y are known and

x ¼ l1c<sup>1</sup> þ l2c<sup>12</sup> (41) y ¼ l1s<sup>1</sup> þ l2s<sup>12</sup> (42)

<sup>2</sup> <sup>þ</sup> <sup>2</sup>l1l2c<sup>2</sup> (43)

(37)

(38)

(39)

(40)

(44)

for didactic reasons.

Kinematics - Analysis and Applications

effector is given:

matrices we have:

where

118

O NT ¼

components of the vector f g xyz

O <sup>N</sup>T ¼

Let us assume that the orientation and position specified in inside the workspace. In this case, we can compute s<sup>2</sup> as:

$$
\sigma\_2 = \pm \sqrt{1 - c\_2^2} \tag{45}
$$

The joint angle can then be computed by:

$$\theta\_2 = \text{Atom2}\left(\frac{\varsigma\_2}{\varsigma\_2}\right) \tag{46}$$

Atan2 <sup>s</sup><sup>2</sup> c2 � � computes tan �<sup>1</sup> <sup>s</sup><sup>2</sup> c2 � �, but uses the signs of both <sup>s</sup><sup>2</sup> and <sup>c</sup><sup>2</sup> to identify the quadrant in which the resulting angle lies.

Note that there are signs to be considered in Eq. (45). The signs are related to multiple solutions.

Once we have obtained θ2, we can obtain θ<sup>1</sup> by manipulating Eqs. (41) and (42). To solve the problem, we should define:

$$k\_1 = l\_1 + l\_2 c\_2 \tag{47}$$

$$k\_2 = s\_2 \tag{48}$$

and use the trigonometric identities

$$
\sigma\_{12} = \sigma\_1 \sigma\_2 - \mathfrak{s}\_1 \mathfrak{s}\_2 \tag{49}
$$

$$\mathfrak{s}\_{12} = \mathfrak{c}\_{1}\mathfrak{s}\_{2} + \mathfrak{s}\_{1}\mathfrak{c}\_{2} \tag{50}$$

to rewrite Eqs. (41) and (42) as:

$$\mathbf{x} = k\_1 \mathbf{c}\_1 - k\_2 \mathbf{s}\_1 \tag{51}$$

$$x = k\_1 s\_1 + k\_2 c\_1 \tag{52}$$

Next, we make a change of variable by writing k<sup>1</sup> and k<sup>2</sup> in the form:

$$r = +\sqrt{{k\_1}^2 + {k\_2}^2} \tag{53}$$

and

$$\gamma = \text{Atom2}(k\_1, k\_2) \tag{54}$$

Using these new variables, we can write Eqs. (51) and (52) as:

$$\frac{\infty}{r} = c\gamma c\theta\_1 - s\gamma c\varsigma\theta\_1 = c(\chi + \theta\_1) \tag{55}$$

$$\frac{\mathcal{Y}}{r} = c\gamma s\theta\_1 + s\gamma\theta\_1 = s(\gamma + \theta\_1) \tag{56}$$

from which

$$\gamma + \theta\_1 = \text{Atom2}\left(\frac{\mathcal{Y}}{r}, \frac{\infty}{r}\right) = \text{Atom2}(\mathcal{Y}, \infty) \tag{57}$$

Then, using the definition of γ, we obtain the expression for θ<sup>1</sup> as:

$$\theta\_1 = \text{Atom2}(y, x) - \text{Atom2}(k\_1, k\_2) \tag{58}$$

As we have θ<sup>1</sup> and θ2, we can find θ<sup>3</sup> by using Eqs. (39) and (40)

$$\sum\_{i=1}^{3} \theta\_i = \text{Atom2}(s\rho, c\rho) = \rho = \theta\_1 + \theta\_2 + \theta\_3 \tag{59}$$

Finally, we can write other solution for θ<sup>2</sup> in terms of the tangent since we have

<sup>c</sup><sup>ψ</sup> <sup>¼</sup> <sup>x</sup><sup>2</sup> <sup>þ</sup> <sup>y</sup><sup>2</sup> <sup>þ</sup> <sup>l</sup><sup>1</sup>

β ¼ Atan2ð Þ y; x (64)

θ<sup>1</sup> ¼ β � ψ (66)

<sup>2</sup> <sup>¼</sup> <sup>x</sup><sup>2</sup> <sup>þ</sup> <sup>y</sup><sup>2</sup> <sup>þ</sup> <sup>l</sup><sup>1</sup>

θ<sup>i</sup> ¼ θ<sup>1</sup> þ θ<sup>2</sup> þ θ<sup>3</sup> ¼ φ (67)

<sup>6</sup>T (69)

(65)

<sup>2</sup> � <sup>2</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffiffi

<sup>x</sup><sup>2</sup> <sup>þ</sup> <sup>y</sup><sup>2</sup> <sup>p</sup> <sup>l</sup>1cos<sup>ψ</sup>

(68)

(70)

<sup>2</sup> � <sup>l</sup><sup>2</sup> 2

2l1l<sup>2</sup>

By solving this equation by cosψ, we obtain Eq. (65). Finally, we use Eqs. (64)

Solving Eq. (67) for θ3, we complete the solution. Note that φ is the specified orientation and we had already computed θ<sup>1</sup> and θ2. Thus, the problem is solved.

Let us consider a six-DOF revolute manipulator. The transformation from

a<sup>11</sup> a<sup>12</sup> a<sup>13</sup> px a<sup>21</sup> a<sup>22</sup> a<sup>23</sup> py a<sup>31</sup> a<sup>32</sup> a<sup>33</sup> pz 0 0 00

We know that this final matrix represents the following successive transforms

By premultiplying the left side of Eq. (68) by the inverse of Eq. (70), we obtain:

both sin of and cosine of the angle θ<sup>2</sup> Now Let us compute a solution for θ<sup>1</sup>

Geometry yields then

where β ¼ Atan2ð Þ y; x .

and (65), to obtain θ<sup>1</sup> as in Eq. (66). Now, we can compute θ<sup>3</sup> from:

frames {O} to the {N} (N = 6) is

For ψ, we use the law of cosines to obtain:

Kinematics for Spacecraft-Type Robotic Manipulators DOI: http://dx.doi.org/10.5772/intechopen.85636

ψ is obtained by using the law of cosines as ℓ<sup>2</sup>

∑ 3 i¼1

3.4 Algebraic solution method involving working matrix algebra

0 <sup>6</sup>T ¼

> 0 <sup>6</sup><sup>T</sup> <sup>¼</sup> <sup>0</sup> 1T<sup>1</sup> 2T<sup>2</sup> 3T<sup>3</sup> 4T<sup>4</sup> <sup>5</sup>T <sup>5</sup>

0 <sup>1</sup>T ¼

Where, for the revolute manipulator

121

Once have θ<sup>1</sup> and θ2, we obtain θ<sup>3</sup> from Eq. (59).

Transcendental equations commonly arise when applying the algebraic approach. The greater the number of degrees of freedom of the robotic manipulator, the more laborious is the solution of the problem.

#### 3.3 Geometric solution

By using this approach, we simply work the geometry of the arm. For the same problem, we may have the position of the joint associated to the frame {N} with respect to the origin of frame {O} by using geometry. Frame {N} refers to the frame attached to the last link of the manipulator. Consider Figure 6 showing the triangle formed by the geometrical configuration of the sides of the links and the line from the origin of the frame {O} to the origin of frame {N}. Applying the law of cosines, we can write:

$$\mathbf{x}^2 + \mathbf{y}^2 = l\_1^2 + l\_2^2 - 2l\_1l\_2 \cos\left(180 - \theta\_2\right) = l\_1^2 + l\_2^2 + 2l\_1l\_2c\_2\tag{60}$$

$$x\_2 = \frac{x^2 + y^2 - \left(l\_1^2 + l\_2^2\right)}{2l\_1l\_2} \tag{61}$$

Again, this solution must be between 1 and � 1 in order to have a valid solution. Assuming that the solution exists, then:

$$s\_2 = \sqrt{1 - \left[\frac{\varkappa^2 + \mathcal{y}^2 - \left(l\_1^{-2} + l\_2^{-2}\right)}{2l\_1l\_2}\right]^2} \tag{62}$$

Finally, we can write another solution for θ<sup>2</sup> in terms of the tangent since we have both sin and cosine of the angle θ<sup>2</sup>

$$\theta\_2 = \text{Atom2}\left\{ \pm \sqrt{1 - \left[ \frac{\mathbf{x}^2 + \mathbf{y}^2 - \left(l\_1^2 + l\_2^2\right)}{2l\_1l\_2} \right]^2}, \frac{\mathbf{x}^2 + \mathbf{y}^2 - \left(l\_1^2 + l\_2^2\right)}{2l\_1l\_2} \right\} \tag{63}$$

Figure 6. Space robot manipulator—geometric solutions for the inverse kinematics.

Finally, we can write other solution for θ<sup>2</sup> in terms of the tangent since we have both sin of and cosine of the angle θ<sup>2</sup>

Now Let us compute a solution for θ<sup>1</sup>

$$
\beta = \text{Atom2}(\mathfrak{y}, \mathfrak{x}) \tag{64}
$$

For ψ, we use the law of cosines to obtain:

$$c\psi = \frac{\varkappa^2 + \mathcal{y}^2 + l\_1^{\;\;2} - l\_2^{\;\;2}}{2l\_1l\_2} \tag{65}$$

Geometry yields then

As we have θ<sup>1</sup> and θ2, we can find θ<sup>3</sup> by using Eqs. (39) and (40)

Transcendental equations commonly arise when applying the algebraic approach. The greater the number of degrees of freedom of the robotic manipula-

By using this approach, we simply work the geometry of the arm. For the same problem, we may have the position of the joint associated to the frame {N} with respect to the origin of frame {O} by using geometry. Frame {N} refers to the frame attached to the last link of the manipulator. Consider Figure 6 showing the triangle formed by the geometrical configuration of the sides of the links and the line from the origin of the frame {O} to the origin of frame {N}. Applying the law of cosines, we can write:

<sup>2</sup> � <sup>2</sup>l1l<sup>2</sup> cos 180 ð Þ¼ � <sup>θ</sup><sup>2</sup> <sup>l</sup><sup>1</sup>

<sup>c</sup><sup>2</sup> <sup>¼</sup> <sup>x</sup><sup>2</sup> <sup>þ</sup> <sup>y</sup><sup>2</sup> � <sup>l</sup><sup>1</sup>

θ<sup>i</sup> ¼ Atan2ð Þ¼ sφ;cφ φ ¼ θ<sup>1</sup> þ θ<sup>2</sup> þ θ<sup>3</sup> (59)

<sup>2</sup> <sup>þ</sup> <sup>l</sup><sup>2</sup>

<sup>2</sup> <sup>þ</sup> <sup>l</sup><sup>2</sup> <sup>2</sup> � �

> <sup>2</sup> <sup>þ</sup> <sup>l</sup><sup>2</sup> <sup>2</sup> � �

vuut (62)

<sup>x</sup><sup>2</sup> <sup>þ</sup> <sup>y</sup><sup>2</sup> � <sup>l</sup><sup>1</sup>

2l1l<sup>2</sup>

<sup>2</sup> <sup>þ</sup> <sup>l</sup><sup>2</sup> <sup>2</sup> � �

9 = ;

2l1l<sup>2</sup>

Again, this solution must be between 1 and � 1 in order to have a valid solution.

<sup>1</sup> � <sup>x</sup><sup>2</sup> <sup>þ</sup> <sup>y</sup><sup>2</sup> � <sup>l</sup><sup>1</sup>

Finally, we can write another solution for θ<sup>2</sup> in terms of the tangent since we

<sup>2</sup> <sup>þ</sup> <sup>ℓ</sup><sup>2</sup> <sup>2</sup> � �

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

2l1l<sup>2</sup> " #<sup>2</sup>

<sup>1</sup> � <sup>x</sup><sup>2</sup> <sup>þ</sup> <sup>y</sup><sup>2</sup> � <sup>l</sup><sup>1</sup>

Space robot manipulator—geometric solutions for the inverse kinematics.

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

2l1l<sup>2</sup> " #<sup>2</sup>

,

<sup>2</sup> <sup>þ</sup> <sup>2</sup>l1l2c<sup>2</sup> (60)

(61)

(63)

∑ 3 i¼1

Kinematics - Analysis and Applications

3.3 Geometric solution

<sup>x</sup><sup>2</sup> <sup>þ</sup> <sup>y</sup><sup>2</sup> <sup>¼</sup> <sup>l</sup><sup>1</sup>

<sup>2</sup> <sup>þ</sup> <sup>l</sup><sup>2</sup>

s<sup>2</sup> ¼

Assuming that the solution exists, then:

have both sin and cosine of the angle θ<sup>2</sup>

vuut

8 < :

θ<sup>2</sup> ¼ Atan2 �

Figure 6.

120

Once have θ<sup>1</sup> and θ2, we obtain θ<sup>3</sup> from Eq. (59).

tor, the more laborious is the solution of the problem.

$$
\theta\_1 = \beta \pm \psi \tag{66}
$$

where β ¼ Atan2ð Þ y; x .

ψ is obtained by using the law of cosines as ℓ<sup>2</sup> <sup>2</sup> <sup>¼</sup> <sup>x</sup><sup>2</sup> <sup>þ</sup> <sup>y</sup><sup>2</sup> <sup>þ</sup> <sup>l</sup><sup>1</sup> <sup>2</sup> � <sup>2</sup> ffiffiffiffiffiffiffiffiffiffiffiffiffiffi <sup>x</sup><sup>2</sup> <sup>þ</sup> <sup>y</sup><sup>2</sup> <sup>p</sup> <sup>l</sup>1cos<sup>ψ</sup> By solving this equation by cosψ, we obtain Eq. (65). Finally, we use Eqs. (64) and (65), to obtain θ<sup>1</sup> as in Eq. (66).

Now, we can compute θ<sup>3</sup> from:

$$\sum\_{i=1}^{3} \theta\_i = \theta\_1 + \theta\_2 + \theta\_3 = \wp \tag{67}$$

Solving Eq. (67) for θ3, we complete the solution. Note that φ is the specified orientation and we had already computed θ<sup>1</sup> and θ2. Thus, the problem is solved.

#### 3.4 Algebraic solution method involving working matrix algebra

Let us consider a six-DOF revolute manipulator. The transformation from frames {O} to the {N} (N = 6) is

$$\begin{aligned} \, \_6^0T = \begin{bmatrix} a\_{11} & a\_{12} & a\_{13} & p\_x \\ a\_{21} & a\_{22} & a\_{23} & p\_y \\ a\_{31} & a\_{32} & a\_{33} & p\_x \\ 0 & 0 & 0 & 0 \end{bmatrix} \end{aligned} \tag{68}$$

We know that this final matrix represents the following successive transforms

$$T\_6^0 T = {}\_1^0 T\_2^1 T\_3^2 T\_4^3 T\_5^4 T\_6^5 T \tag{69}$$

Where, for the revolute manipulator

$$\begin{aligned} \, \_1^0T = \begin{bmatrix} c\_1 & -s\_1 & \mathbf{0} & \mathbf{0} \\ s\_1 & c\_1 & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{1} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{1} \end{bmatrix} \end{aligned} \tag{70}$$

By premultiplying the left side of Eq. (68) by the inverse of Eq. (70), we obtain:

$$
\begin{bmatrix} c\_1 & s\_1 & 0 & 0 \\ -s\_1 & c\_1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} a\_{11} & a\_{12} & a\_{13} & p\_x \\ a\_{21} & a\_{22} & a\_{23} & p\_y \\ a\_{31} & a\_{32} & a\_{33} & p\_x \\ 0 & 0 & 0 & 0 \end{bmatrix} = {}\_2^T T\_3^2 T\_4^3 T\_5^4 T\_6^5 T \tag{71}
$$

θ<sup>1</sup> ¼ �Atan2 px; py

Comparing elements from both sides yields

Kinematics for Spacecraft-Type Robotic Manipulators DOI: http://dx.doi.org/10.5772/intechopen.85636

then, we obtain from tan <sup>θ</sup><sup>3</sup> <sup>¼</sup> <sup>a</sup><sup>31</sup>

10 0 0 00 1 0 �1 �d<sup>3</sup> 0 00 0 1

The left side yields

compute θ<sup>3</sup> as:

transform;

Premultiplying the result of left side again by the inverse of <sup>1</sup>

a11cθ<sup>1</sup> þ a21sθ<sup>1</sup> a12cθ<sup>1</sup> þ a22sθ<sup>1</sup> a13cθ<sup>1</sup> þ a23sθ<sup>1</sup> pxcθ<sup>1</sup> þ pysθ<sup>1</sup> a<sup>31</sup> a<sup>32</sup> a<sup>33</sup> pz a11sθ<sup>1</sup> � a21cθ<sup>1</sup> a12sθ<sup>1</sup> � a22cθ<sup>1</sup> a13sθ<sup>1</sup> � a23cθ<sup>1</sup> pxsθ<sup>1</sup> � pycθ<sup>1</sup> � d<sup>2</sup> 000 1

By comparing the matrix elements (2,1) and (2,2) for both sides, we can

c<sup>1</sup> s<sup>1</sup> 0 0 �s<sup>1</sup> c<sup>1</sup> 0 0 0 010 0 001 � � (74)

pz ¼ 0 (76)

a<sup>31</sup> ¼ sθ<sup>3</sup> and a<sup>32</sup> ¼ cθ<sup>3</sup> (79)

<sup>a</sup><sup>32</sup> ¼)θ<sup>3</sup> ¼ Atan2ð Þ a31; a<sup>32</sup>

Of course, this example is simple. When we have a six-DOF manipulator, the algebraic manipulation to solve the inverse kinematics is hard, and the solution may require using transformation of variables. The MATLAB® software package has applications that easier the hard algebraic manipulation such as the symbolic manipulator and some specific functions for forward and inverse kinematics computation. MATLAB® functions include kinematics, trajectory generation, dynamics, and control. The toolbox also includes Simulink® models to describe the evolution of arm or mobile robot state over time for the sake of control. Regarding this chapter, functions like manipulating and converting between data types such as vectors, rotation matrices, homogeneous transformations, and twists are very

• angvec2r that converts angle and vector orientation to a 3x3 rotation matrix;

• angvec2tr that converts angle and vector orientation to a 4x4 homogeneous

• syms command that allows the definition of symbolic variables for rotation and homogeneous algebra. The symbolic manipulation allows include several

commands to simplify results of matrix algebraic manipulation.

homogeneous transformation matrix.

• eul2tr (phi, theta, psi, options) converts the Euler angles to a (4x4)

important. Some important functions related to this chapter are:

T, we obtain:

cθ<sup>3</sup> �sθ<sup>3</sup> 0 0 sθ<sup>3</sup> cθ<sup>3</sup> 0 0 sθ<sup>3</sup> cθ<sup>3</sup> 1 l<sup>2</sup> 0 0 01

(77)

(78)

d<sup>2</sup> ¼ pxsθ<sup>1</sup> � pycθ<sup>1</sup> � l<sup>2</sup> (75)

a<sup>11</sup> a<sup>12</sup> a<sup>13</sup> px a<sup>21</sup> a<sup>22</sup> a<sup>23</sup> py a<sup>31</sup> a<sup>32</sup> a<sup>33</sup> pz 0 0 00

By equating both sides of Eq. (71), we start identifying solutions by comparing the elements of matrices from both sides. Then, we continue by multiplying both sides, by the inverse of the <sup>1</sup> T and repeat the process to identify more solutions to obtain:

$$
\begin{bmatrix} c\_2 & s\_2 & 0 & -l\_1 c\_2 \\ -s\_2 & c\_2 & 0 & l\_1 s\_2 \\ \mathbf{0} & \mathbf{0} & \mathbf{1} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{1} \end{bmatrix} \begin{bmatrix} c\_1 & -s\_1 & \mathbf{0} & \mathbf{0} \\ s\_1 & c\_1 & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{1} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{1} \end{bmatrix} \begin{bmatrix} a\_{11} & a\_{12} & a\_{13} & p\_x \\ a\_{21} & a\_{22} & a\_{23} & p\_y \\ a\_{31} & a\_{32} & a\_{33} & p\_x \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} \end{bmatrix} = {}^3T\_4^3T\_5^4T\_6^5T\tag{72}
$$

By implementing this process, we can solve the inverse kinematics for the six-DOF revolute joints with all the z-axis parallels.

Consider Table 2 and the associated transformation matrix Eq. (36) for the RPR manipulator of Figure 5. Let us solve the inverse kinematics using the above method of matrix algebra. Al the matrix manipulations were computed by using the MATLAB® symbolic manipulator. The specified orientation and position is given by Eq. (68), by premultiplying the right side of this equation by inverse of Eq. (30), we obtain:

c1 s1 0 0 �s1 c1 0 0 0 0 10 0 001 a11 a12 a13 px a21 a22 a23 py a31 a32 a33 pz cθ<sup>3</sup> �sθ<sup>3</sup> 0 0 0 0 �1 �ð Þ l2 þ d2 sθ<sup>3</sup> cθ<sup>3</sup> 0 0 000 1 a11cθ<sup>1</sup> þ a21sθ<sup>1</sup> a12cθ<sup>1</sup> þ a22cs a13sθ<sup>1</sup> pxcθ<sup>1</sup> þ pysθ<sup>1</sup> a21cθ<sup>1</sup> � a11sθ<sup>1</sup> a22cθ<sup>1</sup> � a12cθ<sup>1</sup> a23cθ<sup>1</sup> � a13sθ<sup>1</sup> pycθ<sup>1</sup> � pxsθ<sup>1</sup> a31 a32 a33 pz 0 0 01 1 3T (73)

Comparing the resulting matrices, on the left and right, for the elements (1,4) and (2,4), respectively, we obtain:


Table 2. The D-H parameters and joint variables. Kinematics for Spacecraft-Type Robotic Manipulators DOI: http://dx.doi.org/10.5772/intechopen.85636

$$\theta\_1 = -\text{Atom2}\left(p\_x, p\_y\right) \tag{74}$$

$$d\_2 = p\_x \mathfrak{e} \theta\_1 - p\_y \mathfrak{e} \theta\_1 - l\_2 \tag{75}$$

Comparing elements from both sides yields

$$p\_x = 0\tag{76}$$

Premultiplying the result of left side again by the inverse of <sup>1</sup> T, we obtain:

$$
\begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & -1 & -d\_3 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} c\_1 & s\_1 & 0 & 0 \\ -s\_1 & c\_1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} a\_{11} & a\_{12} & a\_{13} & p\_x \\ a\_{21} & a\_{22} & a\_{23} & p\_y \\ a\_{31} & a\_{32} & a\_{33} & p\_x \\ 0 & 0 & 0 & 0 \end{bmatrix} = \begin{bmatrix} c\theta\_3 & -s\theta\_3 & 0 & 0 \\ s\theta\_3 & c\theta\_3 & 0 & 0 \\ s\theta\_3 & c\theta\_3 & 1 & l\_2 \\ 0 & 0 & 0 & 1 \end{bmatrix} \tag{77}
$$

The left side yields

c<sup>1</sup> s<sup>1</sup> 0 0 �s<sup>1</sup> c<sup>1</sup> 0 0 0 010 0 001

DOF revolute joints with all the z-axis parallels.

c<sup>1</sup> �s<sup>1</sup> 0 0 s<sup>1</sup> c<sup>1</sup> 0 0 0 0 10 0 0 01

a<sup>11</sup> a<sup>12</sup> a<sup>13</sup> px a<sup>21</sup> a<sup>22</sup> a<sup>23</sup> py a<sup>31</sup> a<sup>32</sup> a<sup>33</sup> pz 0 0 00

By equating both sides of Eq. (71), we start identifying solutions by comparing the elements of matrices from both sides. Then, we continue by multiplying both

By implementing this process, we can solve the inverse kinematics for the six-

Consider Table 2 and the associated transformation matrix Eq. (36) for the RPR

a11cθ<sup>1</sup> þ a21sθ<sup>1</sup> a12cθ<sup>1</sup> þ a22cs a13sθ<sup>1</sup> pxcθ<sup>1</sup> þ pysθ<sup>1</sup> a21cθ<sup>1</sup> � a11sθ<sup>1</sup> a22cθ<sup>1</sup> � a12cθ<sup>1</sup> a23cθ<sup>1</sup> � a13sθ<sup>1</sup> pycθ<sup>1</sup> � pxsθ<sup>1</sup> a31 a32 a33 pz 0 0 01

Comparing the resulting matrices, on the left and right, for the elements (1,4)

manipulator of Figure 5. Let us solve the inverse kinematics using the above method of matrix algebra. Al the matrix manipulations were computed by using the MATLAB® symbolic manipulator. The specified orientation and position is given by Eq. (68), by premultiplying the right side of this equation by inverse of Eq. (30), we

> a11 a12 a13 px a21 a22 a23 py a31 a32 a33 pz

 ¼ 1 T<sup>2</sup> T<sup>3</sup> T<sup>4</sup> T<sup>5</sup>

T and repeat the process to identify more solutions to

a<sup>11</sup> a<sup>12</sup> a<sup>13</sup> px a<sup>21</sup> a<sup>22</sup> a<sup>23</sup> py a<sup>31</sup> a<sup>32</sup> a<sup>33</sup> pz 0 0 00

 ¼ 2 T<sup>3</sup> T<sup>4</sup> T<sup>5</sup> T

cθ<sup>3</sup> �sθ<sup>3</sup> 0 0

sθ<sup>3</sup> cθ<sup>3</sup> 0 0 000 1

0 0 �1 �ð Þ l2 þ d2

T (71)

(72)

 1 3T

(73)

Kinematics - Analysis and Applications

sides, by the inverse of the <sup>1</sup>

c<sup>2</sup> s<sup>2</sup> 0 �l1c<sup>2</sup> �s<sup>2</sup> c<sup>2</sup> 0 l1s<sup>2</sup> 0 01 0 0 00 1

c1 s1 0 0 �s1 c1 0 0 0 0 10 0 001

and (2,4), respectively, we obtain:

The D-H parameters and joint variables.

obtain:

obtain:

Table 2.

$$\begin{bmatrix} a\_{11}c\theta\_1 + a\_{21}s\theta\_1 & a\_{12}c\theta\_1 + a\_{22}s\theta\_1 & a\_{13}c\theta\_1 + a\_{23}s\theta\_1 & p\_x c\theta\_1 + p\_y s\theta\_1\\ a\_{31} & a\_{32} & a\_{33} & p\_x\\ a\_{11}s\theta\_1 - a\_{21}c\theta\_1 & a\_{12}s\theta\_1 - a\_{22}c\theta\_1 & a\_{13}s\theta\_1 - a\_{23}c\theta\_1 & p\_x s\theta\_1 - p\_y c\theta\_1 - d\_2\\ 0 & 0 & 0 & 1 \end{bmatrix} \tag{78}$$

By comparing the matrix elements (2,1) and (2,2) for both sides, we can compute θ<sup>3</sup> as:

$$
\sigma\_{31} = \mathfrak{so}\_3 \text{ and } \mathfrak{a}\_{32} = \mathfrak{c}\theta\_3 \tag{79}
$$

then, we obtain from tan <sup>θ</sup><sup>3</sup> <sup>¼</sup> <sup>a</sup><sup>31</sup> <sup>a</sup><sup>32</sup> ¼)θ<sup>3</sup> ¼ Atan2ð Þ a31; a<sup>32</sup>

Of course, this example is simple. When we have a six-DOF manipulator, the algebraic manipulation to solve the inverse kinematics is hard, and the solution may require using transformation of variables. The MATLAB® software package has applications that easier the hard algebraic manipulation such as the symbolic manipulator and some specific functions for forward and inverse kinematics computation. MATLAB® functions include kinematics, trajectory generation, dynamics, and control. The toolbox also includes Simulink® models to describe the evolution of arm or mobile robot state over time for the sake of control. Regarding this chapter, functions like manipulating and converting between data types such as vectors, rotation matrices, homogeneous transformations, and twists are very important. Some important functions related to this chapter are:


• rpy2tr (roll, pitch, yaw, options) command refers to the LVLH frame referred in this chapter and converts the roll-pitch-yaw angles to homogeneous transform.

References

Congress. 2014

pp. 257-258

[1] Fonseca I, Mint PM. The state-of-theart in space robotics. Journal of Physics: Conference Series. 2015;641:7. DOI: 10.1088/1742-6596/641/1/012025

Kinematics for Spacecraft-Type Robotic Manipulators DOI: http://dx.doi.org/10.5772/intechopen.85636

> [9] Song W, Hu GM. A fast inverse kinematics algorithm for joint

24:350-354. DOI: 10.1016/j. proeng.2011.11.2655

[10] Corke P. Robotics, Vision & Control: Fundamental Algorithms in MATLAB. 2nd ed. Springer. 2011

animation. Procedia Engineering. 2011;

[2] Fonseca I, Pontuschka M, Saotome O, Bainum P, Lima G. The impact of the

operations of robotic manipulators. In: 65th International Astronautical

[3] Fonseca I, Goes L, Seito N, Duarte M, Mint OE. Attitude dynamics and control

manipulator when implementing onorbit servicing. Acta Astronautica. 2017;

[4] Skaar SB, Ruoff CF. Teleoperation and Robotics in Space, Progress in Astronautics and Aeronautics. Vol. 61. Washington DC: American Institute of Aeronautics and Astronautics, Inc; 1994.

[5] Craig J. Introduction to Robotics Mechanics and Control. 3rd ed. New York: Pearson; 2005. pp. 62-123. DOI:

10.13140/2.1.4486.1446

[6] Spong MW, Hutchinson S, Vidyasagar M. Robot Dynamics and Control. 2nd ed. ISBN: 978-0-

Sons, Inc; 2004. pp. 57-89

Institute of Aeronautics and

DOI: 10.2514/4.860119

125

471-64990-8. New York: John Wiley &

[7] Wie B. Space Vehicle Dynamics and Control. 2nd ed. Reston, VA: American

Astronautics, Inc; 2012. pp. 397-381.

[8] Zohar I, Ailon A, Mint RR. Mobile robot characterized by dynamic and kinematic equations and actuator dynamics: Trajectory tracking and related application. Robotics and Autonomous Systems. 2011;59:343-353. DOI: 10.1016/j.robot.2010.12.001

non-inertial base motion in the

of a spacecraft like a robotic

137:490-497. DOI: 10.1016/j. actaastro.2016.12.020

There are several other commands/functions that are very useful for applications approached in this chapter and the reader is advised to see [10].

### 4. Conclusion

In this chapter, we have approached the forward and inverse kinematics for space robot manipulators via the D-H convention. In all cases presented here, the spacecraft nominal attitude was considered stabilized. Otherwise, it would be necessary to consider kinematics in terms of absolute position described in the LVLH or in the spacecraft fixed frame with variable attitude angles. Some examples are solved, and a related reference is provided. Al the algebraic matrix manipulations presented here were obtained by using the MATLAB® Symbolic tool box. Also, the chapter presented other features of the MATLAB® and Simulink® software related to the subject of the chapter. Despite its importance, differential kinematics is out of scope of this chapter and maybe the subject of future work.

### Acknowledgements

Thanks to the Technological Institute of Aeronautics ITA/CTE and PG-EEC, National Institute for Space Research INPE/DMC, the Brazilian Space Agency (AEB), the Pontifical Catholic University of Sao Paulo PUC-SP, and FAPESP/PIPE —Project 00882-4/2017.

### Author details

Ijar Milagre da Fonseca1,2\*, Maurício Nacib Pontuschka3 and Glaydson Luiz Bertoze Lima<sup>1</sup>

1 Technological Institute of Aeronautics – DCTA/ITA, Sao Jose dos Campos, Brazil


\*Address all correspondence to: ijar@ita.br

© 2019 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Kinematics for Spacecraft-Type Robotic Manipulators DOI: http://dx.doi.org/10.5772/intechopen.85636

### References

• rpy2tr (roll, pitch, yaw, options) command refers to the LVLH frame

transform.

Kinematics - Analysis and Applications

4. Conclusion

Acknowledgements

—Project 00882-4/2017.

Author details

124

and Glaydson Luiz Bertoze Lima<sup>1</sup>

\*Address all correspondence to: ijar@ita.br

provided the original work is properly cited.

referred in this chapter and converts the roll-pitch-yaw angles to homogeneous

There are several other commands/functions that are very useful for applica-

In this chapter, we have approached the forward and inverse kinematics for space robot manipulators via the D-H convention. In all cases presented here, the spacecraft nominal attitude was considered stabilized. Otherwise, it would be necessary to consider kinematics in terms of absolute position described in the LVLH or in the spacecraft fixed frame with variable attitude angles. Some examples are solved, and a related reference is provided. Al the algebraic matrix manipulations presented here were obtained by using the MATLAB® Symbolic tool box. Also, the chapter presented other features of the MATLAB® and Simulink® software related to the subject of the chapter. Despite its importance, differential kinematics is out of

Thanks to the Technological Institute of Aeronautics ITA/CTE and PG-EEC, National Institute for Space Research INPE/DMC, the Brazilian Space Agency (AEB), the Pontifical Catholic University of Sao Paulo PUC-SP, and FAPESP/PIPE

1 Technological Institute of Aeronautics – DCTA/ITA, Sao Jose dos Campos, Brazil

2 National Institute for Space Research – DMC/INPE, Sao Jose dos Campos, Brazil

© 2019 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/ by/3.0), which permits unrestricted use, distribution, and reproduction in any medium,

3 Pontifical Catholic University of Sao Paulo PUC-SP, Sao Paulo, Brazil

tions approached in this chapter and the reader is advised to see [10].

scope of this chapter and maybe the subject of future work.

Ijar Milagre da Fonseca1,2\*, Maurício Nacib Pontuschka3

[1] Fonseca I, Mint PM. The state-of-theart in space robotics. Journal of Physics: Conference Series. 2015;641:7. DOI: 10.1088/1742-6596/641/1/012025

[2] Fonseca I, Pontuschka M, Saotome O, Bainum P, Lima G. The impact of the non-inertial base motion in the operations of robotic manipulators. In: 65th International Astronautical Congress. 2014

[3] Fonseca I, Goes L, Seito N, Duarte M, Mint OE. Attitude dynamics and control of a spacecraft like a robotic manipulator when implementing onorbit servicing. Acta Astronautica. 2017; 137:490-497. DOI: 10.1016/j. actaastro.2016.12.020

[4] Skaar SB, Ruoff CF. Teleoperation and Robotics in Space, Progress in Astronautics and Aeronautics. Vol. 61. Washington DC: American Institute of Aeronautics and Astronautics, Inc; 1994. pp. 257-258

[5] Craig J. Introduction to Robotics Mechanics and Control. 3rd ed. New York: Pearson; 2005. pp. 62-123. DOI: 10.13140/2.1.4486.1446

[6] Spong MW, Hutchinson S, Vidyasagar M. Robot Dynamics and Control. 2nd ed. ISBN: 978-0- 471-64990-8. New York: John Wiley & Sons, Inc; 2004. pp. 57-89

[7] Wie B. Space Vehicle Dynamics and Control. 2nd ed. Reston, VA: American Institute of Aeronautics and Astronautics, Inc; 2012. pp. 397-381. DOI: 10.2514/4.860119

[8] Zohar I, Ailon A, Mint RR. Mobile robot characterized by dynamic and kinematic equations and actuator dynamics: Trajectory tracking and related application. Robotics and Autonomous Systems. 2011;59:343-353. DOI: 10.1016/j.robot.2010.12.001

[9] Song W, Hu GM. A fast inverse kinematics algorithm for joint animation. Procedia Engineering. 2011; 24:350-354. DOI: 10.1016/j. proeng.2011.11.2655

[10] Corke P. Robotics, Vision & Control: Fundamental Algorithms in MATLAB. 2nd ed. Springer. 2011

## *Edited by Joseph Mizrahi*

Numerous problems in engineering and biology can be described, characterized, and analyzed in kinematics terms. In classical machinery and robotics the most distinctive characteristic is constrained motion of multi-degree-of-freedom kinematic chains. Robotic arms and manipulators have become essential devices in industrial applications and medicine. This book provides the reader with an updated look at the current trends in kinematics methods and applications. Section 1 deals with kinematics of linkages and includes analysis of cam mechanisms and transformation of rotary motion into oscillation. Section 2 covers compliant mechanisms, whereby elastically deformable parts are part of the mechanism. Finally, Section 3 deals with kinematics of spacecrafts and satellites in the contexts of global navigation systems and of space robot analysis.

Published in London, UK © 2018 IntechOpen © deversteel / iStock

Kinematics - Analysis and Applications

Kinematics

Analysis and Applications

*Edited by Joseph Mizrahi*