# A Novel Comprehensive Model of Suitability Analysis for Matching Area in Underwater Geomagnetic Aided Inertial Navigation

Authors: Hongmei Zhang; Le Yang; Minglong Li
Journal: Mathematical Problems in Engineering
 (2019)

---

## Abstract

Geomagnetic aided inertial navigation is a way to use the geophysical field for navigation. It can locate the carrier position by the correlation between geomagnetic data and running track. It is an effective mean to realize autonomous navigation. Matching area suitability is one of the important factors affecting geomagnetic aided inertial navigation. Through the suitability analysis of matching areas, the areas with obvious geomagnetic features and rich information are selected as matching areas, which can effectively improve the real-time and accuracy of geomagnetic aided navigation. However, matching area suitability analysis for geomagnetic aided navigation is a complex process and needs to consider diverse factors, based on which a decision may be made. The area suitability analysis inherently can be considered as a multicriterion decision analysis (MCDA) problem. This paper presented a novel comprehensive model combining principal component analysis (PCA) and analytical hierarchy process (AHP) to evaluate the suitability of the navigation matching area. Firstly, according to the features of the areas, key feature parameters and the corresponding weights are determined by PCA and AHP, respectively. Then comprehensive evaluation values of the navigation matching areas are calculated through the comprehensive model. Finally, experiments were implemented in Bohai Bay; the correlation-matching algorithm is applied to verify the validity of the model in the areas. The experiment results well indicate the consistency between the comprehensive evaluation value and the matching area suitability. It is reasonable to regard the comprehensive evaluation value as a basis for area suitability analysis.


---

## Body

## 1. Introduction

Geomagnetic aided navigation (GAN) has become a research focus in recent years because of its merits such as all-weather, passive, no-radiation, error bounded, and strong anti-interference ability [1]. The geomagnetic aided inertial navigation system is one of the GAN systems. The principle of geomagnetic aided inertial navigation system is that when a vessel enters the matching areas, a matching algorithm is used to obtain the position information; this information can correct the cumulative error of inertial navigation system (INS) in time to guarantee that the final integrated navigation system output is credible [2–4].Since each matching area includes different geomagnetic information, the matching effect is different in each area. Selecting the areas that have rich geomagnetic features and good suitability can help improve navigation accuracy and stability, and optimize navigation path planning as well. The area suitability analysis is one of the key technologies for GAN, the objective of which is to identify the areas with the highest potential for meeting the needs of navigation consistently at an acceptable cost [5].Suitability is an intrinsic property of the geomagnetic matching area. According to the published literature, there are many features describing the suitability of geomagnetic field, such as roughness, gradient deviation, information entropy, standard deviation, and correlation coefficient [6–8]. Many methods for evaluating the suitability of navigation matching areas based on these features have been proposed. Kang took geomagnetic entropy and geomagnetic variance entropy as reference standards of geomagnetic matching area [9]. In view of the problems of one-sided appraisal when using single feature parameter to evaluate geomagnetic map suitability, some multiattribute decision-making methods have been proposed. Zhu proposed a comprehensive evaluation method using entropy technology to enhance the weight, thus overcoming the poor objectivity of traditional fuzzy evaluation method when confirming the weight of index [10]. Li proposed a method for evaluating geomagnetic matching area by synthesizing several geomagnetic features [11]. All the above methods aggregate the feature parameters, but they are limited to the fusion of several parameters with abundant information. However, in practical application, some parameters’ information is not rich but can still play a positive role in the evaluation of area suitability. Therefore, it is particularly important to consider the feature parameters extensively. In such a context, the idea of multicriterion decision analysis (MCDA) can be introduced to address this problem. The principle of MCDA is to aggregate decision-making information in a certain way, and then the schemes are sorted and optimized to provide reference for decision-makers (DMS) [12–14].Among all the feature parameters for the matching area in GAN, some have great influence on the level of suitability. If these parameters can be qualitatively distinguished, that will be beneficial to the weight decision analysis. Principal component analysis is a commonly used multivariate analysis method; its core idea is dimension reduction. This feature of PCA is very meaningful when there are many analysis indicators, especially considering the large amount of calculation. Through principal component analysis, the analysis process and amount of calculations can be greatly reduced. Due to the advantages such as objectivity, simple calculation, and convenient application, PCA has been widely used in mathematical modeling, mathematical analysis, and other disciplines [15–17]. Li [15] introduced PCA into geomagnetic navigation matching area selection, which took each principal component as the object, and the contribution rate was used as the weight to obtain the comprehensive evaluation value.Analytical hierarchy process (AHP) is a common method in MCDA. In recent years, AHP has been widely used to determine weights [18, 19] and thus can be used to perform weight decision analysis on matching area parameters. The advantage of AHP is that it can provide a quantitative reference for qualitative analysis of human subjective judgments. The disadvantage is that when constructing the judgment matrix, the importance between the parameters depends too much on empirical judgment, not objective enough, so the AHP cannot be directly used to calculate the evaluation value of the original data. It is with this in mind that Wang combined accelerated genetic algorithm and AHP to extract the ranking information of evaluation index sample set and determine the ranking weight of each evaluation index, which could obtain more accurate parameter weights [20].In addition, there are some other weighting methods, such as Delphi method, entropy method, and coefficient of variation method [21–23]. The Delphi method relies entirely on the judgment of experts, and sometimes it is one-sided and difficult to review whether it is correct or not. The entropy method and coefficient of variation method are objective valuation methods, which can avoid the deviation caused by human factors. However, the entropy method neglects the importance of the index itself, and sometimes the weight of the index will differ greatly from the expected value. The coefficient of variation method does not pay enough attention to the specific significance of the index, and there will be some errors. Although AHP is also a subjective weighting method, its disadvantage is that the construction of evaluation matrix is not objective enough, but when combined with PCA, the PCA-AHP model can provide an objective reference. Therefore, based on the above analysis, the combination of PCA and AHP is more suitable in this case.This paper combine the principal component analysis and the analytic hierarchy process to establish the PCA-AHP model, and then all parameters in the matching areas are analyzed and fused by the model, which ensures the integrity of the geomagnetic parameters. After the data in the navigation area are processed by PCA, the load coefficients of each principal component can provide objective basis for AHP to construct the judgment matrix, thereby obtaining more accurate weight information.The structure of this paper is as follows. First, the common geomagnetic feature parameters in the matching areas are described. Then the principle of PCA and AHP is elaborated upon. After that, the PCA-AHP model is proposed. Based on the above, a comprehensive evaluation value of the matching area is calculated. Actual applications verify the proposed method, performances of the proposed method are also discussed, and some conclusions and suggestions are drawn from the experiments and discussion.


