Joint solution method and system for locating radiation sources and constructing spectrum maps based on terrain

Abstract

The present invention belongs to the technical field of electromagnetic spectrum map construction, and specifically relates to a joint solution method and system for radiation source positioning and spectrum map construction combined with terrain. The method comprises the following steps: S1, constructing a rough spectrum map through the signal power values ​​collected in the area, and then determining the radiation source position search domain based on the rough spectrum map; S2, based on the Longley-Rice semi-empirical propagation model, jointly iteratively optimizing the radiation source position and the parameters of the propagation model in the radiation source position search domain; S3, after the iteration converges, forward calculation is performed using the optimized radiation source position and the parameters of the propagation model to obtain the final spectrum map. This method simultaneously meets the two practical needs of radiation source position positioning and spectrum map construction. By utilizing the Longley-Rice propagation model, the present invention takes into account the influence of terrain environmental factors, making the model closer to reality. The present invention takes into account terrain elements in the complete solution process, meeting the requirements for higher-precision spectrum map construction.

Description

Technical Field

[0001] The present invention belongs to the technical field of electromagnetic spectrum map construction, and specifically relates to a joint solution method and system for radiation source positioning and spectrum map construction combined with terrain. Background Art

[0002] There are a large number of radiation source devices in the modern wireless electromagnetic environment, which will form a time-space related energy difference distribution at their working frequency points. Spectrum maps can describe spectrum data from multiple dimensions such as geographic location, frequency, time, and energy, and can integrate and present multi-dimensional electromagnetic spectrum space data. Accurate electromagnetic spectrum maps can provide support for many applications such as spectrum resource utilization, frequency behavior analysis and prediction, and user location spatial positioning. In the field of people's livelihood, with the development of cities and the increase in frequency density, electromagnetic pollution problems are becoming increasingly prominent, and the construction of spectrum maps is beneficial to electromagnetic pollution monitoring and prevention. For example, the existence of "black radio" and "fake base stations" seriously affects normal communications, and spectrum maps can help law enforcement personnel analyze the laws of electromagnetic interference behavior, thereby investigating and punishing illegal frequency equipment and maintaining frequency order. In the field of military applications, mastering electromagnetic control and seizing spectrum advantages is the key to winning combat. By constructing spectrum maps, regional electromagnetic energy distribution can be obtained, and combat teams can use spectrum maps to quickly distinguish between enemies and friends and explore the enemy, thereby greatly improving the team's ability to locate enemy radio equipment, flexibly ensure their own communications, and assess battlefield damage during operations.

[0003] The construction of spectrum maps is usually completed based on the observations of sparsely deployed sensors in the spatial domain. According to whether the propagation prior information is used, the construction methods can be divided into two categories: parametric methods and non-parametric methods. For parametric methods, existing research focuses on solving them in the form of statistical empirical models, which cannot take into account the specific terrain factors of the propagation environment. Summary of the invention

[0004] In view of the above-mentioned deficiencies, the object of the present invention is to provide a method, device and storage medium for jointly solving radiation source positioning and spectrum map construction by taking terrain elements into consideration, realizing joint solution of radiation source position and spectrum map and obtaining higher precision.

[0005] In order to achieve the above object, the present invention provides the following technical solutions:

[0006] In a first aspect, the present invention provides a method for jointly solving radiation source positioning and spectrum map construction, comprising the following steps:

[0007] S1. construct a rough spectrum map through the signal power values ​​collected in the area, and then determine the radiation source location search domain based on the rough spectrum map;

[0008] S2, based on the Longley-Rice semi-empirical propagation model, the radiation source position and the parameters of the propagation model are jointly iteratively optimized in the radiation source position search domain;

[0009] S3. After the iteration converges, the optimized radiation source position and propagation model parameters are used to perform forward calculation to obtain the final spectrum map.

[0010] As an optional embodiment of the first aspect, in S1, the signal power value is obtained by measuring sensor nodes sparsely deployed in an area.

[0011] As an optional embodiment of the first aspect, in S1, the signal power values ​​collected in the area are used to construct a rough spectrum map through a non-parametric method.

[0012] Furthermore, the non-parametric method is an inverse distance weighted method, or a spatial interpolation method, or a kernel learning method, or a graph signal processing method.

