Magnetic Signature Description of Ellipsoid-Shape Vessel Using 3D Multi-Dipole Model Fitted on Cardinal Directions

The article presents a continuation of the research on the 3D multi-dipole model applied to the reproduction of magnetic signatures of ferromagnetic objects. The model structure has been modified to improve its flexibility - model parameters determined by optimization can now be located in the cuboid contour representing the object’s hull. To stiffen the model, the training dataset was expanded to data collected from all four cardinal directions. The robustness of the modified multi-dipole model was verified with various noise levels applied to the synthetic data. A comprehensive numerical verification of the proposed methodology was performed using only data not involved in determining the modified multi-dipole model parameters: the data from intercardinal directions and from different depth were used for cross-validation. An analysis of the influence of initial conditions on the optimization process was carried out. In addition to the gradient optimization method, an evolutionary strategy was also used. Regularization was carried out to search for effective model parameterization. New verification methods were also applied based on the balance of magnetic moments and on the average width of the fit error interval. The results of the performed experiments have shown high robustness of the modified multi-dipole model, even in the face of high noise in the input data. The most significant advantage of the model is its predictive ability, enabling determination of magnetic signatures in any directions and depths with high accuracy.


I. INTRODUCTION
Ferromagnetic hulls of ships disrupt locally natural Earth's magnetic field and create a local magnetic anomaly. This anomaly is defined as the ship's own magnetic field, identified by its complex magnetic signature [1], [2]. That magnetic signature has two main magnetization components: induced and permanent [1]. The first component is related to the reaction of ferromagnetic material placed in the Earth's magnetic field and is dependent on the current geographical position and orientation (course) of the ship [3]. The second type of magnetization depends on ship size, its ''magnetic history'' (production and storage of ship's sheet metal, and ship's The associate editor coordinating the review of this manuscript and approving it for publication was Ladislau Matekovits . building technology), ferromagnetic properties of sheets, or even mechanical strikes and temperature stresses during exploitation [3]. While the first component can be easily calculated, the second component has to be estimated based on measurements -its deterministic calculation is not possible without knowing the magnetic history of the object [1]- [4]. The magnetic signature technology has a practical significance for naval transport, as it allows object detection and classification, as well as performing a safety analysis by predicting ships' own magnetization and analysing its own magnetic risk of being detected by naval mines [3], [4].
The main goal of the research reported in this paper was to develop a methodology (multi-dipole model structure and its parameter identification procedure) that allows predicting and determining magnetic signatures of the ferromagnetic object in any directions and depths in the underwater environment with high accuracy, even in the presence of significant noise in the data. For the purposes of this article, it was assumed that perfect knowledge of ship's position in relation to the measuring devices is available. The influence of the error coming from virtual magnetometers is investigated. Comprehensive analyzes of this type can only be performed by simulation with the synthetic data.

A. PREVIOUS WORKS
The results of previous research by the authors on the multidipole model and its first extensions were presented in two papers [5], [6]. The 3D ellipsoidal shell geometry, used as the case study in those papers, is widely accepted as a close and appropriate representation of a vessel hull [8]- [10]. The model simulation results obtained from virtual measurements done by multi-axis magnetometers working in the virtual measuring range equivalent to the real measurement ranges are used in this paper as raw synthetic data with measurement noise [4], [5], [11]- [14], [20].
In the first paper [6], only the vertical magnetic flux density (MFD) component B Z on N (North, 0 • ) and E (East, 270 • ) directions was considered. In that case, the positions of individual dipoles, of permanent and induced nature, were assumed along three main lines: port (P), keel (K ), and starboard (S), located at the deck height symmetrically to the ship hull. The magnetic moments and positions of all dipoles were calculated by solving the appropriately defined optimization task [14]- [17] with the Levenberg-Marquardt optimization algorithm [6] and the use of synthetic input data gathered along P, K and S lines. While, in the second paper [5] it was assumed that ship's magnetic data were available on four cardinal headings, i.e. directions N (North, 0 • ), W (West, 90 • ), S (South, 180 • ) and E (East, 270 • ). Notice that those directions are defined differently than the well-known geographic cardinal directions -the counterclockwise notation is used [5]. The data related to courses N and E were used for determining multi-dipole model parameters, while the data from the other two courses, S and W , were used to verify the resulting multi-dipole model. Particular point dipoles, of both permanent and induced nature, were located not only along three main lines P, K and S [6], but also inside a rectangle on the horizontal plane XY being the approximation enclosing the cross-section of the assumed ellipsoidal shape of the object [5]. The synthetic data gathered along P, K and S lines, and the gradient-based optimization method with constraints (Trust-Region-Reflective optimization algorithm [21]) were used for determining magnetic moments and positions of all considered dipoles in limited ranges.
In [7], the authors synthesized a model that enabled reproducing the magnetic signature on the basis of real data. The main problem was the measurement error in determining the ship position in relation to the magnetometers. The input data included the information about the values of the magnetic field components Bx, By, Bz and their position in relation to the centre of the modelled object. In this article, it is assumed that precise information about object's position is available. Thus, the only sources of errors are those from magnetometers. Different levels of disturbances can be analysed using the simulation approach.