## 2. Geomagnetic Feature Parameters

The geomagnetic field data is stored in the computer in the form of a grid. The features extracted from the geomagnetic field data are called geomagnetic features. Suitability is an intrinsic property of geomagnetic fields. Since geomagnetic fields are expressed in digital geomagnetic maps, this property must be reflected by the geomagnetic grid content and statistical characteristics. Set the size of the geomagnetic map toM×N (M is the length in the latitudinal direction, N is the length in the longitudinal direction), and f(i,j) is the geomagnetic field value of coordinate (i,j), where (i,j) corresponds to a pair of latitude and longitude coordinates φ,λ. This paper considers several geomagnetic feature parameters as follows.


### 2.1. Average Value

Average valuef¯ represents the average of the magnetic field in the candidate matching areas.(1)f¯=1MN∑i=1M∑j=1Nfi,j.


### 2.2. Standard Deviation

The geomagnetic standard deviation reflects the discretization and the general fluctuation of the geomagnetic field. The larger the standard deviation, the more obvious the fluctuation of the geomagnetic field, which is more favorable for geomagnetic matching. The standard deviationδ is defined as(2)δ=1MN-1∑i=1M∑j=1Nfi,j-f¯2.


### 2.3. Roughness

Geomagnetic roughness reflects the average smooth level and local heave of geomagnetic field in certain area. The bigger the roughness, the richer the geomagnetic information. The roughnessr is defined as(3)r=rx+ry2,where(4)rx=1MN-1∑i=1M∑j=1N-1fi,j-fi,j+12,(5)ry=1M-1N∑i=1M-1∑j=1Nfi,j-fi+1,j2.Here rx denotes the roughness in latitudinal direction and ry denotes the roughness in longitudinal direction.


### 2.4. Kurtosis Coefficient

Kurtosis coefficient reflects the data concentration level. The bigger the kurtosis coefficient, the higher the concentration of the data near the average value, and the more difficult for matching; conversely, the more uniform the distribution, the easier for matching. The kurtosis coefficientCe is defined as(6)Ce=1MN∑i=1M∑j=1Nfi,j-f¯4δ4-3.


### 2.5. Skewness Coefficient

Skewness coefficient reflects the symmetry or skewness of the geomagnetic field, the bigger the value, the higher the asymmetry of the data, and the simpler the matching along the asymmetric direction; the smaller the value, the higher the symmetry of the date, and the easier the mismatching along the symmetric direction. The skewness coefficientCs is defined as(7)Cs=1MN∑i=1M∑j=1Nfi,j-f¯3δ3.


### 2.6. Information Entropy

Information entropy, also known as Shannon entropy, is proposed by Shannon to solve the problem of quantitative measurement of information [8]. Geomagnetic information entropy is used to measure the amount of geomagnetic information. The information entropy H is defined as(8)H=-∑i=1M∑j=1Npi,jlog2⁡pi,j,where(9)pi,j=fi,j∑i=1M∑j=1Nfi,j.


## 2.1. Average Value

Average valuef¯ represents the average of the magnetic field in the candidate matching areas.(1)f¯=1MN∑i=1M∑j=1Nfi,j.


## 2.2. Standard Deviation

The geomagnetic standard deviation reflects the discretization and the general fluctuation of the geomagnetic field. The larger the standard deviation, the more obvious the fluctuation of the geomagnetic field, which is more favorable for geomagnetic matching. The standard deviationδ is defined as(2)δ=1MN-1∑i=1M∑j=1Nfi,j-f¯2.


## 2.3. Roughness

Geomagnetic roughness reflects the average smooth level and local heave of geomagnetic field in certain area. The bigger the roughness, the richer the geomagnetic information. The roughnessr is defined as(3)r=rx+ry2,where(4)rx=1MN-1∑i=1M∑j=1N-1fi,j-fi,j+12,(5)ry=1M-1N∑i=1M-1∑j=1Nfi,j-fi+1,j2.Here rx denotes the roughness in latitudinal direction and ry denotes the roughness in longitudinal direction.


## 2.4. Kurtosis Coefficient

Kurtosis coefficient reflects the data concentration level. The bigger the kurtosis coefficient, the higher the concentration of the data near the average value, and the more difficult for matching; conversely, the more uniform the distribution, the easier for matching. The kurtosis coefficientCe is defined as(6)Ce=1MN∑i=1M∑j=1Nfi,j-f¯4δ4-3.


## 2.5. Skewness Coefficient

Skewness coefficient reflects the symmetry or skewness of the geomagnetic field, the bigger the value, the higher the asymmetry of the data, and the simpler the matching along the asymmetric direction; the smaller the value, the higher the symmetry of the date, and the easier the mismatching along the symmetric direction. The skewness coefficientCs is defined as(7)Cs=1MN∑i=1M∑j=1Nfi,j-f¯3δ3.


## 2.6. Information Entropy

Information entropy, also known as Shannon entropy, is proposed by Shannon to solve the problem of quantitative measurement of information [8]. Geomagnetic information entropy is used to measure the amount of geomagnetic information. The information entropy H is defined as(8)H=-∑i=1M∑j=1Npi,jlog2⁡pi,j,where(9)pi,j=fi,j∑i=1M∑j=1Nfi,j.


## 3. PCA and AHP Model

### 3.1. The Principle of Principal Component Analysis