[0013] As an optional embodiment of the first aspect, in S1, the method for determining the radiation source position search domain is: sorting the signal power values ​​in the approximate spectrum map from large to small, and then selecting the spatial position of the top-ranked signal power value as the radiation source position search domain, and recording the potential spatial position set as pos.

[0014] Furthermore, in S2, within the spatial position set of the radiation source position search domain, based on a heuristic algorithm, the radiation source position and the propagation model parameters are jointly searched and corrected by taking the minimum difference between the signal power calculated by the propagation model and the actual measured position as the objective function.

[0015] Furthermore, the heuristic algorithm is a particle swarm optimization algorithm.

[0016] Furthermore, the error of the propagation model includes a multiplicative parameter error and an additive parameter error, which are corrected using a coefficient k and a constant Const, respectively.

[0017] Const=P TX +P TX_Gain +P RX_Grain -P RX_Loss ;

[0018] Among them, P TX is the transmitting antenna power, P TX_Gain is the transmitting antenna gain, P RX_Gain is the receiving antenna gain, P RX_Loss is the receiving antenna loss.

[0019] Furthermore, the general spectrum map construction method is shown in the following formula:

[0020]

[0021] Where N represents the number of sensor nodes deployed in a spatial area of ​​size k×l, (x i ,y i ) represents the spatial position corresponding to the sensor node, P i ,i∈1,2,...,N represents the signal power value obtained by monitoring, F(·) represents the non-parametric method, and REMcoarse represents the rough spectrum map;

[0022] When there are N sensor nodes in the space, the corresponding received power value is recorded as [P1, P2, ..., P N ], in the estimated radiation source position (X E ,Y E ) is calculated by the Longley-Rice propagation model as [PL1(X E ,Y E ),PL2(X E ,Y E ),...,PL N (X E ,Y E )], the measurement noise is recorded as ε, then the calibration equation is as follows:

[0023]

[0024] set up:

[0025]

[0026] Then the objective function is:

[0027]

[0028] τ=Λ + Q=(Λ T Λ) -1 Λ T Q;

[0029] After the heuristic algorithm converges, the optimal solution (X E ,Y E ); The radiation source position (X E ,Y E ), optimize the model parameters τ obtained by the calibration model, use the Longley-Rice model combined with the regional digital elevation map to calculate the regional spectrum map REMcal; use the Longley-Rice model to calculate the difference between the theoretical received signal strength value and the actual received signal strength value of each sensor node, recorded as [Shad1, Shad2, ..., Shad N]; using the non-parametric method F(·), with (x i ,y i ,Shad i ) is used to calculate the correction matrix REMcor of the constructed area for the measured values; the spectrum map REM output by the final algorithm is REMcal+REMcor.

[0030] In a second aspect, the present invention provides a computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein when the processor executes the computer program, the processor implements the first aspect or any possible method for jointly solving the radiation source positioning and spectrum map construction of the first aspect.

[0031] In a third aspect, the present invention provides a computer storage medium storing a computer program, characterized in that when the computer program is executed by a processor, it implements the first aspect or any possible joint solution method for radiation source positioning and spectrum map construction in the first aspect.

[0032] The beneficial effects of the present invention are:

[0033] (1) The present invention combines the sparsely collected spectrum signal strength values ​​with the digital elevation terrain, and uses the power values ​​measured by the sparsely deployed sensor nodes in the spatial domain to realize spectrum map construction and radiation source location solution. The method simultaneously meets the two practical needs of radiation source location and spectrum map construction.

[0034] (2) The present invention takes into account the influence of terrain and environmental factors by utilizing the Longley-Rice propagation model, making the model closer to reality. It is proposed to summarize the correction terms of the Longley-Rice calculation model into additive parameter factors and multiplicative parameter factors so that they can be iteratively calculated. By utilizing a heuristic algorithm, the joint optimization iteration of spectrum map construction and radiation source positioning is realized, and finally the radiation source positioning and spectrum map construction are jointly solved, and the solution results are combined to construct the final airspace spectrum map. The present invention takes terrain elements into account in the complete solution process, meeting the requirements for higher-precision spectrum map construction. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 It is a flow chart of the joint solution method of radiation source location and spectrum map construction;