B. NEW CONTRIBUTIONS IN CURRENT PAPER
Compared to the previous papers by the authors [5], [6], in this paper a number of extensions to the multi-dipole model determination methodology have been applied. Firstly, it was assumed that particular induced and permanent point dipoles can be located at arbitrary points inside the 3D cuboid contour representing the object's hull, and not only on the 2D rectangular surface defining its deck [5], [6]. Secondly, the multi-dipole model parameters were determined from the ship's magnetic data gathered along the P, K and S lines at depth -20 m for all four cardinal directions N , W , S and E. For the purpose of the extended multi-dipole model robustness analysis, the synthetic data were supplemented with higher noise levels (10 nT and 100 nT) than in the past research (1 nT) [5], Additionally, the resulting multi-dipole model was verified based on all intercardinal directions NE (North-East, 315 • ), SE (SouthEast, 225 • ), SW (SouthWest, 135 • ), NW (NorthWest, 45 • ), and at a different depth (-30 m). The new verification methods based on the balance of magnetic moments and based on the average width of the fit error interval were also applied.
As the model used in the article is highly non-linear, the optimization procedures are exposed to the risk of being stuck at a local minimum. The problem whether the optimization procedure has got stuck and did not find a satisfactory result has been analyzed with the use of multiple simulations. A discussion on when the obtained result can be treated as correct and useful is also presented. Simulations were carried out to check how different sets of initial conditions affect the course and results of the optimization process. Extensive studies were performed using the regularization procedures L1 and L2 to check the model parameterization and performance.

C. PAPER STRUCTURE
The paper is organized as follows. Sections II and III present brief descriptions of, respectively, the 3D multi-dipole model structure and the 3D longitudinal ellipsoidal shell virtual object as a synthetic data source. Section IV gives the overall idea of the 3D multi-dipole model structure creation, fitting, and verification. Section V reports synthetic benchmarking aspects and results. In Section VI, regularizations L1 and L2 are performed to assess the model parameterization. Section VII gives the aggregated verification of the 3D multi-dipole model. Finally, Section VII concludes the paper. Three Appendices, A, B and C include, respectively, tables with indicators characterising all scenarios analysed in the paper, regularization results for S7 scenario and figures with reconstructed magnetic signatures and fields for selected scenarios.

II. BRIEF DESCRIPTION OF 3D MULTI-DIPOLE MODEL
In the past publications [5], [6], the authors proposed to use the multi-dipole model for magnetic signature reproduction/prediction. This model consists of k dipoles, which describe the magnetic signature of the object as the sum composing the total vector. It is assumed that an appropriate number of n induced dipoles (I) and m permanent dipoles (P), where k = n + m, represent the complex magnetic signature of the analysed object.
In this paper the authors apply a new approach, which gives the ability to place dipoles in the entire 3D space representing cuboid approximation of the vessel hull - Fig. 1.
The multi-dipole model provides more flexibility in reconstructing the entire magnetic field disturbance map around the ferromagnetic source object located at the origin of the Cartesian coordinate system. The vector of magnetic flux density generated at an arbitrary point (x, y, z) by the i-th dipole located at point (x i , y i , z i ), can be described as follows [18], [19]: where i = 1, . . . , n, n + 1, . . . , n + m, and µ 0 =4π×10 −7 H/m is the space permeability, B i is the vector of magnetic flux density of i-th magnetic dipole in each orthogonal direction (x, y, z), M i is the vector of i-th magnetic dipole moments in each orthogonal direction (x, y, z), and R i = |R i | denotes the distance vector of the analyzed point (x, y, z) from the i-th dipole with coordinates (x i , y i , z i ).
Notice that the independent variables in the presented multi-dipole model (eq. (1-4)), are magnetic moments m x,i , m y,i , m z,i and locations x i , y i and z i . of the dipoles. The magnetic moments of the dipoles can have permanent m xP,i , m yP,i , m zP,i (P) or induced nature m I 1,i , m I 2,i , m I 3,i (I). Comparing the above multi-dipole model (1)- (4) to that proposed by the authors in the previous paper [5], the release of positions of all dipoles in the z-axis direction can be observed as the most distinctive difference. The transformation of position and magnetic moment of each considered dipole (induced and permanent) to follow ship course changes is needed. A detailed description of this transformation may be found in [5].

III. THE VIRTUAL 3D LONGITUDINAL ELLIPSOIDAL SHELL OBJECT AS SYTHETIC DATA SOURCE
Technically, there are two possible ways to acquire the magnetic signature of any ferromagnetic object of arbitrary type: from measurements on a real object, or from numerical modelling and simulations of the magnetic field surrounding the object. For the simulation studies described in the paper, the access to the data from a virtual, simulated measuring range with three three-axis magnetometers [13] in each direction (working in star configuration) was assumed.
The 3D ellipsoidal shell geometry ( Fig. 2) used as the case study in this and previous papers [5], [6], is widely accepted as a close and appropriate representation of a naval vessel hull [8]- [10]. All basic geometric and magnetic parameters of the 3D longitudinal ellipsoidal shell geometry are given in Table 1.
Hence, it is possible to get a complex magnetic signature of the object in the form of a complete map of the magnetic field at each point around it, and in the form of partial magnetic signature related to the sensor location. It is additionally assumed that the size of the measurement area is from -100 m to 100 m, and the measurement resolution is 1 m. The data are acquired along P, K and S lines related to object geometry.

IV. OVERALL IDEA OF 3D MULTI-DIPOLE MODEL STRUCTURE CREATION, FITTING, AND VERIFICATION
The purpose of creating a 3D multi-dipole model is the ability to reproduce the magnetic signature of the modelled ferromagnetic object at any direction and depth. It is assumed that there is limited knowledge about magnetic flux density components obtained from measurements during one or more passes of the modelled object through the measuring range. An equivalent of this approach is the research with synthetic data, and cutting out the paths from the fields generated in simulators based on the FEM method [1], [2]. To bring the synthetic data closer to real data features, and to check model's resistance to imperfections, noise is added to the data.
The model has a specific structure and parameters. The structure of the model is shaped by the designer who determines the number of dipoles and their possible locations. The number of dipoles depends on e.g. the distance from the object, and its shape and magnetic properties. As can be seen from the present and previous articles [5], [6], the distribution of dipoles is an important research problem.
Until now, the principle has been in force that the greater the freedom in dipole distribution and the flexibility of the model, the better the model. The values of model parameters are determined by optimization with constraints resulting from the dimensions of the modelled object. In the first verification phase, it is checked how well the model parameters have been selected based on the data on which the learning took place. However, the real test is the cross-validation phase, in which it is checked how well the model determines magnetic flux density values on the directions and depths that were not used in the model training phase. The general workflow of model construction is presented in Fig. 3.