Principal component analysis (PCA) aims to reveal the intrinsic relationship between multifeatures and large sample size by using mathematical dimension reduction ideas. PCA can be used to transform multivariate date into a few uncorrelated synthetic variables and reduce the dimension of observation space, while retaining the most important information. In the progress of multicriterion decision analysis, the complexity of decision-making will increase greatly if the number of variables is large, so it is best to use as few variables as possible. In many cases, there is a certain correlation between the variables. When there is a certain correlation between the two variables, it can be explained that the information reflected by the two variables has a certain overlap. PCA is about excluding redundant variables (closely related variables) and creates as few new variables as possible, so that these new variables are uncorrelated, and these new variables keep as much as possible of the original information.The detailed steps are as follows:(1)
Establish the evaluation matrix according to the geomagnetic field feature parameters of the matching area (evaluation index). Suppose the number of matching area ism and the evaluation index is n; then the evaluation matrix is(10)X=x11…x1n⋮⋱⋮xm1⋯xmn=xijm×ni=1,2,…,m;j=1,2,…,n,where xij represents the evaluation value of thejth indicator of theith matching area.(2)
Normalize the evaluation matrix. Different evaluation indexes often have different dimensions and dimension units, which affect the results of data analysis. The purpose of normalization is to eliminate the dimension effect between indexes and to solve the comparability between data indexes. After data normalization, the original data are in the same order of magnitude, which is suitable for comprehensive comparison and evaluation.(11)zij=xij-xj¯σj,(12)xj¯=1m∑i=1mxij,(13)σj=1m∑i=1mxij-xj¯2,where xj¯ is the average value of jth index and σj is the variance of the jth index. The normalized evaluation matrix is(14)Z=zijm×ni=1,2,…,m;j=1,2,…,n.(3)
Calculate the correlation coefficient matrixR.(15)R=rkjn×nk,j=1,2,…,n.(4)
Solve the eigenvalues and corresponding eigenvectors of correlation coefficient matrixR. Suppose that λ1≥λ2≥⋯≥λn≥0 are the eigenvalues of R and Lg1,Lg2,…,Lgi are the corresponding eigenvectors. Thus, the principal component can be obtained, which is represented by Fi.(16)Fi=∑j=1nLgiziji=1,2,…,n.(5)
Determine the number of principal components. In order to reduce the workload and the loss of information as much as possible, only the firstk principal components are kept. The k value can be determined through the cumulative contribution rate α(k), and the criteria are as follows:(17)αg=λg∑g=1nλg,(18)αk=∑g=1kαg=∑g=1kλg∑g=1nλg.The greater the contribution rate, the stronger the information about the original variable contained in the principal component. Generally, thek principal components whose cumulative contribution rate is above 80% and whose eigenvalues are greater than 1 are taken as the principal components for the final selection.(6)
After determining the number of principal components, the eigenvectors corresponding to each principal component eigenvalue are calculated according to the eigenvalues of the correlation coefficient matrix, and the influence of each index on the principal component can be known according to the descending order of the eigenvector coefficients. The larger the value, the greater the influence of the indicator on the principal component. If the coefficient is less than 0.1, its influence on the principal component is insignificantly small and can be neglected; the index represented by the coefficient can be eliminated.According to the above steps, a complete index system that can comprehensively reflect the information contained in all indexes can be determined and can eliminate the correlation between the indexes and facilitate the operation of AHP.


### 3.2. The Principle of Analytical Hierarchy Process

The analytic hierarchy process (AHP) was put forward by Saaty in the middle of 1970s [24]. It is a combination of a qualitative and quantitative, systematic and hierarchical analysis method. Because of its practicality and effectiveness in dealing with complex decision-making issues, it soon gained worldwide attention. The principle of AHP is to decompose the problem according to its feature parameters and the ultimate goal to be achieved, and then classify and combine the factors at different levels according to their interrelated effects and membership relations, so as to form a structural model with multilevel analysis. Finally, the problem is attributed to the determination of the relative importance weights of the lowest level relative to the highest level or the arrangement of merits and demerits.The detailed steps are as follows.(1) Establish a Hierarchical Structure Model. The general structural model is divided into three layers, the top is the target layer, the bottom is the plan layer, and the middle is the criterion layer.(2) Construct All the Judgment Matrixes in Each Level. The judgment matrix is constructed by comparing the weights of the criteria with respect to the target by mutual comparison. In the analytic hierarchy process, in order to quantify the importance of each element in the matrix, a matrix judgment scale (1-9 scale) is introduced, as shown in Table 1.Table 1
Matrix judgment scale.
Scale Meaning 1 Equal importance  3 Weak importance  5 Essential or strong importance  7 Very strong importance  9 Absolute importance  2,4,6,8 The intermediate values of the adjacent judgments mentioned above If the importance ratio of elementi to element j is aij, then the importance ratio of element j to element i is aji=1/aij.(3) Hierarchical Single Ranking. Hierarchical single ranking refers to the importance ranking of the elements at the current level for a certain element on the previous level. The relative weights of the compared elements for the current level elements are calculated by the judgement matrix.For the judgement matrixA, calculate the eigenvalues and eigenvectors that satisfy (19)AW=λmax′W,where λ′max is the largest eigenvalue of matrixA andW is the normalized eigenvector corresponding to λ′max. The component ωi ofW is the weight of the corresponding elements in single order. The detailed steps are as follows:(a)
Normalize each column vector ofA.(20)ωij∨=aij∑i=1naij.(b)
Sum each row ofωij∨.(21)ωi∨=∑j=1nωij∨.(c)
Normalizeωi∨.(22)ωi=ωi∨∑i=1nωi∨,(23)ω=ω1,ω2,…ωnT.(d)
λ′max can be calculated from (24)λ′max=1n∑i=1nAωiωi.(4)
Consistency index calculation:Calculating the weight vector from the judgment matrix requires that the judgment matrix has a general consistency, and the consistency indexCI needs to be calculated by (25)CI=λ′max-nn-1.In order to test whether the judgment matrix has satisfactory consistency, it is necessary to find the standard for measuring the consistency indexCI of the matrix A, and Saaty introduced the random consistency index RI, as shown in Table 2 [16].Table 2
Random consistency index.
n 1 2 3 4 5 6 7 8 9 10 11 RI 0 0 0.58 0.90 1.12 1.24 1.32 1.41 1.45 1.49 1.51Consistency ratioCR is defined as (26)CR=CIRI.When CR<0.1, the degree of inconsistency of A is within the allowable range and the eigenvector of A can be used as the weight vector.(5)
Hierarchical total ranking.Hierarchical total ranking is a process of ranking weights to determine the relative importance of all elements in a layer to the overall objective. Suppose the rank of theM elements B1,B2,…,Bm in layer B to total objective is b1,b2,…,bm; the hierarchical single ranking of the N elements in layer C to layer B is c1j,c2j,…,bnjj=1,2,…,m. Then the hierarchical total ranking of layer C is(27)C1:b1c11+b2c12+⋯bmc1mC2:b1c21+b2c22+⋯bmc2m⋯Cn:b1cn1+b2cn2+⋯bmcnm.That is, the weight of the ith element of the C layer to the total objective is(28)Wi¯=∑j=1mbjcij.Suppose that the hierarchical single ranking consistency index of the factors in the C-layer to the B-layer isCIj, and the random consistency index is RIj; then the total order consistency ratio is(29)CR′=a1CI1+a2CI2+⋯amCIma1RI1+a2RI2+⋯amRIm.When CR′<0.1, it means that the hierarchical total ranking passes the consistency check.