[0036] Figure 2 This is the calculation flow chart for the Longley-Rice model;

[0037] Figure 3 It is a schematic diagram of the experimental verification site and result presentation;

[0038] Figure 4This is a simplified process example diagram of the present invention in an open air communication scenario;

[0039] Figure 5 is a flow chart of a method for jointly solving radiation source positioning and spectrum map construction in one embodiment of the present invention;

[0040] Figure 6 It is a spectrum monitoring unit node design structure diagram. DETAILED DESCRIPTION

[0041] Embodiment 1

[0042] like Figure 1 As shown, a joint solution method for locating a radiation source and constructing a spectrum map includes the following steps:

[0043] S1. construct a rough spectrum map through the signal power values ​​collected in the area, and then determine the radiation source location search domain based on the rough spectrum map;

[0044] S2, based on the Longley-Rice semi-empirical propagation model, the radiation source position and the parameters of the propagation model are jointly iteratively optimized in the radiation source position search domain;

[0045] S3. After the iteration converges, the optimized radiation source position and propagation model parameters are used to perform forward calculation to obtain the final spectrum map.

[0046] Compared with the existing method, the solution of the present invention takes terrain factors into consideration by using a semi-empirical model. In unfamiliar areas, the radiation source location and spectrum map construction can be jointly solved to obtain a spectrum map and radiation source location with high accuracy, which has strong application value.

[0047] Embodiment 2

[0048] In this embodiment, for step S1, the signal power value is measured by sensor nodes sparsely deployed in the area, and can also be collected by other devices that can collect signal power values.

[0049] Embodiment 3

[0050] In this embodiment, for step S1, the signal power values ​​collected in the area are used to construct a rough spectrum map by a non-parametric method. The non-parametric method can be any one of an inverse distance weighted method, a spatial interpolation method, a kernel learning method, and a graph signal processing method.

[0051] Embodiment 4

[0052] In this embodiment, for step S1, the method for determining the radiation source position search domain is: sorting the signal power values ​​in the general spectrum map from large to small, and then selecting the spatial position of the top-ranked signal power value as the radiation source position search domain.

[0053] Assume that N sensor nodes are deployed in a space area of ​​size k×l, and the corresponding spatial position is (x i ,y i ), the observed signal power monitoring value corresponds to P i ,i∈1,2,...,N. First, we can use (x i ,y i ,P i ), the spatial spectrum map is completed by the non-parametric method F(·), so as to obtain a rough spectrum map, which is represented by the symbol REMcoarse. This process can be expressed as follows:

[0054]

[0055] According to the principle that the higher the signal power, the more likely it is to be the location of a potential radiation source, all signal power values ​​in REMcoarse are sorted from large to small, and then the spatial positions where the top c% intensity values ​​with relatively high power are located are selected as the radiation source position search domain. The value of c can be determined independently according to actual conditions. After the radiation source position search domain is determined, the potential spatial position set is recorded as pos.

[0056] For step S2, the Longley-Rice model is used as a propagation model to characterize the spatial signal propagation. Within the spatial position set pos of the radiation source position search domain, based on a heuristic algorithm, the minimum difference between the signal power calculated by the propagation model and the actual measured position is used as the objective function to perform a joint search and correction of the radiation source position and the propagation model parameters.

[0057] The Longley-Rice model is a semi-empirical model with an application frequency range of 20MHz to 40GHz. It uses the path geometry of the terrain and the diffraction and scattering characteristics of the troposphere to predict the median transmission loss. The total path loss calculated by this model can be expressed as A free (d,H)+A ref (d,H), where A free (d,H) represents the free space loss from transmitting to receiving, A ref (d,H) represents a semi-empirical calculation method combining digital elevation along the path from the transmitting source to the receiving device. The overall calculation procedure of Longley-Rice is as follows Figure 2As shown in the figure, its input variables include the terrain parameters between the transmitting and receiving paths, the heights of the transmitting and receiving antennas h1 and h2, the operating frequency of the transmitter fMHz, the radio wave polarization mode ipol, the monthly average atmospheric refractive index en0 in the radar wave propagation area, the ground relative dielectric constant ε and the ground conductivity σ. The terrain parameters are obtained from the digital elevation map of a given resolution, and the digital elevation matrix of the calculation area is represented as DEM. In actual use, the heights of the transmitting and receiving antennas h1 and h2, the operating frequency of the transmitter fMHz, the radio wave polarization mode ipol, the monthly average atmospheric refractive index en0 in the radar wave propagation area, the ground relative dielectric constant ε and the ground conductivity σ can be obtained through prior related knowledge such as environmental measurements and historical data. According to the distance from the prediction point to the emission source, the propagation calculation mode is divided into the line-of-sight zone, the diffraction zone and the scattering zone, where the line-of-sight calculation threshold distance is dLs and the scattering calculation threshold distance is dLs. x , and then perform the corresponding fading calculation according to the calculation formula of the corresponding mode.