A. DATA FOR MODEL FITTING
Model parameters are selected by optimization so that the model output corresponds to the reference data with the lowest possible error. This phase of model synthesis is called training, or fitting. A very important factor is how much data is used for the matching phase and what is their diversity or otherwise defining the representativeness of signature description. On the one hand, there is a natural expectation that the more input data, the better. On the other hand, however, it is known from the theory of estimation and learning of artificial neural networks that there is a phenomenon of overfitting. Another disadvantage of using too much data is that that data should be physically obtained, e.g. by conducting a measurement campaign based on readings from real magnetometers. From this point of view, the fewer data required to develop a model, the better. The ideal situation would be to determine how much data and what origin is sufficient to build a model of a certain quality. This important task is probably too difficult to put into a mathematical formula, but on the basis of the performed analyzes it is possible to develop the heuristics similarly describing favorable conditions. Generating data in a synthetic way means that we have unlimited data resources, both in terms of quantity and diversity. However, the real advantage of the presented approach is the ability to develop a model based on limited information. The developed approach aims at the possibility of determining the signature on the basis of limited information, i.e. obtained only from the measurement. It has been assumed that the source data is generated in the form of B x , B y and B z components of the magnetic flux density. The obtained data fields are cut through, thus creating data paths on four cardinal directions N (0 • ), W (90 • ), S (180 • ) and E (270 • ). The data from the paths can also be obtained directly from measurement. It was also assumed that three data paths along P, K, and S lines on each course would be used to fit the 3D multi-dipole model.

B. NOISING DATA
Typical noise of fluxgate magnetometers is 6 pTrms/Hz 0.5 [13]. Due to possible industrial magnetic disturbances and the occurrence of variations in the Earth's magnetic field, a reference magnetometer is used in the measurements [5]. In this way the typical noise of magnetic field measurements is of the order of nT.
In the previous article [5], a very good fit of the reference data to the model data was achieved when the data was noisefree. Adding the noise at the level typical for measuring devices meant that the data generated from the model was already burdened with a significant error. This points to either insufficient data used to fit the model, or an incorrect model structure. In any case, the robustness of the model can be VOLUME 10, 2022  questioned in the sense of the lack of its resistance/strength when used for data with imperfections and/or uncertainties. Based on the experience from the previous article [5], some recommendations were formulated with respect to increasing the amount of data in the model learning phase and diversifying their origin in terms of geographical directions.
Another, even more spectacular change introduced to the present model is enabling the poles to be located inside a 3D cuboid covering the ellipsoid, and not only on the surface of a 2D rectangle describing the virtual ship deck. In addition to the introduced modifications, it was decided that the effect of higher noise was also worth analysing. Theoretically, such noise is not actually present due to the properties of the magnetometers, but the noise virtually superimposed on the synthetic data may allow the model robustness to be assessed in the face of above-standard input data disturbances. Therefore, different noise profiles with Gaussian distribution and levels of 1 nT, 10 nT, and 100 nT were prepared. Fig. 4 shows the input data for the calculation procedure with superimposed noise. It can be seen that at extreme values of magnetic flux density, the signal-to-noise ratio is still high, and the noise values are quite small. But firstly, a large part of the data  has values significantly lower than the extreme values, and secondly, the example from the previous article has shown that relatively little noise may cause a lot of perturbation. The data with 1 nT noise can be considered close to reality (i.e. closer than the completely noise-free data), while the data with 10 nT and 100 nT noise were only used to assess the robustness of the model. In the model learning phase, partial information about magnetic flux density was used to obtain an analogy of measurements from three magnetometers. Three paths: P, K and S, were cut according to the idea presented in Fig. 5. Fig. 6 shows an approximate fragment of the data range to illustrate differences in individual noise variants. Fig. 7 shows the zoomed fragment of Fig. 6, covering the range from 20 m to 40 m along the x-axis.

C. CONSTRAINED OPTIMIZATION-BASED MODEL FITTING
The important part of the whole model synthesis procedure is computing the data-fitting model parameters using optimization routines. The optimization procedure selects parameter values and applies them to the assumed model structure.
Next, the obtained model output values are compared with the reference input data. In general, the optimization procedure tries to find such parameters of the model that make the model output fit to the reference data.

D. MODEL VALIDATION
The most important advantage and functionality of the multidipole model is that after determining its parameters, the new magnetic signature can be calculated for an arbitrary course and depth greater than that for which the parameters were determined.
Different verification procedures are available to confirm that the multi-dipole model is constructed correctly and may operate in various conditions, not only in those well fitted with the reference data. Using another set of data, called the verification data, the model can be tested for different directions and different depths. The final verification test consists in comparing the magnetic planes generated by the model with the reference planes cut into paths (see Appendix C).