### 3.3. Establishment of PCA-AHP Model

The purpose of establishing the PCA-AHP model is to get a comprehensive evaluation value of navigation area suitability. The flowchart of the model is shown in Figure1.Figure 1
PCA-AHP model.There are four main methods to synthesize multi-index: linear synthesis method, geometric synthesis method, mixed synthesis method, and model synthesis method. Simply speaking, they are weighted summation, weighted geometric average, linear weighting, and geometric synthesis [25–28]. In this paper, the comprehensive evaluation values of the matching areas can be determined by (30)Ti=∑j=1nzijW¯i.The larger Ti is, the better the suitability performance of the ith candidate matching area is. According to this, each matching area can be sorted.


## 3.1. The Principle of Principal Component Analysis

Principal component analysis (PCA) aims to reveal the intrinsic relationship between multifeatures and large sample size by using mathematical dimension reduction ideas. PCA can be used to transform multivariate date into a few uncorrelated synthetic variables and reduce the dimension of observation space, while retaining the most important information. In the progress of multicriterion decision analysis, the complexity of decision-making will increase greatly if the number of variables is large, so it is best to use as few variables as possible. In many cases, there is a certain correlation between the variables. When there is a certain correlation between the two variables, it can be explained that the information reflected by the two variables has a certain overlap. PCA is about excluding redundant variables (closely related variables) and creates as few new variables as possible, so that these new variables are uncorrelated, and these new variables keep as much as possible of the original information.The detailed steps are as follows:(1)
Establish the evaluation matrix according to the geomagnetic field feature parameters of the matching area (evaluation index). Suppose the number of matching area ism and the evaluation index is n; then the evaluation matrix is(10)X=x11…x1n⋮⋱⋮xm1⋯xmn=xijm×ni=1,2,…,m;j=1,2,…,n,where xij represents the evaluation value of thejth indicator of theith matching area.(2)
Normalize the evaluation matrix. Different evaluation indexes often have different dimensions and dimension units, which affect the results of data analysis. The purpose of normalization is to eliminate the dimension effect between indexes and to solve the comparability between data indexes. After data normalization, the original data are in the same order of magnitude, which is suitable for comprehensive comparison and evaluation.(11)zij=xij-xj¯σj,(12)xj¯=1m∑i=1mxij,(13)σj=1m∑i=1mxij-xj¯2,where xj¯ is the average value of jth index and σj is the variance of the jth index. The normalized evaluation matrix is(14)Z=zijm×ni=1,2,…,m;j=1,2,…,n.(3)
Calculate the correlation coefficient matrixR.(15)R=rkjn×nk,j=1,2,…,n.(4)
Solve the eigenvalues and corresponding eigenvectors of correlation coefficient matrixR. Suppose that λ1≥λ2≥⋯≥λn≥0 are the eigenvalues of R and Lg1,Lg2,…,Lgi are the corresponding eigenvectors. Thus, the principal component can be obtained, which is represented by Fi.(16)Fi=∑j=1nLgiziji=1,2,…,n.(5)
Determine the number of principal components. In order to reduce the workload and the loss of information as much as possible, only the firstk principal components are kept. The k value can be determined through the cumulative contribution rate α(k), and the criteria are as follows:(17)αg=λg∑g=1nλg,(18)αk=∑g=1kαg=∑g=1kλg∑g=1nλg.The greater the contribution rate, the stronger the information about the original variable contained in the principal component. Generally, thek principal components whose cumulative contribution rate is above 80% and whose eigenvalues are greater than 1 are taken as the principal components for the final selection.(6)
After determining the number of principal components, the eigenvectors corresponding to each principal component eigenvalue are calculated according to the eigenvalues of the correlation coefficient matrix, and the influence of each index on the principal component can be known according to the descending order of the eigenvector coefficients. The larger the value, the greater the influence of the indicator on the principal component. If the coefficient is less than 0.1, its influence on the principal component is insignificantly small and can be neglected; the index represented by the coefficient can be eliminated.According to the above steps, a complete index system that can comprehensively reflect the information contained in all indexes can be determined and can eliminate the correlation between the indexes and facilitate the operation of AHP.


## 3.2. The Principle of Analytical Hierarchy Process

