A Node Localization Algorithm Based on Multi-Granularity Regional Division and the Lagrange Multiplier Method in Wireless Sensor Networks

With the integrated development of the Internet, wireless sensor technology, cloud computing, and mobile Internet, there has been a lot of attention given to research about and applications of the Internet of Things. A Wireless Sensor Network (WSN) is one of the important information technologies in the Internet of Things; it integrates multi-technology to detect and gather information in a network environment by mutual cooperation, using a variety of methods to process and analyze data, implement awareness, and perform tests. This paper mainly researches the localization algorithm of sensor nodes in a wireless sensor network. Firstly, a multi-granularity region partition is proposed to divide the location region. In the range-based method, the RSSI (Received Signal Strength indicator, RSSI) is used to estimate distance. The optimal RSSI value is computed by the Gaussian fitting method. Furthermore, a Voronoi diagram is characterized by the use of dividing region. Rach anchor node is regarded as the center of each region; the whole position region is divided into several regions and the sub-region of neighboring nodes is combined into triangles while the unknown node is locked in the ultimate area. Secondly, the multi-granularity regional division and Lagrange multiplier method are used to calculate the final coordinates. Because nodes are influenced by many factors in the practical application, two kinds of positioning methods are designed. When the unknown node is inside positioning unit, we use the method of vector similarity. Moreover, we use the centroid algorithm to calculate the ultimate coordinates of unknown node. When the unknown node is outside positioning unit, we establish a Lagrange equation containing the constraint condition to calculate the first coordinates. Furthermore, we use the Taylor expansion formula to correct the coordinates of the unknown node. In addition, this localization method has been validated by establishing the real environment.


Introduction
Precision agriculture is one of the most promising application domains where wireless sensor networks (WSN) may deliver a feasible or even optimal solution. Generally, wireless sensor networks consist of a large number of densely deployed small sensor nodes with sensing, computation, and wireless communication capabilities. Sensor nodes do not incorporate an infrastructure. They build up a network autonomously, without any external guidance or supervision. Precision agriculture is a crop and livestock production management system that uses a wireless sensor network to monitor equipment field positions to collect information. Precision agriculture technologies include equipment

Multi-Granularity Region Partition Based on RSSI
In this paper, a multi-granularity [22] region partition method is proposed. Firstly, the location method is chosen based on the received signal strength. Due to the fluctuation of strength value, Gauss fitting is introduced to estimate the value of RSSI. Secondly, it introduces a Thiessen polygon [23] using each anchor node as the center. The whole area is divided into a plurality of sub-regions, and the sub-regions of the adjacent nodes are combined to form multiple triangles, locking the unknown node into the final area.

Ranging Method Based on RSSI
The distance measurement method based on the RSSI theory model includes the log normal shielding model (Shadowing Model) and the free space propagation path-loss model. The free space propagation path-loss model is not easy to use in actual environments, however, because it belongs to the ideal state transfer model. It gives the energy consumption of the signal transmission distance in an infinite vacuum, with no influence from other factors. However, in practical application, the transmission distance of the signal is not only non-linear, but also the interference from the signal is not insignificant. Considering the influence factors of various kinds of reflection, scattering, and occlusion in the practical application environment, the attenuation of the channel is similar to the log normal distribution, and the Shadowing Model is more in line with the practical application.
Sensors 2016, 16, 1934 4 of 26 The Shadowing Model formula is as follows: where Loss indicates the signal path energy consumption, d indicates signal transmission distance, the unit is meter, n indicates the path-loss factor of the actual environment, f indicates the radio signal power, and the unit is MHz. The free space propagation path-loss model formula is as follows: where P(d) indicates the signal path-loss when the actual measurement distance is d; the unit is dBm. P(d 0 ) indicates the signal path-loss when the actual measurement distance is d 0 . The path-loss refers to the absolute power and n indicates the path-loss factor. The loss factor has different values in different environments. ε indicates the shadowing factor, the standard deviation ranges from 4 to 10, and the unit is dB. In this paper, the loss model of distance is d 0 = 1m, that is, ε ∼ N(0, δ 2 ).
where P t indicates the signal transmitting power; the unit is dBm. P(d) indicates the signal path-loss when the actual measurement distance is d. RSSI indicates the signal strength value when the receiving node distance is d 0 . Using Equations (2) and (3), the following formula can be obtained: From the above formula, when d 0 = 1m, the relationship between the intensity and the distance is as follows:

RSSI Data Processing Method
Each node repeatedly measures intensity and then collects a large number of test data as far as possible to remove the error and the noise data. To obtain the optimal intensity data, the chosen intensity data will use the wireless signal loss model to estimate the distance of the nodes.

Experimental Data Acquisition
The sensor network system used is from the WIRELESS DRAGON TECHNOLOGY COMPANY (Chengdu, China). The node uses the CC2530 module, as shown in Figure 1 In the experiment, the node receives the signal from the gateway transmission equipment. The node records the signal strength values and sends the data packets to the gateway transmission equipment. Intensity data is from different locations in the gateway node; in the whole process of the experiments, the actual distance is around 25 m from gateway to node. If it exceeds the range, because signal attenuation is too big, the measured RSSI value is not accurate. The application background of node localization has no special value.
The Table 1 shows the strength value of four nodes in different positions. In the experiment, the intensity data is 1000 sets. Statistics are performed on multiple datasets for each measurement node.

RSSI Data Processing Method
Each node repeatedly measures intensity and then collects a large number of test data as far as possible to remove the error and the noise data. To obtain the optimal intensity data, the chosen intensity data will use the wireless signal loss model to estimate the distance of the nodes.