V. SYNTHETIC BENCHMARKING
The research reported in this article is closely related to the effects and results obtained previously [5]. In order to ensure the continuity of research and the possibility to compare the results achieved previously and now, the same conditions for calculations and verification were assumed. In particular, the presently used input data came from the ellipsoid model with the same dimensions and magnetic properties. The data paths referred to the same P, K and S lines. The same optimization procedure [5], [18] and constraints were used, the data were compared according to the same fitting (FIT) and cross validation (CV) criteria, and the results are presented in tables with the same data layout (see Appendix A). Wherever possible, the same computation conditions were provided, while the studies in which more data was used to train the model or for verification to ensure comparability were performed with averaged values.

A. NUMERICAL SIMULATION SCENARIOS WITH CONTINUATION
As in the previous article [5], several modifications were made to earlier solutions. In order to present and analyse the impact of individual modifications, they were divided into a number of simulation scenarios, collected in Table 2. The new scenarios with important changes analysed in this paper were given new numbers, being the continuation of those used for previously analysed scenarios [5], [6]. The versions of the same scenario with superimposed noise were additionally marked with the letter N and a number denoting the noise level. For instance, S6 means the scenario with an important change and without noise, while S7N10 means the scenario with another important change in the approach to model building, and with data noise at the level of 10 nT. In the previous article [5], the scenario S3 with the noise level of 1 nT was labelled S4 to achieve numbering uniformity. In this article, it is called S3N1.
The base scenario is S3 [5], to which all modifications have been made. The scenario S3 contains released positions of dipoles along the Y coordinate and the data from the 3D magnetometer, but only from two directions N (0 • ) and E (270 • ). The modification applied in scenario S5 is related to the release of dipole positions along the Z coordinate, while that applied in scenario S6 concerns the use of data from four cardinal directions. Scenario S7 combines together the modifications of scenarios S5 and S6. Scenarios S5N1, S6N1 and S7N1 are equivalents of scenarios S5, S6, and S7, respectively, with added noise of 1 nT. The remaining scenarios refer to scenario S7 with different noise levels: 10 nT in scenario S7N10 and 100 nT in S7N100.
In all scenarios, the validation was performed at the depths of -20 m and -30 m. In the scenarios where the model was trained on data from cardinal directions, the data for verification at the depth of -20 m was taken from intercardinal directions NW (45 • ), SW (135 • ), SE (225 • ), and NE (315 • ). The verification data at the depth of -30 m was used to assess how the model with parameters set at -20 m predicts the value of magnetic flux density at a different depth. However, the data at -30 m was not used for comparisons between the scenarios because the flux density at -30 m and -20 m has different values and entering them into comparisons with data set at -20 m would falsify the results.

B. OPTIMIZATION PROBLEM FORMULATION
Due to the expansion of the research range (adding new directions and new possible dimension of dipole placement), it was necessary to structure the objective function. The task of the optimization procedure is still to minimize the distance between the reference and model data on P, K and S lines in all considered directions. The optimization problem for a chosen model structure with m permanent dipoles and n induced dipoles is defined as follows: l ∈ {x, y, z} , and min i , max i are the vectors of minimal and maximal constraint values for the decision variables related to the i-th considered dipole.
The objective function J G (5) defines the difference in matching the reference and model data in all considered directions, for paths P, K and S (9), over the length of 200 m, with the resolution of one meter for magnetic field components B x , B y , B z (8). Depending on the simulation scenario, the objective function takes a different form. Inside the criterion function, there is the sum of squares of model and source data differences for individual magnetic field components.
The number of decision variables depends on whether the location of dipoles is only imposed on the plane, or it is possible to place dipoles at arbitrary locations in the 3D space. In the former case [5], five parameters (m x , m y , m z , x, y) were necessary to describe one dipole, while in the currently proposed case, six parameters (m x , m y , m z , x, y, z) are needed ( (5) and (7)).
To solve the above optimization task (5)-(10), simultaneously ensuring continuity and possibility to compare results, a nonlinear least-squares (nonlinear data-fitting) algorithm with TRR (Trust-Region-Reflective) optimization [21] was used. This algorithm has the ability to take into account constraints (6), which allows to control the position of dipoles and to impose range restrictions on magnetic moments. The mentioned constraints were used in a uniform manner to all n induced dipoles and m permanent dipoles (n + m dipoles in total) -see Table 3.
The structure of the model is determined, among others, by the number of dipoles. For all the scenarios from Table 2, the model structure included 10 permanent and 10 induced dipoles. This article does not analyse the minimal, nor the optimal number of dipoles sufficient to describe and reconstruct the magnetic signature. Some model parameterization evaluation aspects are presented in Chapter VI where regularization was performed. Still, it is certainly an interesting aspect of the application of the multi-dipole model.
The optimization with conservative box constraints ( Table 3) for all dipole locations and magnetic moments was used to solve the dipole location problem. Limiting the space of occurrence of dipoles to the cuboid contour of ellipsoid seems natural, and the limitation of magnetic moments has been treated liberally, i.e. in the range of + /-10 6 A · m 2 , based on the observation of many typical naval ferromagnetic objects modelled with the use of the multi-dipole model. The unconstrained version of optimization for this problem leads to longer computation time and much larger errors, at the same stop conditions for the optimization procedure.

C. RMSE AND MAE INDICATORS
The results of individual simulation scenarios can be compared qualitatively in the form of graph observations, and quantitatively in the form of RMSE (Root Mean Square Error) and MAE (Mean Absolute Error) fit indicators on paths, as well as by comparing the maximum modelling errors. The root mean square error is given as and the mean absolute error is given as where model i is the vector of N signature values at i-th position coordinate counted by the model, and ref i is the vector of N reference signature values at the same position.