The analytic hierarchy process (AHP) was put forward by Saaty in the middle of 1970s [24]. It is a combination of a qualitative and quantitative, systematic and hierarchical analysis method. Because of its practicality and effectiveness in dealing with complex decision-making issues, it soon gained worldwide attention. The principle of AHP is to decompose the problem according to its feature parameters and the ultimate goal to be achieved, and then classify and combine the factors at different levels according to their interrelated effects and membership relations, so as to form a structural model with multilevel analysis. Finally, the problem is attributed to the determination of the relative importance weights of the lowest level relative to the highest level or the arrangement of merits and demerits.The detailed steps are as follows.(1) Establish a Hierarchical Structure Model. The general structural model is divided into three layers, the top is the target layer, the bottom is the plan layer, and the middle is the criterion layer.(2) Construct All the Judgment Matrixes in Each Level. The judgment matrix is constructed by comparing the weights of the criteria with respect to the target by mutual comparison. In the analytic hierarchy process, in order to quantify the importance of each element in the matrix, a matrix judgment scale (1-9 scale) is introduced, as shown in Table 1.Table 1
Matrix judgment scale.
Scale Meaning 1 Equal importance  3 Weak importance  5 Essential or strong importance  7 Very strong importance  9 Absolute importance  2,4,6,8 The intermediate values of the adjacent judgments mentioned above If the importance ratio of elementi to element j is aij, then the importance ratio of element j to element i is aji=1/aij.(3) Hierarchical Single Ranking. Hierarchical single ranking refers to the importance ranking of the elements at the current level for a certain element on the previous level. The relative weights of the compared elements for the current level elements are calculated by the judgement matrix.For the judgement matrixA, calculate the eigenvalues and eigenvectors that satisfy (19)AW=λmax′W,where λ′max is the largest eigenvalue of matrixA andW is the normalized eigenvector corresponding to λ′max. The component ωi ofW is the weight of the corresponding elements in single order. The detailed steps are as follows:(a)
Normalize each column vector ofA.(20)ωij∨=aij∑i=1naij.(b)
Sum each row ofωij∨.(21)ωi∨=∑j=1nωij∨.(c)
Normalizeωi∨.(22)ωi=ωi∨∑i=1nωi∨,(23)ω=ω1,ω2,…ωnT.(d)
λ′max can be calculated from (24)λ′max=1n∑i=1nAωiωi.(4)
Consistency index calculation:Calculating the weight vector from the judgment matrix requires that the judgment matrix has a general consistency, and the consistency indexCI needs to be calculated by (25)CI=λ′max-nn-1.In order to test whether the judgment matrix has satisfactory consistency, it is necessary to find the standard for measuring the consistency indexCI of the matrix A, and Saaty introduced the random consistency index RI, as shown in Table 2 [16].Table 2
Random consistency index.
n 1 2 3 4 5 6 7 8 9 10 11 RI 0 0 0.58 0.90 1.12 1.24 1.32 1.41 1.45 1.49 1.51Consistency ratioCR is defined as (26)CR=CIRI.When CR<0.1, the degree of inconsistency of A is within the allowable range and the eigenvector of A can be used as the weight vector.(5)
Hierarchical total ranking.Hierarchical total ranking is a process of ranking weights to determine the relative importance of all elements in a layer to the overall objective. Suppose the rank of theM elements B1,B2,…,Bm in layer B to total objective is b1,b2,…,bm; the hierarchical single ranking of the N elements in layer C to layer B is c1j,c2j,…,bnjj=1,2,…,m. Then the hierarchical total ranking of layer C is(27)C1:b1c11+b2c12+⋯bmc1mC2:b1c21+b2c22+⋯bmc2m⋯Cn:b1cn1+b2cn2+⋯bmcnm.That is, the weight of the ith element of the C layer to the total objective is(28)Wi¯=∑j=1mbjcij.Suppose that the hierarchical single ranking consistency index of the factors in the C-layer to the B-layer isCIj, and the random consistency index is RIj; then the total order consistency ratio is(29)CR′=a1CI1+a2CI2+⋯amCIma1RI1+a2RI2+⋯amRIm.When CR′<0.1, it means that the hierarchical total ranking passes the consistency check.


## 3.3. Establishment of PCA-AHP Model

The purpose of establishing the PCA-AHP model is to get a comprehensive evaluation value of navigation area suitability. The flowchart of the model is shown in Figure1.Figure 1
PCA-AHP model.There are four main methods to synthesize multi-index: linear synthesis method, geometric synthesis method, mixed synthesis method, and model synthesis method. Simply speaking, they are weighted summation, weighted geometric average, linear weighting, and geometric synthesis [25–28]. In this paper, the comprehensive evaluation values of the matching areas can be determined by (30)Ti=∑j=1nzijW¯i.The larger Ti is, the better the suitability performance of the ith candidate matching area is. According to this, each matching area can be sorted.


## 4. Experiment and Analysis

An experiment was implemented in a water area in Bohai Bay. An underwater reference for simulation is illustrated in Figure2. The size of the geomagnetic map is about 20km∗10km. Through the geomagnetic background field modeling, the geomagnetic map with ±5nT accuracy is obtained. The PCA and PCA-AHP were used to obtain the comprehensive evaluation values of each matching area, and then the accuracy was verified by experiments.Figure 2
Underwater reference geomagnetic map for simulation.In order to test the methods mentioned above, seven candidate matching areas of size 15×15 (grid precision: 200m) are selected, which are marked with red squares in Figure2. The feature parameters of the candidate matching areas are calculated and listed in Table 3.Table 3
The parameters in the candidate matching areas.
Area Parameters Standard Deviation/nT x-Roughness/nT y-Roughness/nT Roughness/nT Entropy/bit Kurtosis Coefficient Skewness Coefficient 1 18.513 1.204 2.203 2.511 7.053 -1.454 0.318  2 9.886 1.273 1.233 1.772 7.053 -0.776 -0.298  3 30.721 1.776 2.488 3.057 7.082 -0.821 -0.091  4 46.993 1.907 4.366 4.764 7.053 -1.044 0.041  5 34.242 1.211 3.256 3.474 7.053 -1.077 -0.073  6 21.615 1.204 1.749 2.123 7.082 -0.122 0.547  7 13.303 0.993 1.522 1.817 7.082 -0.345 0.348


### 4.1. Principal Component Analysis

According to the (10)-(15), the normalized evaluation matrix Z and the correlation matrix R can be obtained.(31)Z=-0.500-0.482-0.182-0.258-0.802-1.4310.685-1.162-0.278-1.067-0.945-0.8020.065-1.3760.4361.2110.0780.2501.069-0.034-0.6831.6831.5991.7921.837-0.802-0.526-0.2410.706-0.4610.7790.638-0.802-0.599-0.623-0.263-0.482-0.596-0.6191.0691.5091.452-0.900-1.106-0.804-0.9031.0691.0170.786.(32)R=1.0000.7630.9690.975-0.227-0.332-0.1630.7631.0000.6940.771-0.118-0.298-0.4070.9690.6941.0000.993-0.412-0.490-0.1870.9750.7710.9931.000-0.397-0.484-0.249-0.227-0.118-0.412-0.3971.0000.7770.485-0.332-0.298-0.490-0.4840.7771.0000.407-0.163-0.407-0.187-0.2490.4850.4071.000.Solve the eigenvaluesλn and eigenvector Lgi of the correlation matrixR and calculate the contribution rate α(g) and cumulative contribution rate α(k) by formula (17) and (18), as listed in Table 4.Table 4
Eigenvalues,α(g) and α(k).
Eigenvalueλn EigenvectorLgii=1,2,3,4,5,6,7 corresponding to λg α(g) α ( k ) L g 1 L g 2 L g 3 L g 4 L g 5 L g 6 L g 7 4.192 -0.440 -0.387 -0.461 -0.470 0.274 0.318 0.218 0.599 0.599  1.567 -0.323 -0.235 -0.188 -0.191 -0.587 -0.481 -0.439 0.224 0.823  0.813 0.066 -0.460 0.210 0.101 -0.278 -0.291 0.754 0.116 0.939  0.286 0.194 -0.508 0.217 0.137 -0.351 0.662 -0.279 0.041 0.980  0.133 0.206 -0.560 0.163 0.049 0.602 -0.376 -0.335 0.019 0.999  0.009 -0.785 0.005 0.467 0.377 0.135 0.061 -0.035 0.001 1  1.68E -16 -0.021 -0.114 -0.646 0.754 0.017 -0.011 0.014 0 1As can be seen from Table4, the eigenvalues of the first two principal components are greater than 1 and the cumulative contribution rate is greater than 80%. This means that the first two principal components contain most of the information of the original matrix and can be used to evaluate the suitability of the matching areas.According to formula (18) and Table 4, it can be found that the absolute value of standard deviation, x-roughness, y-roughness, and roughness has greater load factors in the first principal component. The absolute values of entropy, kurtosis coefficient, and skewness coefficient have greater load factors in the second principal component.