Experimental Data Acquisition
The RSSI data can be approximated by a normal distribution. The curve of the peak shows the locations of the optimal RSSI values and the corresponding distance d is regarded as the optimal distance. From the experimental data, we can see that the RSSI data are in line with the normal distribution. In this paper the Gauss fitting method is used to select the optimal RSSI value. The Gauss fitting function is as follows: where n−1 . In the above formula, a is constant and greater than 0, b is the average strength value, and c equals the standard deviation. In this paper, we use four nodes to carry out multiple sets of measurements. Figure 2 gives the RSSI fitting curve of Node 4 at 1 m, 2 m, 3 m, and 4 m, respectively, where the abscissa is the RSSI value and the ordinate is the probability value.

Wireless Signal Transmission Loss Model
By obtaining the optimal intensity value, the distance is estimated using the wireless signal loss model. According to n = P(d) −RSSI , it can be seen that the environmental factor indicates the degree of loss in the actual transmission, and the numerical value of the signal varies with the change in the distance. Table 2 gives the n calculation results when the distance is 1.4 m, 2.2 m, 3.2 m, and 4.5 m.
In the formula, , the formula is as follows:  From the above table, we can get the result of the multi-group n value. The average value of each node is obtained, and the formula is as follows: In the formula, n node 1 indicates the environmental factor. The average n value is used and the formula is as follows: n = n node1 + n node2 + n node3 + n node4 4 .
Combining the environmental factor n, the wireless transmission loss model, the strength P(d) = 103(dBm), and the distance d 0 = 1m, the formula is as follows:

Region Division Method
A multi-granularity partition method is proposed in this paper. Firstly, according to the characteristics of Thiessen polygons, positioning areas were preliminarily divided using the anchor node location and communication. Secondly, the overlapping portions are acquired by combining the adjacent nodes triangle with a polygon area. By gradually narrowing the scope, we will lock target nodes into the possible regional positioning. Finally, the coordinate is calculated using the localization algorithm. Definition 1. Take a plane with two different points A and B. C is the perpendicular bisector of the line from A to B. The plane is divided into two parts: L right and L le f t , where A ∈ L le f t ,∃P, |PA| < |PB|, that is, the point-to-point distance from P to A is less than the point-to-point distance from P to B, namely the space located inside the plane L le f t is closer to point P, as shown in Figure 6. Definition 2. Take a plane with N different points {P 1 , P 2 , ..., P n }. For any point P i , there is a closest point distance P 2 , and the distance from P 2 to P i is less than the distance from P 2 to other points, that is, the region V p contains P i in the intersection of N − 1 half plane. The N − 1 half plane is determined by the perpendicular bisector of P i and other points, where region V p is composed of a plurality of vertical bisectors of the polygon.
By the above definition, the plane is divided into N regions V P i ; each region V p has a center, and the edges of V p are two adjacent lines of a perpendicular bisector. These two points are known Thiessen polygons, with the edges of V p between the focus and vertex. To sum up: if (x, y) ∈ V P i , P i is the closest point of the distance from (x, y).

Regional Primary Division
In this paper, the distance vector is established. According to the characteristics of Thiessen polygons, the area was preliminarily divided to narrow it down. In Figure 3, this paper establishes an overall positioning area of 10 × 10 m 2 , placing a plurality of anchor nodes. The blue point shows the position of anchor nodes; the blue line shows the perpendicular bisector of the anchor line. The closed region is a polygon using the anchor node as the center. Because each node belongs to the respective area, the whole positioning area is divided into a plurality of smaller areas. In the node communication radius, the target node U communicates with anchor nodes K 1 , K 2 , ..., K n . Sequence distance d 1 , d 2 , ..., d n are acquired between the target node U and anchor node K 1 , K 2 , ..., K n .
In Figure 3, the nearest anchor node K 1 of the target node U is calculated, and the sequence distance d k 2 , d k 3 , ..., d k s between K 1 and the adjacent nodes K 2 , K 3 , ..., K s is calculated (i ≤ 6). The formula is as follows: where (x 1 , y 1 ) and (x 2 , y 2 ) are the respective coordinates of the anchor nodes. d is the distance between K 1 and K 2 . The vertical line of the connecting line from K 1 to K 2 , K 3 , ..., K s is established, and a closed area is composed of a plurality of vertical lines. The area is a polygonal area using an anchor node K 1 as the center, called V K 1 .  In Figure 4, the nearest anchor node 1 K of the target node U is calculated, and the sequence K and the adjacent nodes ). The formula is as follows:

Region Dividing Based on Second Division
After dividing, the polygon area is still large, so a sequence table is established to choose the anchor node closest to the unknown node. At this time, the unknown node region is wide, and the final judgment coordinate sequences are not entirely accurate. Therefore, in this paper, we consider overlapping region of polygons and triangles. As far as possible we will target the minimum range, and then calculate the coordinates of the target node. The main division process is as follows: According to the region division it has been judged that U is in  To determine U: if it is in V K 1 , the area will continue division; if it is not in V K 1 , it will select sequence distance ranking second anchor node K 2 . If U is within V K 2 , it is found one V K i , which the target node U is in this region V K i . According to the nature of the Thiessen polygon that can be drawn, if the target node U is within this region V K 1 , it needs to meet the following conditions, where d (U,K 1 ) is the distance from U to K 1 , d (U,K 1 ) is the distance between U and K 1 , K 2 , ..., K n :