[0058] The role of the propagation model is to describe the spatial power attenuation trend from a statistical perspective by combining the terrain. However, due to the complexity of the actual propagation environment, there are differences between the propagation model and the actual environment. Therefore, the error of the propagation model includes multiplicative parameter error and additive parameter error, which are corrected by the coefficient k and the constant Const respectively.

[0059] The calculation formula of constant Const is:

[0060] Const=P TX +P TX_Gain +P RX_Grain -P RX_Loss ;

[0061] Among them, P TX is the transmitting antenna power, P TX_Gain is the transmitting antenna gain, P RX_Gain is the receiving antenna gain, P RX_Loss is the receiving antenna loss.

[0062] Assume that there is a single radiation source base station in the space, and there are N sensor nodes in total. The corresponding received power value is recorded as [P1, P2, ..., P N ], in the estimated radiation source position (X E ,Y E ) is calculated by the Longley-Rice propagation model as [PL1(X E ,Y E ),PL2(X E ,Y E ),...,PL N (X E ,Y E)], the measurement noise is recorded as ε, then the calibration equation is as follows:

[0063]

[0064] To express the calculation process more clearly, assume:

[0065]

[0066] Then the objective function is:

[0067]

[0068] τ=Λ + Q=(Λ T Λ) -1 Λ T Q.

[0069] Since the Longley-Rice model combined with terrain is a calculation process rather than an analytical solution calculation formula, the present invention uses a heuristic algorithm to solve the radiation source position within the position search domain pos, with the minimum gap between the model and the observed value as the objective function, and obtains the model calibration variable τ through the calculated radiation source position.

[0070] After the heuristic algorithm converges, the optimal solution (X E ,Y E ). Using the radiation source position obtained by optimization (X E ,Y E ), optimize the model parameters τ obtained by the calibration model, and use the Longley-Rice model combined with the regional digital elevation map to calculate the regional spectrum map REMcal. At the same time, the Longley-Rice model is used to calculate the difference between the theoretical received signal strength value and the actual received signal strength value of each sensor node, recorded as [Shad1, Shad2, ..., Shad N ] Using the nonparametric method F(·), (x i ,y i ,Shad i ) is used to calculate the correction matrix REMcor of the constructed area for the measured values. The spectrum map REM output by the final algorithm is REMcal+REMcor.

[0071] like Figure 5As shown, in step 101, the input of the construction method is the sparse measurement points in space and the digital elevation map DEM of the region; in step 102, the power values ​​observed by the spatial sensor are used to obtain the rough spectrum map REMcoarse through a non-parametric construction method. The power values ​​of the rough spectrum map are sorted from large to small, and the top c% of the high power areas are selected and defined as the potential spatial position set, which is recorded as pos; in step 103, a heuristic algorithm is used to iteratively optimize the radiation source position and the parameters related to the propagation model in the pos set with the minimum difference between the observed and model calculated power values ​​as the objective function, and finally minimize the objective function value; in step 104, after the algorithm converges, the final radiation source position and the propagation model correction parameters are brought into the forward calculation to obtain the final spectrum map REM.