### 4.2. Index Weights Determination by AHP

According to the principal components and key parameters determined by PCA, a hierarchical structure is established as shown in Table5.Table 5
Hierarchical structure.
A(Objective level) B(Criterion level) C(Index level) Matching area suitability evaluation B1 The first principal C11 standard deviation C12 x-roughness C13 y-roughness C14 roughness B2 The second principal C21 entropy C22 kurtosis coefficient C23 skewness coefficientThe judgment matrixAA is established according to the eigenvalue corresponding to each principal component and rules in Table 1. The larger the eigenvalue, the more important the corresponding principal component. The judgment matrices AB1, AB2 are established according to the feature parameters’ importance in every principal component and the rules in Table 1. AA, AB1, AB2 are shown in Tables 6, 7, and 8.Table 6
Judgment matrixAA.
A A B1 B2 B1 1 5  B2 1/5 1 λ m a x = 2, CI=0, CR=0<0.1ω=0.8330.167TTable 7
Judgment matrixAB1.
A B 1 C11 C12 C13 C14 C11 1 3 1 1  C12 1/3 1 1/3 1/3  C13 1 3 1 1  C14 1 3 1 1 λ m a x = 3.999, CI=-2.501×10-4, CR=-2.778×10-4<0.1ω=0.30.10.30.3TTable 8
Judgment matrixAB2.
A B 2 C21 C22 C23 C21 1 3 3  C22 1/3 1 1  C23 1/3 1 1 λ m a x = 2.999, CI=-3.334×10-4, CR=-5.748×10-4<0.1ω=0.60.20.2TFrom Tables6, 7, and 8, the comprehensive weight can be obtained, as shown in Table 9.Table 9
Comprehensive weight.
Level B1 B2 Comprehensive weight 0.833 0.167 C11 0.3  0.250  C12 0.1  0.083  C13 0.3  0.250  C14 0.3  0.250  C21  0.6 0.100  C22  0.2 0.033  C23  0.2 0.033


### 4.3. Comprehensive Evaluation Value

According to reference [15] and (16)-(18), the comprehensive evaluation values obtained by PCA can be determined by (33); the results are shown in Table 10.(33)Ti′=∑g=1kαgFii=1,2,…,m.Table 10
Comprehensive evaluation values determined by PCA.
Matching area Comprehensive evaluation value Ranking 1 0.3237 4  2 1.0504 2  3 -0.5860 5  4 -2.2633 7  5 -0.6564 6  6 0.8220 3  7 1.3096 1The comprehensive evaluation value of each matching area determined by PCA-AHP can be obtained by (30); the results are shown in Table 11.Table 11
Comprehensive evaluation values determined by PCA-AHP.
Matching area Comprehensive evaluation value Ranking 1 -0.3798 5  2 -0.9400 7  3 0.3747 2  4 1.3553 1  5 0.3718 3  6 -0.2048 4  7 -0.5771 6It can be seen from Tables10 and 11 that the comprehensive evaluation values of the matching regions obtained by PCA and PCA-AHP are quite different. The validity of the two methods will be verified by experiments.


### 4.4. Consistency Analysis of Comprehensive Evaluation Value and Area Suitability

Generally speaking, the greater the comprehensive evaluation value, the better the area suitability. In order to confirm whether the comprehensive evaluation value can reflect the area suitability, other simulation experiments are carried out in these areas. In the above seven areas, a track with the most stable matching performance under zero noise is simulated. When designing a route, it is necessary to avoid track passing through the highly symmetrical areas. After choosing the track, the navigation experiment is carried out with the parameters described in Table12. Random noises with errors of 1 nT, 2 nT, 3 nT, 5 nT, and 8 nT are added to the measured geomagnetic intensity series.Table 12
Matching experiment parameters.
Parameter Value Initial deviation in latitudinal direction 1500m Initial deviation in longitudinal direction 1500m Initial deviation of INS angle 1° INS angle drift error 0.01°/hIf there is no noise in the measured track, the above areas can achieve good matching accuracy. Therefore, in order to reflect the suitability of each area, geomagnetic noise with errors ranging from 1.0 nT to 8.0 nT is added to the measured data. The matching result is shown in Figure3.Figure 3
Performance analysis of different matching areas.From Figure3, it can be clearly seen that the matching suitability of area 4 is the most stable; even if the noise reaches 8 nT, it still maintains a good matching accuracy. Secondly, area 5 keeps a good matching accuracy before the noise reaches 5 nT. Area 3 and area 1 mismatch at 5 nT, but the matching accuracy of area 5 is slightly higher than that of area 1. When the noise is 3 nT, the accuracies of areas 2, 6, and 7 all begin to deteriorate, but the precision of area 2 decreases rapidly, and the other two areas keep good consistency.By analyzing the variation of matching accuracy with the increase of measured geomagnetic errors, the matching suitability of each area is as follows:(34)4>5>3>6>1>7>2.Comparing (34) with the comprehensive evaluation values in Tables 10 and 11, it can be found that the comprehensive evaluation values obtained by PCA-AHP model are more accurate than the comprehensive evaluation values obtained by PCA and have good consistency with the area matching suitability.In the PCA-AHP model:(1)
The ranking of comprehensive evaluation value and matching performance for areas 3 and 5 are slightly different. The comprehensive evaluation values of areas 3 and 5 are 0.3747 and 0.3718, and the difference is very small. Because the standard deviation and roughness of area 5 are larger than those of area 3, the matching accuracy of area 5 is slightly better than that of area 3. The reason for the inconsistency is that only seven feature parameters are used to describe candidate matching areas, and the suitability is not evaluated comprehensively. If there were more parameters, the suitability of the matching area can be more comprehensively evaluated, and the probability of inconsistency would be lower. Therefore, the inconsistency in ranking of areas 3 and 5 is acceptable.(2)
The geomagnetic variation standard deviation corresponds best to the geomagnetic matching performance. The standard deviation of area 4 is 46.993 nT, which is much higher than 9.886 nT in area 2. The experimental results also show that the matching stability of area 4 is much higher than that of area 2. The geomagnetic standard deviation of the remaining areas is also in good agreement with the results described in (34). It shows that the standard deviation is extremely related to the area matching suitability.(3)
The geomagnetic roughness is also the largest in area 4, reaching 4.76 nT, which is higher than that in the remaining areas. Area 2 has the smallest roughness of 1.772 nT. The remaining geomagnetic roughness changes are also in good agreement with (34). It shows that the geomagnetic roughness is related to the area matching suitability.(4)
In the experiment water area, the change of geomagnetic entropy is maintained between 7.053 bits and 7.082 bits, which is not obvious and cannot reflect the change of matching suitability very well. Nevertheless, geomagnetic information entropy reflects the richness of geomagnetic features and is an indispensable feature parameter for judging geomagnetic matching suitability; when the regional area difference is large, it is still recommended to be adopted.(5)
The variation of kurtosis coefficient and skewness coefficient is irregular, which cannot directly reflect the change of matching suitability.In addition, among the seven matching areas, some areas have rich feature information, some feature information is not rich enough, and some areas have similar feature information, with little difference. From the consistency analysis, PCA-AHP can accurately give the evaluation value. Therefore, it is believed that this method can work elsewhere when the conditions are different.The PCA has a large difference between the comprehensive evaluation values and the test results. There may be two reasons as follows.(1)
The abandoned principal component with small contribution rate may contain important information.(2)
The method of calculating the comprehensive evaluation value or determining the weight is inappropriate.