D. EVALUATING THE FINAL VALUE OF THE COST FUNCTION
The main threat in the optimization process is getting stuck in a local minimum and not reaching the global minimum. Hence it is worth discussing when the result can be considered correct and useful. The objective function (5)   field components (Bx, By, Bz). This gives a total of 201 * 3 * 4 * 3 = 7236 residues, which are squared before entering the optimization procedure. Based on this data, the optimization procedure selects the values of 120 model parameters. The question is whether the obtained optimization result is correct and useful as the global minimum or represents only a local minimum. After the optimization process is completed, the result in the form of the resnorm index should be referred to the residual number and subjected to the square root operation. For example, let us assume that after completing the optimization process, the resnorm indicator is 20000 nT 2 for 7236 points, i.e. an average of about 2,76 nT 2 for each point. After calculating the square root, an average of 1.66 nT remains for each point. It is also worth knowing the distribution of this error by plotting together the reference and model waveforms. Note that the module values of the reference waveforms reach even as much as several thousand nT. Also note that an additional error may occur due to a 1 nT noise superimposed on the data. The authors do not put a clear numerical limit that determines the usefulness of the results. The results of resnorm amounting to several tens of thousands with fairly even distribution of errors are considered an absolutely acceptable rating.

E. THE IMPACT OF INITIAL CONDITIONS ON THE OPTIMIZATION PROCESS
The optimization problem stated in this paper is highly nonlinear. It is mainly due to the applied model, which can lead to different results after adopting different initial conditions [22]. To ensure completely the same starting conditions, research was conducted from an arbitrarily selected initial point common to all 10 induced dipoles and another point common to 10 permanent dipoles. To check whether the course of the optimization process considered in the paper depends on the choice of initial conditions, four variants were tested: • S7_IC1D Scenario S7 with random selection of initial conditions for 1 dipole,   • S7_ICAD Scenario S7 with random selection of initial conditions for all dipoles, • S7N1_IC1D Scenario S7N1 with random selection of initial conditions for 1 dipole, • S7N1_ICAD Scenario S7N1 with random selection of initial conditions for all dipoles.
The initial conditions were determined randomly taking into account the constraints from Table 3. For variants S7_ICAD and S7N1_ICAD, the initial conditions were drawn for only one dipole and duplicated to the other dipoles. This operation was performed separately for one induced and one permanent dipole. In variants S7_ICAD and S7N1_ICAD, the values of all dipoles were determined independently of each other. The calculations were repeated a dozen times for each variant.
The results in the form of the course of changes of a given optimization process and the final values of the criterion function are presented in Fig. 8 -11. Although the optimization processes for variants with random initial conditions for one dipole have radically different waveforms, all of them were successful. The final values of the criterion function are comparable within the S7 and S7N1 scenarios, which indicates that the optimization effect is not very sensitive to the character of initial conditions. However, the optimization process was more predictable and orderly for variants with all dipoles having random initial conditions.
The variants with random initial conditions for one dipole required more iterations than those with all randomly selected dipoles to reach the final value. As expected, the final values of the criterion function were significantly lower for noisefree scenarios.
When analyzing the situation in which the same simulation scheme gives different final values, the results should be analyzed in the form of an assessment of the entire population of solutions. It is possible to adopt in this case certain measures such as averages, distributions, the best element, etc.