[0072] In a specific experiment, Figure 3 The actual terrain and spectrum map reconstruction results of the experimental verification conducted in the mountainous environment of Yuelu District, Changsha City in December 2022 are shown. In the real terrain environment of the L-band, a single radiation source is set up and 22 L-band spectrum monitoring sensors are deployed. The system settings of the sensors are as follows Figure 6 As shown in the figure, the spectrum monitoring data is aggregated to the data center using the wireless data transmission module of the sensor. The spectrum analyzer is used as a calibration instrument in the region to meet the speed requirements of the quasi-Lee criterion for radio measurement. It is run along multiple paths along the main road at a speed of less than 10m / s for a distance for testing. The IDW non-parametric method and the particle swarm optimization heuristic optimization algorithm can effectively construct the spectrum map of the region.

[0073] First, a commonly used inverse distance weighted (IDW) non-parametric method is used to construct a spatial spectrum map. As shown in the following formula:

[0074]

[0075] Among them, the node power to be predicted is d represents the Euclidean distance between the known sensor node and the predicted node, and α is the inverse distance weighting coefficient, which can usually be taken as 2 to 5, and in this case is taken as 3. By taking the first 30% area of ​​the power value, the radiation source position search domain pos is obtained. At this stage, the algorithm input is (x i ,y i ,P i ), the output is the search domain pos of the radiation source position.

[0076] A particle swarm optimization heuristic method is used to iteratively calculate the radiation source position [X, Y] and the model correction parameter [Const, k]. An example of the particle swarm optimization algorithm is shown in Algorithm 1 (in this example, a single radiation source scenario). The number of particles is denoted by ψ, and the initial set of optimized particles is set to the search domain pos. f(·) represents the objective function in the formula, w(·) represents the inertia weight, c1, c2 represent the acceleration constants, r 1,j ,r 2,j Follows uniform distribution, y a,j (t) represents the position of particle a that minimizes the objective function in the current iteration number t. Represents the best position of the jth dimension as of iteration t. The formulas involved in the update part of the particle swarm optimization method include:

[0077]

[0078] x a,j (t+1)=x a,j (t)+v a,j (t+1),

[0079] In the formula, v a,j (t+1) and v a,j (t) represents the updated value of particle a in the jth parameter dimension at t and t+1 respectively; where x represents the updated value of particle x in the t+1th iteration a,j In this example, the two-dimensional plane position of one radiation source is optimized, so j is 2. In a multi-radiation source scenario, the present invention can be applied by expanding the dimension.

[0080]