## 4.1. Principal Component Analysis

According to the (10)-(15), the normalized evaluation matrix Z and the correlation matrix R can be obtained.(31)Z=-0.500-0.482-0.182-0.258-0.802-1.4310.685-1.162-0.278-1.067-0.945-0.8020.065-1.3760.4361.2110.0780.2501.069-0.034-0.6831.6831.5991.7921.837-0.802-0.526-0.2410.706-0.4610.7790.638-0.802-0.599-0.623-0.263-0.482-0.596-0.6191.0691.5091.452-0.900-1.106-0.804-0.9031.0691.0170.786.(32)R=1.0000.7630.9690.975-0.227-0.332-0.1630.7631.0000.6940.771-0.118-0.298-0.4070.9690.6941.0000.993-0.412-0.490-0.1870.9750.7710.9931.000-0.397-0.484-0.249-0.227-0.118-0.412-0.3971.0000.7770.485-0.332-0.298-0.490-0.4840.7771.0000.407-0.163-0.407-0.187-0.2490.4850.4071.000.Solve the eigenvaluesλn and eigenvector Lgi of the correlation matrixR and calculate the contribution rate α(g) and cumulative contribution rate α(k) by formula (17) and (18), as listed in Table 4.Table 4
Eigenvalues,α(g) and α(k).
Eigenvalueλn EigenvectorLgii=1,2,3,4,5,6,7 corresponding to λg α(g) α ( k ) L g 1 L g 2 L g 3 L g 4 L g 5 L g 6 L g 7 4.192 -0.440 -0.387 -0.461 -0.470 0.274 0.318 0.218 0.599 0.599  1.567 -0.323 -0.235 -0.188 -0.191 -0.587 -0.481 -0.439 0.224 0.823  0.813 0.066 -0.460 0.210 0.101 -0.278 -0.291 0.754 0.116 0.939  0.286 0.194 -0.508 0.217 0.137 -0.351 0.662 -0.279 0.041 0.980  0.133 0.206 -0.560 0.163 0.049 0.602 -0.376 -0.335 0.019 0.999  0.009 -0.785 0.005 0.467 0.377 0.135 0.061 -0.035 0.001 1  1.68E -16 -0.021 -0.114 -0.646 0.754 0.017 -0.011 0.014 0 1As can be seen from Table4, the eigenvalues of the first two principal components are greater than 1 and the cumulative contribution rate is greater than 80%. This means that the first two principal components contain most of the information of the original matrix and can be used to evaluate the suitability of the matching areas.According to formula (18) and Table 4, it can be found that the absolute value of standard deviation, x-roughness, y-roughness, and roughness has greater load factors in the first principal component. The absolute values of entropy, kurtosis coefficient, and skewness coefficient have greater load factors in the second principal component.


## 4.2. Index Weights Determination by AHP

According to the principal components and key parameters determined by PCA, a hierarchical structure is established as shown in Table5.Table 5
Hierarchical structure.
A(Objective level) B(Criterion level) C(Index level) Matching area suitability evaluation B1 The first principal C11 standard deviation C12 x-roughness C13 y-roughness C14 roughness B2 The second principal C21 entropy C22 kurtosis coefficient C23 skewness coefficientThe judgment matrixAA is established according to the eigenvalue corresponding to each principal component and rules in Table 1. The larger the eigenvalue, the more important the corresponding principal component. The judgment matrices AB1, AB2 are established according to the feature parameters’ importance in every principal component and the rules in Table 1. AA, AB1, AB2 are shown in Tables 6, 7, and 8.Table 6
Judgment matrixAA.
A A B1 B2 B1 1 5  B2 1/5 1 λ m a x = 2, CI=0, CR=0<0.1ω=0.8330.167TTable 7
Judgment matrixAB1.
A B 1 C11 C12 C13 C14 C11 1 3 1 1  C12 1/3 1 1/3 1/3  C13 1 3 1 1  C14 1 3 1 1 λ m a x = 3.999, CI=-2.501×10-4, CR=-2.778×10-4<0.1ω=0.30.10.30.3TTable 8
Judgment matrixAB2.
A B 2 C21 C22 C23 C21 1 3 3  C22 1/3 1 1  C23 1/3 1 1 λ m a x = 2.999, CI=-3.334×10-4, CR=-5.748×10-4<0.1ω=0.60.20.2TFrom Tables6, 7, and 8, the comprehensive weight can be obtained, as shown in Table 9.Table 9
Comprehensive weight.
Level B1 B2 Comprehensive weight 0.833 0.167 C11 0.3  0.250  C12 0.1  0.083  C13 0.3  0.250  C14 0.3  0.250  C21  0.6 0.100  C22  0.2 0.033  C23  0.2 0.033