Region Dividing Based on Second Division
After dividing, the polygon area is still large, so a sequence table is established to choose the anchor node closest to the unknown node. At this time, the unknown node region is wide, and the final judgment coordinate sequences are not entirely accurate. Therefore, in this paper, we consider overlapping region of polygons and triangles. As far as possible we will target the minimum range, and then calculate the coordinates of the target node. The main division process is as follows: According to the region division it has been judged that U is in V K 1 . At the same time, it can get the sequence distance d 1 , d 2 , ..., d n between U and K 2 , K 3 , ..., K 5 . The distance sequence d K 2 , d K 3 , ..., d K 5 between K 1 and K 2 , K 3 , ..., K 5 . K 1 , K 2 , K 3 , ..., K 5 is composed of a number of triangles, where K 1 will be regarded as a fixed point in the triangle, and the rest of the vertices are K 2 , K 3 , ..., K 5 . The probability of target node U being in ∆K 1 K 2 K 3 is assessed according to the following description: In Figure 4, U 1 and U 2 are the target node, U 1 is in the triangle area ∆K 1 K 2 K 3 , and U 2 is outside the triangle area ∆K 4 K 5 K 6 . To determine whether a point is in the triangle, two common methods are used: area method and vector directionmethod. In this paper, the method of area is used to judge whether the target node is in the triangle, and the formula is as follows: Sensors 2016, 16, 1934 9 of 27 In Figure 5, 1 U and 2 U are the target node, 1 U is in the triangle area outside the triangle area To determine whether a point is in the triangle, two common methods are used: area method and vector directionmethod. In this paper, the method of area is used to judge whether the target node is in the triangle, and the formula is as follows: (a) (b) In Figure 6, the red dotted line is composed of triangles, while the area is represented by the continued division in Figure 5. Because each polygon area can be divided into a number of triangular regions, they can be overlapping.  In Figure 5, the red dotted line is composed of triangles, while the area is represented by the continued division in Figure 4. Because each polygon area can be divided into a number of triangular regions, they can be overlapping. In Figure 6, the red dotted line is composed of triangles, while the area is represented by the continued division in Figure 5. Because each polygon area can be divided into a number of triangular regions, they can be overlapping.

Analysis of Regional Division Method
The region partition process includes main two parts. Firstly, the whole region is divided into a plurality of sub-regions using the Thiessen polygon method. Secondly, a combination of adjacent nodes forms a triangular region. Lastly comes the judgment of the final area for unknown nodes. The area partition method pseudo codes (Algorithm 1) are as follows:

Analysis of Regional Division Method
The region partition process includes main two parts. Firstly, the whole region is divided into a plurality of sub-regions using the Thiessen polygon method. Secondly, a combination of adjacent nodes forms a triangular region. Lastly comes the judgment of the final area for unknown nodes. The area partition method pseudo codes (Algorithm 1) are as follows: Input: A set of all anchor nodes: N = {N 1 , N 2 ..., N n } j unknown nodes evenly distributed in the location region; Output: Localization algorithm for unknown node P 1. A set of anchor nodes N marked serial number for n anchor nodes 2. FOR (i = 1 to j) 3. RUN P_RSSI(i), obtaining the strength vector set between N and P P_RSSI = {R 1 , R 2 , ..., R n }. 4.
O_d(i)_top 6, taking the top six anchor nodes A, B, C, D, E, and F, using the vertical line with the O midline forming a closed polygon area V O . 12.
IF P_d O < P_d k (k = A, B, C, D, E, F), judge whether the distance between P and O is less than the other adjacent anchor nodes, O, B, C, D, E, and F. 13.
RUN the V O neighboring anchor nodes A, B, C, D, E, and F following the order of O_d(i)_top 6 forming a number of triangles ∆Oij(i, j = A, B, C, D, E, F). 14.
IF , that is, the point-to-point distance from P to A is less than the point-to-point distance from P to B , namely the space located inside the plane left L is closer to point P , as shown in Figure 3.

Regional Primary Division
In this paper, the distance vector is established. According to the characteristics of Thiessen polygons, the area was preliminarily divided to narrow it down. In Figure 4, this paper establishes an overall positioning area of 2 10 10 m  , placing a plurality of anchor nodes. The blue point shows the position of anchor nodes; the blue line shows the perpendicular bisector of the anchor line. The closed region is a polygon using the anchor node as the center. Because each node belongs to the respective area, the whole positioning area is divided into a plurality of smaller areas. In the node

Node Localization Algorithm Based on Lagrange Multiplication and Taylor Formula
This paper uses a Shadowing Model to show the relationship between the signal intensity and the distance. It uses the Gaussian fitting processing intensity data, and acquires the distance. Firstly, the distance sequence is introduced to select the nearest anchor node. According to the characterization of the Thiessen polygon, the anchor nodes are divided into polygonal regions. Secondly, the neighboring nodes are chosen around anchor nodes. The polygon region division is composed of the triangle using multiple adjacent nodes, so as to determine the unknown node positioning area. Thirdly, by calculating the distance from the virtual reference point to the anchor nodes, the distance vector sequence from reference point to anchor node is established. According to the characteristics of vector similarity, several recent reference points are chosen. Finally, the centroid algorithm is used to calculate the coordinates of unknown nodes. On the other hand, considering the location of edge nodes, if unknown nodes are not in the overlapping region, a close reference point is chosen using the characteristics of vector similarity. Furthermore, the coordinates of the initial node are used in the Lagrange multiplication [24] and Taylor's formula [25]. By iterative obtainedcoordinates, it may acquire rational coordinates of the unknown node.

Node Location Method in Positioning Unit
The localization idea is acquiring location by dividing the overall region. This can improve the nodes in the edge region, namely non-overlapping regions. Nodes in the interior of the unit are as follows. Firstly, according to the nature of the Thiessen polygon, the unknown node is within this area. Finally, the appropriate virtual reference point is selected to compute the coordinates of the target node.
According to Figure 7, the location of the unknown node in the unit is as follows:

Node Location Method in Positioning Unit
The localization idea is acquiring location by dividing the overall region. This can improve the nodes in the edge region, namely non-overlapping regions. Nodes in the interior of the unit are as follows. Firstly, according to the nature of the Thiessen polygon, the unknown node is within this area. Finally, the appropriate virtual reference point is selected to compute the coordinates of the target node.
According to Figure 7, the location of the unknown node in the unit is as follows: In the location area, distance d is estimated by RSSI. Furthermore d is given in ascending order and the ranked the first anchor node must be selected. According to the Thiessen polygon theorem, the polygonal regions of the anchor node are divided. As shown in Figure 7, O represents the unknown node. A, B, C, D, and E represent the anchor nodes. The distance sequence of O, from small to large, is: A→B→C→D→E.
To determine whether the unknown nodes are in the polygon, the specific method is as follows. If the unknown nodes are judged to be within this region, complete the following steps. Otherwise, ranked the second anchor node is selected to judge whether or not the unknown nodes are in this triangle. As shown in Figure 7, O is in the region using A as the center, that is, within A V .
Coordinates are acquired in the polygon area. The d sequence between the unknown nodes and the nearby anchor node is acquired. The d values are in ascending order. The number is ranked from 1 to 7 (  The unknown nodes are in the triangle area. By calculating the distance from each virtual reference node to the anchor nodes, corresponding distance vectors are formed. As shown in Figure 7, O is in VA and in ΔABC (where a, b, c, and d are virtual reference nodes). These reference points are pairwise linear points or points of intersection. Their coordinates can be obtained.
According to the distance vector similarity, we can select 3 the highest similar degree virtual reference node. As shown in Figure 7, the midpoint may be obtained by calculating from the virtual In the location area, distance d is estimated by RSSI. Furthermore d is given in ascending order and the ranked the first anchor node must be selected. According to the Thiessen polygon theorem, the polygonal regions of the anchor node are divided. As shown in Figure 7, O represents the unknown node. A, B, C, D, and E represent the anchor nodes. The distance sequence of O, from small to large, is: A→B→C→D→E.
To determine whether the unknown nodes are in the polygon, the specific method is as follows. If the unknown nodes are judged to be within this region, complete the following steps. Otherwise, ranked the second anchor node is selected to judge whether or not the unknown nodes are in this triangle. As shown in Figure 7, O is in the region using A as the center, that is, within V A .
Coordinates are acquired in the polygon area. The d sequence between the unknown nodes and the nearby anchor node is acquired. The d values are in ascending order. The number is ranked from 1 to 7 (i ≤ 6). The ranked the first anchor node of d is acquired. According to the d sequence, the other anchor nodes form multiple triangles. As shown in Figure 7, O is in V A , and neighboring nodes of V A are B, C, D, and E; thus, in accordance with the principle of distance sequence, various triangular combinations exist: ∆ABC, ∆ABD, ∆ABE, ∆ACD, ∆ACE, and ∆ADE.
The unknown nodes are in the triangle area. By calculating the distance from each virtual reference node to the anchor nodes, corresponding distance vectors are formed. As shown in Figure 7, O is in VA and in ∆ABC (where a, b, c, and d are virtual reference nodes). These reference points are pairwise linear points or points of intersection. Their coordinates can be obtained.
According to the distance vector similarity, we can select 3 the highest similar degree virtual reference node. As shown in Figure 7, the midpoint may be obtained by calculating from the virtual reference nodes a, b, c, and d to anchor node A and then finding the area of the overlap. Finally, the coordinates of unknown nodes are calculated by a centroid algorithm.

Determining the Location Area of Node
As shown in Figure 8a, P represents the unknown node. A, B, and D represent three anchor nodes. The polygon region within the solid lines represents a polygon region with D as the center, called V D . The dotted area represents ∆ABD. The shaded area represents the overlap of two regions, V D and ∆ABD. a and b represent the midpoint from D to the neighboring nodes A and B, respectively. C represents the intersection point. V D is intersected by two edges, forming a solid line area. The A, B, and D coordinates are known: (x A , y A ), (x B , y B ), and (x D , y D ). If a and b are the midpoints, you can get the virtual reference point a, b coordinates; the formula is as follows: Also, the intersection C is line DA L and DB L perpendicular.
As shown in Figure 8b, P represents the unknown node. A, B, D, and E represent four anchor nodes. The area within the solid lines represents a polygon region with D as the center, called D V .
The dotted area represents ΔADE. The shading represents region D V and ΔADE in the overlapping region. The distance sequence from P to A, B, D, and E is ordered from small to large as follows: D→E→A→B, thus judging that P is in ΔADE, where a and b represent the midpoint between D and A or E, respectively.
where the values of the parameters are, respectively, (a) One of overlapping area (b) Other overlapping area  As shown in Figure 8b, P represents the unknown node. A, B, D, and E represent four anchor nodes. The area within the solid lines represents a polygon region with D as the center, called V D . The dotted area represents ∆ADE. The shading represents region V D and ∆ADE in the overlapping region. The distance sequence from P to A, B, D, and E is ordered from small to large as follows: D→E→A→B, thus judging that P is in ∆ADE, where a and b represent the midpoint between D and A or E, respectively. C represents the intersection by the perpendicular bisector of the line through L DA and L DB .
The overlap region does not select coordinates of the intersection, as shown in Figure 8b. e represents line L AE and L DA the perpendicular bisector intersection. f represents line L AE and L DA the perpendicular bisector intersection. Selecting intersections e and f determines the overlapping region more accurately. It may guarantee a point P being in ∆ADE. The coordinates of the intersection are as follows: where the values of the parameters are, respectively, b ae = (x a y a +1)(y A −y D ) Some virtual reference points can be obtained for the overlapping area. The coordinates of the reference points are for the midpoint Equation (13) obtained, as shown in Figure 9. Anchor nodes are A, D, and E; virtual reference point h is the midpoint between D and B. The small red circles are virtual reference points and the black triangle represents the unknown node.

Selecting Virtual Reference Node
There is an unknown node P whose coordinates are ) , ( y x ; reference nodes P are A, B, and C, and their coordinates are , respectively. P is in ABC Δ ; the distance can be obtained by the following equation: In the above equation group, the coordinates ) , ( y x are judged by the distance d.

Selecting Virtual Reference Node
There is an unknown node P whose coordinates are (x, y); reference nodes P are A, B, and C, and their coordinates are (0, 0), (1, 1), (0, 2), respectively. P is in ∆ABC; the distance can be obtained by the following equation: In the above equation group, the coordinates (x, y) are judged by the distance d.

Lemma 1.
In the positioning unit, there are a number of unknown nodes P i (i = 0, 1, 2, ..., n). The distance vector between any P i and the reference node is corresponding to the coordinates of that point, moreover the value of the coordinates is only one.
Definition 3. In the regional location, there is a sample point, denoted as x i . The sample point can receive the other node strength RSSI value and form a RSSI vector set sorted in descending order, i.e., R 1 = {r 11 , r 12 , ..., r 1n }, where r 1 denotes the sample x 1 collecting the strength value from node 1 and R 1 denotes the sample x 1 collecting the strength value of all nodes.

Definition 4.
In the location area, there are multiple samples x 1 , x 2 , ..., x n , and each sample point can receive the other node strength RSSI value. From the many groups an RSSI vector set is formed and sorted in descending order, i.e., R = {R 1 , R 2 , ..., R n }, where R 1 represents samples x 1 collecting the intensity values from all the nodes and R represents the samples x 1 , x 2 , ..., x n collecting the strength set from all nodes.
According to Definitions 3 and 4, the RSSI vector set is as follows: According to Equation (16), we can get the corresponding distance vector set: where d 11 represents the distance vector between the sample point x 1 and the node 1; D 1 represents the distance vector of the sample points x 1 between and all nodes; and D represents the distance vector set from all the sample points x 1 , x 2 , ..., x n to all nodes. The distance from each unknown node to the anchor nodes produces a set of vectors. Using fuzzy mathematics, we can examine the close degree of unknown nodes, called vector similarity. In this paper, we establish the distance vector between the unknown node and the anchor node, and the distance vector between the reference node and the anchor node. The vector cosine similarity [26] is used to select the nearest reference node.
Vector cosine similarity is used to measure the degree of individual differences; the parameter is the cosine of the angle between the two vectors. The formula is as follows: If A and B are two n-dimensional vectors, A is [A 1 , A 2 , ..., A n ], B is [B 1 , B 2 , ..., B n ], the cosine value of angle θ between A and B is calculated as follows: As shown in Figure 10, A, B, D, and E represent anchor nodes, and the black triangle represents unknown node P. At the same time, it has been determined that P is in V D and also in ∆DAE. Next the virtual reference point is acquired. A red point represents a virtual reference point; reference points h, b, and i are close to the virtual reference point.

 
As shown in Figure 10, A, B, D, and E represent anchor nodes, and the black triangle represents unknown node P. At the same time, it has been determined that P is in D V and also in DAE Δ . Next the virtual reference point is acquired. A red point represents a virtual reference point; reference points h, b, and i are close to the virtual reference point. The main selection process is as follows. Firstly, P collects the RSSI from A, D, and E and forms the vector set in descending order. . The formula is as follows: The θ cos value of the formula is close to 1, which indicates that the angle is close to 0, that is, the two vectors are more similar. The distance similarity between virtual reference nodes and The main selection process is as follows. Firstly, P collects the RSSI from A, D, and E and forms the vector set in descending order. Then it may estimate each RSSI value corresponding to the distance value d using a wireless signal loss model. This can also form a distance vector {d PD , d PE , d PB } using A, D, and E. Secondly, in selected areas we can compute the distance from all virtual reference points to A, D, and E. Using the Euclidean distance formula we can obtain each virtual reference point in the A, D, E distance vector, called {d iD , d iE , d iB }. Before the distance vector is formed, each distance value is sorted according to the corresponding distance. In this paper, the rules for anchor nodes are sorted according to the distance, where d iD represents the distance between virtual reference point i and anchor node D. Lastly, the distance vector sequence similarity is used to compute the similarity degree between {d PD , d PE , d PB } and {d iD , d iE , d iB }. The formula is as follows: The cosθ value of the formula is close to 1, which indicates that the angle is close to 0, that is, the two vectors are more similar. The distance similarity between virtual reference nodes and unknown nodes are computed. Furthermore, similarity is given in descending order and the virtual reference point is used when ranking the top three. The P distance vector and h, b, i distance vector are similar, that is, h, b, i was close to P's three virtual reference nodes.
Finally, the centroid algorithm is used to calculate the coordinates of unknown nodes P. Figure 8 is an example and the formula is as follows:

Determining the Range of Node Coordinates
As shown in Figure 11, P represents the unknown node, A represents the nearest anchor node from the P, and P is in the region V A . B represents the nearest virtual reference point to P.
Sensors 2016, 16,1934 15 of 27 unknown nodes are computed. Furthermore, similarity is given in descending order and the virtual reference point is used when ranking the top three. The P distance vector and h, b, i distance vector are similar, that is, h, b, i was close to P's three virtual reference nodes. Finally, the centroid algorithm is used to calculate the coordinates of unknown nodes P. Figure 8 is an example and the formula is as follows:

Determining the Range of Node Coordinates
As shown in Figure 11, P represents the unknown node, A represents the nearest anchor node from the P, and P is in the region A V . B represents the nearest virtual reference point to P. The steps for determining the coordinate range are as follows: Step 1. According to the scope of the regional positioning, the boundary equation is acquired, as well as the computing distance from the center of the polygon region to the boundary. As in Figure 11, The steps for determining the coordinate range are as follows: Step 1. According to the scope of the regional positioning, the boundary equation is acquired, as well as the computing distance from the center of the polygon region to the boundary. As in Figure 11, the boundary of location area are x = 0, x = 10, y = 0 and y = 10, respectively, and we can calculate the distance from A to x = 0, x = 10, y = 0 and y = 10 in ascending rank order.
Step 2. By choosing the distance value, we can get the distance of the corresponding boundary equations, and the value of this equation is considered as one of the candidate coordinates. For example, the minimum distance value d A2 , which represents the distance between A and x = 10. Then x = 10 will be considered as a candidate.
Step 3. Taking A and B coordinate values and selecting candidate values using the above step, all values are compared, such as There are three coordinates being compared, i.e., x A , x b , x. The comparative value of y are two, i.e., y A , y b . The number of comparative value is determined by candidate value.
Step 4. Given the comparative values of x, y, and the x min = x A , x max = 10, y min = x A , y min = x b , the final coordinates are determined, that is, (x min ≤ x ≤ x max , y min ≤ y ≤ y max ).

Regional Primary Division
In this paper, given the two-variables function f (x, y) = (x i − x) 2 + (y i − y) 2 − d 2 i , the minimum value is computed, namely, the minimal difference between the actual distance and the calculated distance.
As shown in Figure 12, P represents the unknown node, A represents the nearest anchor node to B, and b represents the virtual reference point. The rectangular region represents P's location region, that is, the coordinate range of P. The calculated coordinate process is shown as follows: Sensors 2016, 16,1934 16 of 27 Given the two-variables function is computed, namely, the minimal difference between ≥ y x f , that is, ; the formula is as follows: The node coordinate calculation model is established by the Lagrange multiplier method as follows: In the above formula, where y x , represents the coordinates of the unknown nodes, and ε represents the sum of the range ) , ( y x , that is, Using the Lagrange multiplier method, the general equation is obtained as follows: According to the above Equations (24) and (25), the partial derivative of x, yis computed Given the two-variables function f (x, y) = ( y) is computed, namely, the minimal difference between (x i − x) 2 + (y i − y) 2 and d 2 i . In order to simplify the calculations, this paper only discusses the extreme value of the function f (x, y) ≥ 0, that is, (x i − x) 2 + (y i − y) 2 ≥ d 2 ; the formula is as follows: In the above formula, (x A , y A ) are the nearest anchor node coordinates from unknown nodes; the distance between A and P is expressed as and (x, y) are the unknown node's coordinates.
The node coordinate calculation model is established by the Lagrange multiplier method as follows: L(x, y) = f (x, y) + λ * ϕ(x, y).
In the above formula, f (x, y) represents the two-variables function about anchor node coordinates, unknown node coordinate, measuring distance; ϕ(x, y) represents the unknown node area, namely node coordinate constraint condition. L(x, y) represents the positioning function model and λ represents the parameter. By the constraint of node coordinates (x min ≤ x ≤ x max , y min ≤ y ≤ y max ), the constraint equation is established as follows: where x, y represents the coordinates of the unknown nodes, and ε represents the sum of the range (x, y), that is, ε = (x max − x min ) + (y max − y min ). Using the Lagrange multiplier method, the general equation is obtained as follows: According to the above Equations (24) and (25), the partial derivative of x, y is computed respectively. Let the partial derivative equal zero and combine the constraint Equation (24) establishing a set of equations, as follows: By solving the equations, we can get the value x, y, λ. At this time x, y are the extreme point coordinates, that is, the extreme points in the function z = f (x, y) in constrained conditions. The above figure is an example, using the Lagrange multiplier method; the full equations are as follows: Solving the above equations, x, y, λ are acquired. Firstly, by using the method of Lagrange, the initial coordinates (x, y) are obtained; (x n , y n ) indicates the coordinates of the anchor nodes. f (x, y) is determined as follows: The formula is expanded in the (x, y) Taylor series; the modifying step size k and only the first-order partial derivative are retained, and the higher order functions are ignored: From the above calculations, the coordinates of the range are (x min ≤ x ≤ x max , y min ≤ y ≤ y max ). The fixed step size k, h range is determined, that is, (x min − x ≤ h ≤ x max − x, y min − y ≤ k ≤ y max − y), and then the threshold is acquired, Combining the two equations above, calculation step correction value, and judgment h, k meets the threshold. If it meets the threshold, the current coordinates are the final coordinates of the unknown nodes. If it does not satisfy the threshold, iteration step size h, k is modified until it meets the threshold, finally returning to the current coordinates.

Simulation Results
Anchor nodes provide different coverage and thus different simulation environments are established. The experimental platform is MATLAB 7.0 and MyEclipse 9.0. First of all, a 10 × 10 m 2 square area simulation environment is established. The anchor node number is 8, 16, 24, or 32. Unknown nodes are evenly distributed throughout the whole positioning region. Assume that all unknown nodes can communicate with the anchor nodes.
As shown in Figure 13, the red dot represents a randomly generated unknown node; the blue point represents the anchor node coordinates. Taking the number of anchor nodes to be 32 for an example, shown in Figure 14, the regional positioning is as below, where a blue point represents a location within all anchor nodes and solid lines enclose a polygon region. meets the threshold, finally returning to the current coordinates.

Simulation Results
Anchor nodes provide different coverage and thus different simulation environments are established. The experimental platform is MATLAB 7.0 and MyEclipse 9.0. First of all, a 2 10 10 m  square area simulation environment is established. The anchor node number is 8, 16, 24, or 32. Unknown nodes are evenly distributed throughout the whole positioning region. Assume that all unknown nodes can communicate with the anchor nodes. As shown in Figure 13, the red dot represents a randomly generated unknown node; the blue point represents the anchor node coordinates. Taking the number of anchor nodes to be 32 for an example, shown in Figure 14, the regional positioning is as below, where a blue point represents a location within all anchor nodes and solid lines enclose a polygon region. Figure 13. One of the simulation location area distribution maps. Figure 13. One of the simulation location area distribution maps. This paper uses original strength value fluctuation of 5% and 10% in simulating. Taking the original coordinates of unknown nodes (6.5, 4.9) as an example, shown in Figure 15, the blue represents all the anchor nodes in the location area. Solid lines represent each anchor node in the center partition polygon area, in which a, b, c, and d represent anchor nodes with respective coordinates of (7.0, 4.0), (7.0, 6.0), (5.0, 6.0), and (5.0, 4.0). The dotted lines form a four-triangle area. The small red dot represents the original position of the unknown nodes, with coordinates of (6.5, 4.9), The black squares represent intensity fluctuation, being 5% at the estimated coordinates. Green dots represent the estimated node coordinates, when the intensity fluctuation of 10%. With intensity fluctuations, the node localization results are Tables 3 and 4. This paper uses original strength value fluctuation of 5% and 10% in simulating. Taking the original coordinates of unknown nodes (6.5, 4.9) as an example, shown in Figure 15, the blue represents all the anchor nodes in the location area. Solid lines represent each anchor node in the center partition polygon area, in which a, b, c, and d represent anchor nodes with respective coordinates of (7.0, 4.0), (7.0, 6.0), (5.0, 6.0), and (5.0, 4.0). The dotted lines form a four-triangle area. The small red dot represents the original position of the unknown nodes, with coordinates of (6.5, 4.9), The black squares represent intensity fluctuation, being 5% at the estimated coordinates. Green dots represent the estimated node coordinates, when the intensity fluctuation of 10%. With intensity fluctuations, the node localization results are Tables 3 and 4. represents all the anchor nodes in the location area. Solid lines represent each anchor node in the center partition polygon area, in which a, b, c, and d represent anchor nodes with respective coordinates of (7.0, 4.0), (7.0, 6.0), (5.0, 6.0), and (5.0, 4.0). The dotted lines form a four-triangle area. The small red dot represents the original position of the unknown nodes, with coordinates of (6.5, 4.9), The black squares represent intensity fluctuation, being 5% at the estimated coordinates. Green dots represent the estimated node coordinates, when the intensity fluctuation of 10%. With intensity fluctuations, the node localization results are Tables 3 and 4.     In this paper, the error formula is as follows:     In this paper, the error formula is as follows: where (x 0 , y 0 ) are the coordinates, the actual measurement coordinates is (x, y), and the communication radius is R. As shown in Figures 16 and 17, the abscissae of the two charts reflect the number of anchor nodes: N = 8, N = 16, N = 24, N = 32, respectively.
Sensors 2016, 16,1934 19 of 27      Figure 18 shows the localization algorithm, SBL [18], three-centroids in a triangle with a sequence algorithm [19], and three-centroids in a triangle [14]; the error trends are correlated with the number of anchor nodes. The position error curve basically showed a decline. More than 20 unknown nodes are randomly generated, and calculate the coordinates of these nodes, when the number of anchor nodes are   Figure 18 shows the localization algorithm, SBL [18], three-centroids in a triangle with a sequence algorithm [19], and three-centroids in a triangle [14]; Taking into account the different coverage of node communication, the following simulation environment is the region 2 100 100 m  . The unknown node is evenly distributed in the positioning area. As shown in Figure 19, the red dots indicate the random distribution of the unknown nodes; blue points indicate the anchor nodes. Taking into account the different coverage of node communication, the following simulation environment is the region 100 × 100 m 2 . The unknown node is evenly distributed in the positioning area. As shown in Figure 19, the red dots indicate the random distribution of the unknown nodes; blue points indicate the anchor nodes. Taking into account the different coverage of node communication, the following simulation environment is the region 2 100 100 m  . The unknown node is evenly distributed in the positioning area. As shown in Figure 19, the red dots indicate the random distribution of the unknown nodes; blue points indicate the anchor nodes.  Taking the number of anchor nodes to be N = 16, N = 24, N = 32 respectively, the location of the unknown node is computed. With the fluctuation of 5% or 10%, the unknown node is judged to be in a different location. According to the area of the unknown node, this may select the appropriate positioning algorithms, namely the unknown node inside or outside positioning unit method. Finally, part of the results for the positioning of unknown nodes are shown in Tables 5-7.  In the 100 × 100 m 2 positioning area, if a different number of anchor nodes are distributed, each anchor node area's coverage is different, so the regional positioning division scope changes accordingly.
From Figures 20 and 21, we see that the error of the unknown nodes decreases with the increase in the number of anchor nodes when the intensity data fluctuation is 5% and 10%. When the number of anchor nodes is 16, 24, or 32, the average location error of the unknown nodes is 12.55 m, 7.49 m, or 4.11 m. In the 100 × 100 m 2 regional location, because of the greater area of anchor nodes' partitioning, the unknown node localization becomes larger, thus affecting the coordinates of unknown nodes. Compared to the effect of location, the positioning error is increased.

Experimental Results
The sensor node communication range is less than 25 m. In order to further study the feasibility, we use our localization algorithm in the real environment. In this paper, the ZigBee wireless sensor network system is used, which is provided by the Wireless Dragon Technology Company (Chengdu, China). The experimental equipment includes four anchor nodes, three unknown nodes, and gateway transmission equipment.

Experimental Results
The sensor node communication range is less than 25 m. In order to further study the feasibility, we use our localization algorithm in the real environment. In this paper, the ZigBee wireless sensor network system is used, which is provided by the Wireless Dragon Technology Company (Chengdu, China). The experimental equipment includes four anchor nodes, three unknown nodes, and gateway

Experimental Results
The sensor node communication range is less than 25 m. In order to further study the feasibility, we use our localization algorithm in the real environment. In this paper, the ZigBee wireless sensor network system is used, which is provided by the Wireless Dragon Technology Company (Chengdu, China). The experimental equipment includes four anchor nodes, three unknown nodes, and gateway transmission equipment. Figure 22 shows the experimental environment for the laboratory: in the 5 × 5 m 2 area of the four anchor nodes, the black squares are the anchor nodes, for which the coordinates are (2, 1), (5,2), (3,3), and (1, 4), respectively, and the red squares are the unknown nodes' positions. Firstly, the anchor nodes are placed in a fixed position. The RSSI data are analyzed by Gaussian fitting. The peak value is selected as the most appropriate strength. Taking 3 meter intensity data between Node 1,Node 4 and the unknown nodes for example, shown in Table 8, these data calculate the probability density of each intensity, and the Gaussian fitting strength obtains the appropriate RSSI value, as shown in Figure 23.   The main experimental process and steps are as follows: Firstly, the anchor nodes are placed in a fixed position. The RSSI data are analyzed by Gaussian fitting. The peak value is selected as the most appropriate strength. Taking 3 meter intensity data between Node 1,Node 4 and the unknown nodes for example, shown in Table 8, these data calculate the probability density of each intensity, and the Gaussian fitting strength obtains the appropriate RSSI value, as shown in Figure 23. Firstly, the anchor nodes are placed in a fixed position. The RSSI data are analyzed by Gaussian fitting. The peak value is selected as the most appropriate strength. Taking 3 meter intensity data between Node 1,Node 4 and the unknown nodes for example, shown in Table 8, these data calculate the probability density of each intensity, and the Gaussian fitting strength obtains the appropriate RSSI value, as shown in Figure 23.    Secondly, the wireless signal loss model is used to estimate the distance. Taking the original unknown node to be (1, 3) as an example, the strengths of the four anchor nodes are collected and the estimated distance values are as shown in Table 9. Thirdly, for the regional division method, the whole positioning region is divided into polygons, and then the adjacent nodes are selected to be combined into a plurality of triangles. If the original coordinates of the unknown nodes are (1, 3), (2, 2), (3, 0), the nodes' positioning unit areas are as shown in Figure 24: the blue solid line shows the anchor nodes dividing the central polygon region; the gray dotted area consists of neighboring nodes, shown as multiple triangles; the blue squares represent anchor nodes' coordinates, and the red squares represent the original coordinates of the unknown nodes.  Thirdly, for the regional division method, the whole positioning region is divided into polygons, and then the adjacent nodes are selected to be combined into a plurality of triangles. If the original coordinates of the unknown nodes are (1, 3), (2, 2), (3, 0), the nodes' positioning unit areas are as shown in Figure 24: the blue solid line shows the anchor nodes dividing the central polygon region; the gray dotted area consists of neighboring nodes, shown as multiple triangles; the blue squares represent anchor nodes' coordinates, and the red squares represent the original coordinates of the unknown nodes. As shown in Figure 25, in the division of the location of the region, the blue squares represent the anchor nodes' coordinates. The red squares represent the original coordinates of the unknown nodes. The green squares represent the unknown nodes' coordinates. By measuring signal strength in the real environment, the model is fit to select a reasonable strength value. According to the positioning algorithm, the model selects the corresponding calculation method. Finally, we get unknown nodes in the 2

5m
 region, for which the average localization error is 1.089 m. As shown in Figure 25, in the division of the location of the region, the blue squares represent the anchor nodes' coordinates. The red squares represent the original coordinates of the unknown nodes. The green squares represent the unknown nodes' coordinates. By measuring signal strength in the real environment, the model is fit to select a reasonable strength value. According to the positioning algorithm, the model selects the corresponding calculation method. Finally, we get unknown nodes in the 5 × 5 m 2 region, for which the average localization error is 1.089 m. Figure 24. Sketch map of regional divisions.
As shown in Figure 25, in the division of the location of the region, the blue squares represent the anchor nodes' coordinates. The red squares represent the original coordinates of the unknown nodes. The green squares represent the unknown nodes' coordinates. By measuring signal strength in the real environment, the model is fit to select a reasonable strength value. According to the positioning algorithm, the model selects the corresponding calculation method. Finally, we get unknown nodes in the 2

5m
 region, for which the average localization error is 1.089 m.  Because indoor interference factors cannot be absolutely avoided, the signal intensity error cannot be completely eliminated. Estimating the distance value will result in errors and the final node coordinate calculation data will be affected. The original coordinates of the unknown node are (1, 3), (2, 2), (3, 0), and (3,2). Using this algorithm, the average coordinates and positioning error are calculated; the specific data are shown in Table 10.

Conclusions
In this paper, a ranging method based on RSSI-combined with multi-granularity regional division aiming at the area of the unknown nodes, Lagrange multiplication, and a Taylor expansion location algorithmis proposed and the algorithm is verified by simulations and experiments.
This paper researches the node localization algorithm based on RSSI. Intensity data collected are analyzed in the experiment. It is found that data change according to the different distances. The RSSI data values are basically in accord with normal distribution. Using a Gaussian fitting function, the intensity data collected are fitted, as far as possible, to eliminate interference and abnormal situations and choose a fitting peak that gives the optimal intensity values. Moreover, according to the multi-granularity localization, this paper establishes the dividing region method, that is, the overlapping area of the two polygons is regarded as the final positioning region, combining the Thiessen polygon and the triangle.
Different solutions are adopted to find the locations of the unknown nodes. When unknown nodes are inside positioning unit, a candidate virtual reference node is chosen using vector similarity. Three virtual reference points from the unknown node are selected and using the centroid algorithm we can calculate the final coordinates. When unknown nodes are outside of the positioning unit, the first step is to determine where the unknown node of the polygon area is. We also used the vector similarity method to select the appropriate virtual reference node. The Lagrange method was used to establish the calculation formula and substitute the preliminary coordinates into the Taylor series expansion. By gradually modifying the coordinates, reasonable coordinate values are finally acquired.