[0081] Finally, after the iteration converges, the estimated radiation source position [X, Y] and the correction parameter set [Const, k] are obtained. The regional spectrum map REMcal is calculated using the Longley-Rice model combined with the regional digital elevation map. At the same time, the Longley-Rice model is used to calculate the difference between the theoretical received signal strength value and the actual received signal strength value of each sensor node, recorded as [Shad1, Shad2, ..., Shad N Using the IDW method, (x i ,y i ,Shad i ) is used to calculate the correction matrix REMcor of the constructed area for the measured values. The sum of REMcal and REMcor gives the spectrum map REM output by the final algorithm.

[0082] In order to quantify the performance of the present invention, the RMSE indicator is used to analyze the accuracy of spectrum map construction:

[0083]

[0084] Where L represents the total number of evaluation points, represents the algorithm reconstruction signal strength, P i Represents the actual measured signal power value.

[0085] In this example, the spectrum map reconstructed by this method has an RMSE accuracy of 4.4dB on the test evaluation data using a total of 22 sensor nodes, and the RMSE accuracy reconstructed using the Kriging interpolation method is 6.5dB. In this example, the two-dimensional plane positioning accuracy of the radiation source of this method is 70m. Using the transmitter location estimation method (LIvE) and the received signal differential positioning method (RSSD) as comparison methods, their two-dimensional positioning accuracies are 146m and 150m respectively.

[0086] To demonstrate the algorithm intuitively, Figure 4 The simplified process example of the present invention in an open air communication scenario is shown. The asterisks in the figure represent the deployed sensor nodes. In the schematic diagram, only [X, Y, Const, k] is optimized, and the calculation process does not involve the terrain.

[0087] The non-parametric construction method IDW and the heuristic algorithm particle swarm optimization selected in this example description can be replaced by algorithms of the same type according to requirements such as computational efficiency and reconstruction accuracy.

[0088] The present invention also proposes a computer device, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the computer program, the above-mentioned joint solution method for radiation source positioning and spectrum map construction is implemented.

[0089] The present invention also proposes a computer storage medium, which stores a computer program. When the computer program is executed by a processor, the above-mentioned radiation source positioning and spectrum map construction joint solution method is implemented.

[0090] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention is described in detail with reference to the aforementioned embodiments, those skilled in the art can still modify the technical solutions described in the aforementioned embodiments or replace some of the technical features therein by equivalents. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A method for jointly solving radiation source positioning and spectrum map construction, characterized in that: The following steps are involved: S1. construct a rough spectrum map through the signal power values ​​collected in the area, and then determine the radiation source location search domain based on the rough spectrum map; S2, based on the Longley-Rice semi-empirical propagation model, the radiation source position and the parameters of the propagation model are jointly iteratively optimized in the radiation source position search domain; In S2, within the spatial position set of the radiation source position search domain, based on the heuristic algorithm, the minimum difference between the signal power calculated by the propagation model and the actual measured position is used as the objective function to perform a joint search and correction of the radiation source position and the propagation model parameters; S3. After the iteration converges, the optimized radiation source position and propagation model parameters are used to perform forward calculation to obtain the final spectrum map.

2. The method for jointly solving radiation source positioning and spectrum map construction according to claim 1 is characterized in that: In S1, the signal power value is measured by sensor nodes sparsely deployed in the area.

3. The method for jointly solving radiation source positioning and spectrum map construction according to claim 1 is characterized in that: In S1, the signal power values ​​collected in the area are used to construct a rough spectrum map using a non-parametric method.

4. The method for jointly solving radiation source positioning and spectrum map construction according to claim 3 is characterized in that: The non-parametric method is an inverse distance weighted method, or a spatial interpolation method, or a kernel learning method, or a graph signal processing method.

5. The method for jointly solving radiation source positioning and spectrum map construction according to claim 1, characterized in that: In S1, the method for determining the radiation source position search domain is: sort the signal power values ​​in the rough spectrum map from large to small, and then select the spatial position where the signal power value ranked first is located as the radiation source position search domain, and record the potential spatial position set as POS .

6. The method for jointly solving radiation source positioning and spectrum map construction according to claim 5 is characterized in that: The heuristic algorithm is a particle swarm optimization algorithm.

7. The method for jointly solving radiation source positioning and spectrum map construction according to claim 5, characterized in that: The error of the propagation model includes multiplicative parameter error and additive parameter error, which are respectively calculated using coefficients k , and constants Const Make corrections.

8. The method for jointly solving radiation source positioning and spectrum map construction according to claim 7 is characterized in that: ； in, is the transmitting antenna power, is the transmit antenna gain, is the receiving antenna gain, is the receiving antenna loss.

9. The method for jointly solving radiation source positioning and spectrum map construction according to claim 8, characterized in that: The general spectrum map construction method is shown in the following formula: ； in, N The representative size is The number of sensor nodes deployed in the spatial area, represents the spatial position corresponding to the sensor node, Represents the signal power value obtained by monitoring, represents a nonparametric method, REMcoarse represents a rough spectrum map; When there are N sensor nodes in the space, the corresponding received power value is recorded as , in estimating the radiation source position The result calculated by the Longley-Rice propagation model is , the measurement noise is recorded as , then the calibration equation is as follows: ； set up: ， ， ， Then the objective function is: ， ； After the heuristic algorithm iteratively converges, the solution that minimizes the objective function is obtained. ; The radiation source position obtained by optimization , optimize the model parameters obtained by the calibration model , using the Longley-Rice model combined with regional digital elevation maps to calculate regional spectrum maps REMcal ; The difference between the theoretical received signal strength value and the actual received signal strength value of each sensor node is calculated using the Longley-Rice model, which is recorded as ; Using nonparametric methods ,by Calculate the correction matrix for the construction area for the measurements REMcor ; The spectrum map output by the final algorithm REM for REMcal + REMcor .

10. A computer device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that: When the processor executes the computer program, the method for jointly solving the radiation source positioning and spectrum map construction as described in any one of claims 1 to 9 is implemented.

11. A computer storage medium storing a computer program, characterized in that: When the computer program is executed by a processor, the method for jointly solving the radiation source positioning and spectrum map construction as described in any one of claims 1 to 9 is implemented.

Images