F. CMA-ES AS ALTERNATIVE, GRADIENT FREE OPTIMIZATION METHOD
All the above presented optimization calculations used the lsqnonlin function from the Matlab package. This optimization function is based on the gradient trust-regionreflective algorithm. Good results of parameter estimation were obtained, as confirmed by very low values of resnorm and cross-validation indicators. However, it may be interesting to know whether this good performance is the real feature of gradient algorithms and whether optimization algorithms based on other methods can be equally effective. The issues which are worth comparing include the accuracy of calculations, calculation time, and possible 'getting stuck' in a local minimum. The CMA-ES (Covariance Matrix Adaptation -Evolution Strategy) evolution strategy was chosen as an alternative optimization method. The CMA-ES [23] is developed for numerical optimization of non-linear or non-convex optimization problems based on stochastic and derivative free approach. Evolution strategies are optimization techniques based on the idea of biological evolution. Instead of considering one solution, populations of potential solutions are analysed with mutation and selection as search operators. The multivariate normal distribution is used to propose new individuals. The covariance matrix represents dependencies between the variables. The Covariance Matrix Adaptation is realized with an approximation of the inverse Hessian matrix in the quasi-Newton method. According to [24], CMA-ES performed best among 31 considered algorithms of Genetic and Evolutionary Computation for difficult objective functions and larger budgets (time, number of function evaluations).  The CMA-ES algorithm is not included in Matlab, the external code is available in [25].
In order to obtain a comprehensive overview of CMA-ES properties for the considered computational problem, a dozen calculations were additionally performed for scenario S7 and a dozen for scenario S7N1. The population size was set to 30. The results are presented in Figures 12 and 13.
The CMA-ES algorithm did not get stuck at a local minimum in any of the runs. The final values of the objective function (obtained as a result of algorithm's stagnation) reveal very good fit. As expected, noticeably better results were achieved for the no-noise scenario. The results of the CMA-ES optimization with regard to the achieved final values are comparable to those obtained with the gradient method. However, the calculation time of a single run (10E5 iterations) was about 4 days, while for comparison, the average calculation time using the gradient TRR method with 10E5 iterations was about 6 hours. The CMA-ES method is therefore useful to get acquainted with the solution values that can be achieved in the optimization process, but it is simply very time consuming for regular calculations with many variants.

VI. REGULARIZATION
The model may have too many information variables (and related parameters), which may result in the so-called overfitting, i.e. a perfect match to the training data but with poor predictive capabilities for other data sets. On the other hand, the model may have too few information variables, which in turn will prevent accurate prediction due to an excessively simplified form of the model in relation to the phenomenon that it describes. Regularization, also referred to as penalized regression, is used to control the trade-off between bias and variance in the model. The term associated with the parameter values is introduced into the objective function, which causes that not only the residua, but also the parameter values themselves are minimized (penalized). Depending on the type of regularization, the parameter values may asymptotically decrease to zero, but never reach it (L2 -Ridge regularization) or decrease exactly to zero (L1 -Lasso regularization). This shrinkage of parameters may facilitate the selection of model variables. The criterion for regularization may be the AIC (Akaike information criterion) [26] or BIC (Bayesian information criterion) [27], but the most common and obvious metric is the value of cross validation. When, as a result of the regularization procedure, the parameter values are very small, it may be considered that their contribution to the model is small and these parameters can be removed from the model for better parameterization. The multi-dipole model presented in the article makes use of two types of dipole parameters: related with dipole position in the xyz space and with the values of the magnetic moments Mx, My, Mz. The concept of minimizing the impact of certain quantities on the model in the form of decreasing the parameter values, implemented as part of the regularization, will be performed only for magnetic moments. Reducing the values of the parameters related to the xyz dipole position will result in the concentration of all dipoles close to each other and near the origin. Such a concentration of dipoles may make it difficult to reproduce the field distribution, as a result of which the position of the dipoles will not be additionally penalized. Regularization of the second type of parameters, i.e. the magnetic moments, allows theoretically the optimizer to weaken the influence of the least significant dipoles and remove them from the model. The basis of the regularization is the definition of quality criteria that will be an extension of the criterion given by (5). The L1-Lasso [28] regularization problem for a chosen model structure with m permanent dipoles and n induced dipoles is defined as follows: The L2-Ridge [29], [30] regularization problem for a chosen model structure with m permanent dipoles and n induced dipoles is defined as follows: The regularizations L1 and L2 were performed for the simulation scenarios S7 and S7N1. Figures 14 -21 show the regularization results for the scenario S7N1, while the results for S7 are presented in Appendix B (Figures 25-32 and  Tables 18-22). An important issue related with the regularization is standardization. It is actually necessary when the model includes variables of various types with significantly different ranges of values. Typically, standardization is done using the standard deviation, followed by regularization on standardized variables. In the case discussed in the article, only the variables related to the magnetic moments remain as penalized parameters after eliminating the variables related to the dipole position. They all have the same allowable range and therefore there was no need to standardize them. In the model used in the article, based on the authors' experience, 10 permanent and 10 induced dipoles were used. Each dipole has Mx, My, Mz components, which gives additional 60 parameters in the objective function due to regularization.    The key issue in the regularization process is selecting the value for the hyperparameter λ which determines the weight of the term related to the parameters taken into account in the criterion function. As already mentioned, the cross-validation value was selected as a criterion to assess the usefulness of regularization. Figure 15 shows the result of the regularization processes L1 and L2 for the scenario 7N1. The value of the cross-validation index as a function of λ is presented over the entire range and zoomed in (0,1) range to enable presentation of the selection of most favourable values of λ. Figures 14 and 15 allow to determine the minimum value   of the L1CV and L2CV index. The L1CV' indicator was selected as still having small value of CV but significantly higher value of λ than L1CV. It means that the model parameters for L1CV' will be more shrunken but CV will remain at the same level. Figures 16 -18 show the waveforms of the parameters Mx, My, Mz for L1-Lasso regularization. The CV curve is superimposed (with CVmin and Cvmin' marked) to show the relationship between parameter values, λ and CV. Figures 19 -21 show the waveforms of parameters Mx, My, Mz for L2-Ridge regularization. The CV curve is superimposed (with CVmin and Cvmin' marked) to show the  relationship between parameter values, λ and CV. Tables 4-6 present the parameter values for L1-Lasso regularization: for λ =1E-8 (λstart) and for CVmin and CVmin'.
The data in Tables 4-6 confirms proper operation of the L1 regularization. With the increase of λ (while maintaining a similar value of the cross-validation coefficient), the values of the model parameters decrease, some of them to zero. In theory, the values close to zero can be removed from the model. However, in the present case, a multi-dipole model is used and whole dipoles are suitable for removal rather than their individual components. Tables 7-8 present the parameter values for L2-Ridge regularization: for λ =1E-8 (λstart) and for L2CVmin.
The data in Tables 7-8 confirm proper operation of the L2 regularization. With the increase of λ (while maintaining a similar value of the cross-validation coefficient), the values of the model parameters decrease, at the same time aligning the parameters rather than resetting them.
The summary of the regularization results for scenarios S7N1 and S7: -regularizations L1 and L2 significantly reduce the parameter values, which can be seen by comparing the parameter values for λstart and λ for CVmin and CVmin', -L1-Lasso regularization ensures the minimum CV for much larger lambda values than L2-Ridge, -L1 regularization forces the decrease of some parameters to zero, -L2 regularization tends to evenly distribute values between parameters, -assuming the criterion that a parameter can be removed from the model when all its moment components are less than 1000, on the basis of the L1 regularization presented in Table 6, dipoles P4 and P7 can be excluded from the model, and when slightly extending the criterium, also dipoles P5 and P9. It is noteworthy that such an operation was possible only for λ in which the CV' was achieved, -despite different versions of the L1 vs L2 regularization and different CV values achieved (which is understandable as the regularization objective function differed), very similar values of the sum of the moment modules were achieved, -the results and conclusions are similar for the scenarios S7N1 and S7 -a flat CV characteristic without clear minimum indicates that the regularization did not significantly improve the quality of the model, it could only offer elimination of individual dipoles, -parameterization of the model initially proposed by the authors turned out to be basically correct.
It should be noted that the data used in this article is synthetic and has been simulated. In the noise scenarios, it was white regular noise imposed on the reference data. Therefore, the conditions for conducting the experiments were favourable. Regularization is presented in this article as a useful tool to assist in determining model parameterization. However, the real need to use the regularization and the real profit from its use may be more clearly observed for objects for which real measurements are available.

VII. VERIFICATION
The verification process consists of three stages: balance of the sum of magnetic moments, analysis of fitting on paths in the form of charts and RMSE or MAE indicators, and analysis of matching of magnetic fields on basic directions for which input data was available. The most important stage of verification consists in using the multi-dipole model to determine magnetic signatures of the object for directions and depths other than those for which the object parameters were determined. The comparison with data from real measurements or with synthetic data from simulation environments is a form of assessing the applicability of the method. This ultimately proves the effectiveness of the approach proposed in the article and the usability of the multi-dipole model.

A. BALANCE OF THE SUM OF MAGNETIC MOMENTS
The balance of the sum of magnetic moments is the simple quantity test proposed by the authors to verify the correctness of calculations performed with the developed multi-dipole model. The test compares the resultant magnetic moments of all dipoles of the multi-dipole model (Table 9) with the magnetic moments of the 3D ellipsoidal model being the source of the synthetic magnetic data (Figure 2 and Table 1). Those resultant magnetic moments of the multi-dipole model are determined by solving the optimization task (5)-(10) for the simulation scenario S7 (Table 2) and for the number of induced dipoles n = 10 and permanent dipoles m = 10 (k = n + m = 20).
The vector of permanent magnetic moment of the 3D ellipsoid model has three components with the following values: m xP = 80 424 A · m 2 , m yP = −4 398 A · m 2 , m zP = −31 415 A · m 2 (see Table 1). On the other hand, the sum of permanent magnetic moments of all dipoles determined for the multi-dipole model (Table 9) gives for each orthogonal direction (x, y, z) the following values: 81 854 A · m 2 , −4 458 A · m 2 , and 103 877 A · m 2 . Hence, the relative errors between the permanent moments of both models: m xP along the x-axis and m yP along the y-axis, are at the levels of 1.8% and 1.3%, respectively. Whereas, the induced magnetic moments of the 3D ellipsoid model are as follows (Table 1)  The equivalent induced magnetic moments of the 3D ellipsoid model regarding the transformation described in detail in [5], which are related to the course changes, are determined as follows: m I1 = 44 810 A · m 2 , m I2 = 365 634 A · m 2 , and m I3 = −123 115 A · m 2 . While, for the multi-dipole model ( Table 9) the sum of induced magnetic moments of all determined dipoles gives for all induced components [5] the following values: 44 802 A · m 2 , 365 620 A · m 2 , and − 257 745 A · m 2 . Consequently, the relative errors between the induced components m I1 and m I2 of both models are at the levels of 0.02% and 0.004%, respectively. The sum of permanent and induced vertical component values (along the z-axis) of all dipoles determined for the multi-dipole model magnetic moment is equal to 103 877 A · m 2 − 257 747 A · m 2 = −153 870 A · m 2 , while the corresponding sum for the 3D ellipsoid model is equal to −31 415 A · m 2 −123 115 A · m 2 = −154 532 A · m 2 . Thereby, the relative error between the vertical components of the magnetic moment (the sum of induced and permanent moments) of both models is at the level of 0.4%.
The abovementioned relative errors between basic magnetic components of both models are at the minimal level, which quantitatively confirms the correctness of the proposed multi-dipole model structure and the methodology introduced to determine its parameters. Additionally, the authors present qualitative confirmation of this fact in the form of graphical visualization (Figures 22 and 23) of magnetic moments and locations of the designed dipoles (Table 9) against the background of the 3D ellipsoidal shape (Table 1, Figure 2). The sizes of the dipoles have been referred to the greatest value of the magnetic moment component from the set of permanent and induced dipoles, and multiplied by 5 for the clarity of their graphical presentation: where: m dI ,l,i is the i-th dimensionless equivalent component (l = {1, 2, 3}) of the induced (I) magnetic moments (i = 1, . . . , n) [5], m dP,l,i is the i-th dimensionless component (l = {x, y, z}) of the permanent (P) magnetic moments (i = 1, . . . , m) [5]. It can be seen in Figure 22 and Figure 23 that all induced dipoles (marked in blue) are located along the main axis  (x -axis) of the object, and the largest dipole is located at the end of the 3D ellipsoidal shape. Moreover, two of the determined permanent dipoles are located outside the 3D ellipsoidal shape but inside its outer cubic approximation.

B. PRESENTATION OF RESULTS
The direct result of using the optimization procedure mentioned in the previous section are the multi-dipole model parameters, which may be presented in the form of an appropriate list -e.g. Table 9. They can be applied to similar considerations as those made in the previous section.
However, it is still difficult to interpret these results with regard to the ferromagnetic object magnetic signature reproduction task. Therefore, this chapter presents only the aggregated results of all numerical experiments performed -see   reproduction of magnetic signatures by the proposed model for each analysed case described in detail in Table 2.
In scenario S5 (Table 2), the calculations (fitting procedure) were conducted on directions N (0 • ), E (270 • ), while the results are presented on directions W (90 • ), S (180 • ). These results can be treated as an assessment of the possibility to reproduce magnetic signatures in directions for which no calculations were carried out (cross validation). The quality of optimization (fitting) was checked on paths, while the quality of model validation was checked on planes. To assess the quality of the model in all cardinal directions (scenarios S6, S6N1, S7, 7N1, 7N10, 7N100 described in Table 2   (Table 2, Figure 24). Where the FIT task is considered, noisy input data is compared with a smooth representation of the 3D ellipsoidal model. In contrast, in CV comparisons the smooth reference data without noise is compared with the smooth output data from the model.
Comparing model reproductions with the reference data on which the model was not trained bears the name of cross validation (CV). The aim of the research undertaken in this article was to introduce changes in model structure and input data sources, and to assess the improvement brought by these changes. Cross validation data were used to compile the model quality comparison in each version written in different scenarios.
The quality of multi-dipole models can be directly compared based on the results presented in Tables in Appendix A.  The selected results are also presented in the graphical form (figures in Appendix C). An extra effort has been made to make the comparison more clear. Hence, the overall results are presented in Figure 24.
Using the CV data only from the depth of -20 m, the average width of the intervals between the errors dB x , dB y , and dB z (differences between reference data and resulting multi-dipole model data) were calculated according to the following formula:  Table 2), and d cv2 = {90 • , 180 • } for scenarios S5 and S5N (see Table 2). The average value of  error range, counted on CV values for each scenario, is treated as the comparison indicator -described as the scenario score (15).

C. DISCUSSION OF RESULTS-MODEL QUALITY ASSESSMENT
The discussion of the results was carried out on the data from tables in Appendix A (Table 5-Table 11, for scenarios S5, S5N1, S6, S6N1, S7, S7N1, S7N10, S7N100) and figures in Appendix C (Figures 33-64, fitting and cross validation results for selected scenarios S7 and S7N1). Figure 24 was also used in this discussion. The model modification consisting in the possibility of positioning dipoles in the 3D space brought ambiguous results: some improvement was observed in scenarios S5N1  vs S3N1 and S7N1 vs S6N1, but also some deterioration in scenarios S4 vs S3 [5] and S7 vs S6. A huge improvement in model quality was brought by adding the input data from all cardinal directions, which can be observed when comparing scenarios S6 with S3, and S6N1 with S3N1 [5]. The combined use of the 3D multi-dipole model structure and the data from four cardinal directions used in scenario S7 caused a very low error level, accompanied by very high noise immunity. Applying the 10 times higher noise to the input data in scenario S7N10 caused a much lower cross-validation level than that in scenario S4 (S3N1). Only the 100 times higher noise level in scenario S7N100 resulted in exceeding the cross-validation score obtained in scenario S4 (S3N1).
Worth emphasizing are the values of the cross-validation index close to zero in scenarios S6 and S7, as well as in S6N1 and S7N1. This means that for uninterrupted data and for real-level noise data, the cross-validation values were so low that there was no room for significant improvement.

VIII. CONCLUSION
Since enriching the availability of measurement data and making the model structure more flexible resulted in a   significant improvement in the quality of the multi-dipole model (results reported in the previous paper), it was decided to develop further the model, as described in this paper.   Thus, the data for the model fitting phase came from four cardinal geographic directions, and the model structure allowed the dipoles to be located in the three-dimensional space of the cuboidal contour approximating the 3D ellipsoid object model. It was also decided to extend model verification to eight paths and two different depths. The noise tests of the previous version ended with the noise level of real magnetometers, while this article considers 10 and 100 times greater noise to test model robustness.
In the present case, making the model structure more flexible did not bring clear improvement in model quality.                           Depending on the simulated scenario, either improvement or deterioration of the results was observed. The cause may be   that the input test data comes from the plane, i.e. the 2D source, and therefore data diversity for the structure of the 3D multi-dipole model is insufficient. Further research on this issue is planned. The second modification, consisting in training the model with a significantly larger portion of data, brought a great improvement in all analysed scenarios, both with and without noise. Applying the data from four cardinal directions caused that the modelling error was practically close to zero. A combined application of the above modifications resulted in slight increase in the error value, but at a very low level, which proves a significant improvement in modelling compared to the results presented in the previous article. Finally, the strength of the model to measurement data noise was checked.  The model has adequate stiffness and solidity, which can be seen in that the noise 10 times greater than that used in the previous article causes less error, and only the 100 times greater noise causes a larger error than previously. The noise level used to test the strength of the model is many times higher than that in real measuring devices. Such high noise levels were adopted only to analyse the robustness of the model. After summarizing the use of data from cardinal directions and allowing for spatial distribution of dipoles, the obtained model has enormous resistance to interference and gives very low error values, verified at different depths. Therefore, it can be clearly said that the applied changes brought a significant quality improvement in the model when used with noisefree data. However, the greatest qualitative progress of the model refers to its robustness, verified by learning on abovestandard noisy data.
The work undertaken in the article also improved the state of knowledge on the optimization process and model parameterization. The impact of the initial conditions turned out to be large in the transition phase, but with a sufficiently large number of iterations, all experiments ended with results at an acceptable low level of the criterion value. The alternative optimization method based on evolutionary strategies turned out to be effective in terms of determining the parameters of the model, however, it was many times slower than the gradient method. The regularization used to evaluate the parameterization of the model indicated the possibility of a slight model reduction, but in general it confirmed the correct  form used in the article. The application of regularization can be particularly useful for studies making use of real measurements.

APPENDIX A RESULTS OF EXPERIMENTS-TABLES
Tables 10 trough 17 present synthetic results of fitting and cross-validation (numeric values of considered indicators) on various directions and depths, for scenarios S5, S5N1, S6, S6N1, S7, S7N1, S7N10, and S7N100.
All those scenarios are described in Section V, and their fitting and cross-validation details are given in Table 2. For reader's convenience, the fitting procedure outcomes are marked red, and outcomes of the cross-validation procedure -black.

APPENDIX B RESULTS OF REGULARIZATION FOR S7
See Figures 25-32 and Tables 19-22.