## 4.3. Comprehensive Evaluation Value

According to reference [15] and (16)-(18), the comprehensive evaluation values obtained by PCA can be determined by (33); the results are shown in Table 10.(33)Ti′=∑g=1kαgFii=1,2,…,m.Table 10
Comprehensive evaluation values determined by PCA.
Matching area Comprehensive evaluation value Ranking 1 0.3237 4  2 1.0504 2  3 -0.5860 5  4 -2.2633 7  5 -0.6564 6  6 0.8220 3  7 1.3096 1The comprehensive evaluation value of each matching area determined by PCA-AHP can be obtained by (30); the results are shown in Table 11.Table 11
Comprehensive evaluation values determined by PCA-AHP.
Matching area Comprehensive evaluation value Ranking 1 -0.3798 5  2 -0.9400 7  3 0.3747 2  4 1.3553 1  5 0.3718 3  6 -0.2048 4  7 -0.5771 6It can be seen from Tables10 and 11 that the comprehensive evaluation values of the matching regions obtained by PCA and PCA-AHP are quite different. The validity of the two methods will be verified by experiments.


## 4.4. Consistency Analysis of Comprehensive Evaluation Value and Area Suitability

Generally speaking, the greater the comprehensive evaluation value, the better the area suitability. In order to confirm whether the comprehensive evaluation value can reflect the area suitability, other simulation experiments are carried out in these areas. In the above seven areas, a track with the most stable matching performance under zero noise is simulated. When designing a route, it is necessary to avoid track passing through the highly symmetrical areas. After choosing the track, the navigation experiment is carried out with the parameters described in Table12. Random noises with errors of 1 nT, 2 nT, 3 nT, 5 nT, and 8 nT are added to the measured geomagnetic intensity series.Table 12
Matching experiment parameters.
Parameter Value Initial deviation in latitudinal direction 1500m Initial deviation in longitudinal direction 1500m Initial deviation of INS angle 1° INS angle drift error 0.01°/hIf there is no noise in the measured track, the above areas can achieve good matching accuracy. Therefore, in order to reflect the suitability of each area, geomagnetic noise with errors ranging from 1.0 nT to 8.0 nT is added to the measured data. The matching result is shown in Figure3.Figure 3
Performance analysis of different matching areas.From Figure3, it can be clearly seen that the matching suitability of area 4 is the most stable; even if the noise reaches 8 nT, it still maintains a good matching accuracy. Secondly, area 5 keeps a good matching accuracy before the noise reaches 5 nT. Area 3 and area 1 mismatch at 5 nT, but the matching accuracy of area 5 is slightly higher than that of area 1. When the noise is 3 nT, the accuracies of areas 2, 6, and 7 all begin to deteriorate, but the precision of area 2 decreases rapidly, and the other two areas keep good consistency.By analyzing the variation of matching accuracy with the increase of measured geomagnetic errors, the matching suitability of each area is as follows:(34)4>5>3>6>1>7>2.Comparing (34) with the comprehensive evaluation values in Tables 10 and 11, it can be found that the comprehensive evaluation values obtained by PCA-AHP model are more accurate than the comprehensive evaluation values obtained by PCA and have good consistency with the area matching suitability.In the PCA-AHP model:(1)
The ranking of comprehensive evaluation value and matching performance for areas 3 and 5 are slightly different. The comprehensive evaluation values of areas 3 and 5 are 0.3747 and 0.3718, and the difference is very small. Because the standard deviation and roughness of area 5 are larger than those of area 3, the matching accuracy of area 5 is slightly better than that of area 3. The reason for the inconsistency is that only seven feature parameters are used to describe candidate matching areas, and the suitability is not evaluated comprehensively. If there were more parameters, the suitability of the matching area can be more comprehensively evaluated, and the probability of inconsistency would be lower. Therefore, the inconsistency in ranking of areas 3 and 5 is acceptable.(2)
The geomagnetic variation standard deviation corresponds best to the geomagnetic matching performance. The standard deviation of area 4 is 46.993 nT, which is much higher than 9.886 nT in area 2. The experimental results also show that the matching stability of area 4 is much higher than that of area 2. The geomagnetic standard deviation of the remaining areas is also in good agreement with the results described in (34). It shows that the standard deviation is extremely related to the area matching suitability.(3)
The geomagnetic roughness is also the largest in area 4, reaching 4.76 nT, which is higher than that in the remaining areas. Area 2 has the smallest roughness of 1.772 nT. The remaining geomagnetic roughness changes are also in good agreement with (34). It shows that the geomagnetic roughness is related to the area matching suitability.(4)
In the experiment water area, the change of geomagnetic entropy is maintained between 7.053 bits and 7.082 bits, which is not obvious and cannot reflect the change of matching suitability very well. Nevertheless, geomagnetic information entropy reflects the richness of geomagnetic features and is an indispensable feature parameter for judging geomagnetic matching suitability; when the regional area difference is large, it is still recommended to be adopted.(5)
The variation of kurtosis coefficient and skewness coefficient is irregular, which cannot directly reflect the change of matching suitability.In addition, among the seven matching areas, some areas have rich feature information, some feature information is not rich enough, and some areas have similar feature information, with little difference. From the consistency analysis, PCA-AHP can accurately give the evaluation value. Therefore, it is believed that this method can work elsewhere when the conditions are different.The PCA has a large difference between the comprehensive evaluation values and the test results. There may be two reasons as follows.(1)
The abandoned principal component with small contribution rate may contain important information.(2)
The method of calculating the comprehensive evaluation value or determining the weight is inappropriate.


## 5. Conclusions

A single geomagnetic feature cannot fully reflect the performance of geomagnetic matching area. To address this problem, a method of multigeomagnetic feature fusion evaluation based on PCA-AHP model is proposed in this paper. The experiment results show that there is an adequate level of consistency between the comprehensive evaluation value and the area matching suitability. It is effective to evaluate the matching suitability of candidate matching areas by using PCA-AHP model. This method can achieve comprehensive evaluation of candidate area suitability by considering multiple geomagnetic features. The conclusions provide a quantitative basis for the analysis of candidate area suitability. In practical application, this method can help select the best matching area from the candidate matching areas, which can be used for the matching area selection of geomagnetic aided navigation.



---

*Source: 9740812-2019-05-20.xml*