NODE POSITIONING METHOD FOR MARINE ENVIRONMENTAL MONITORING WIRELESS SENSOR NETWORKS

Abstract

A node positioning method for marine environmental monitoring wireless sensor networks is provided, which jointly estimates node positions and path loss factors by considering the node real-time movement and the path loss and absorption effect of underwater communication; transforms an original non-convex problem into a non-negative constrained least squares framework, and finds the optimal solution of marine node positions by two stages of interior point method and block coordinate update. In the first stage, the problem is re-expressed by using the penalty function according to the interior point method to obtain the approximate solution; and in the second stage, the original problem is transformed into a generalized trust region sub-problem, and the approximate solution obtained by the interior point method is used as the initial estimation, the accurate estimation values of the marine node position and path loss factor are obtained by iterative solutions combined with the block coordinate update.

Description

CROSS-REFERENCE TO RELATED APPLICATIONS

This application claims priority to Chinese Patent Application No. 202310724196.8, filed on Jun. 19, 2023, the contents of which are hereby incorporated by reference for all purposes.

TECHNICAL FIELD

The application relates to a technical field of marine wireless sensor network node localization, and in particular to a high-precision and high-robustness positioning method for marine environmental monitoring wireless sensor networks.

BACKGROUND

In recent years, the marine environment is continuously deteriorated and marine resources are in shortage increasingly, it has brought great challenges to monitoring and managing the marine environment, and the effective monitoring and management of the marine environment has become particularly important. Marine environmental monitoring wireless sensor networks (MEM-WSNs), as a new marine monitoring technology, has a broad application prospect in marine science and resource management. However, the monitoring of the marine environment depends on monitoring data collected in accurate node coordinates, and the lack of position data makes the collected monitoring data meaningless. Therefore, obtaining accurate position information of marine nodes is the basis and guarantee to realize marine environmental monitoring. MEM-WSNs includes anchor nodes deployed on the sea surface and unknown nodes deployed underwater. Among them, the known position information of the anchor nodes is obtained through the built-in GPS or Beidou positioning module; however, the position information of unknown underwater nodes needs to be obtained indirectly, that is, the received signal strength indication (RSSI) and other data information are exchanged with the anchor nodes through an underwater acoustic communication to calculate their own position coordinates.

MEM-WSNs nodes mainly rely on underwater acoustic channels to communicate underwater, which has unique technical characteristics and difficulties: first, the attenuation of the underwater acoustic channels increases with the frequency index of acoustic waves, resulting in narrow bandwidth, low communication rates and short transmission distances of the underwater acoustic communication: second, in the marine environment, a signal transmission not only has attenuation loss, but also has absorption loss, which have a great impact on the positioning accuracy: in addition, the stress mode of nodes is complex, and the inherent dynamic characteristics in the ocean, such as vortex, internal wave and thermocline, lead to great uncertainty in the motion mode of nodes, that is, the topological structure of the network is highly dynamic, which makes it difficult to accurately locate MEM-WSNs nodes.

China patent application CN113242512A discloses a node positioning scheme using the method of sound ray compensations to correct the ranging between ordinary nodes and anchor nodes, introduces the crossover in the genetic algorithm, and embeds the crossover into the particle swarm algorithm to obtain a node position scheme for the optimal position of particles and particle swarms. China patent application CN115038165A discloses a target position estimation method for constructing a ranging model of received signal strength based on Snell's law and a ray tracing theorem according to the layered propagation effect of acoustic signals underwater. China patent application CN104302001A discloses a scheme in which each node updates its own movement mode according to the water flow movement to predict the future position of the node.

For the positioning of MEM-WSNs nodes, the existing positioning methods all have at least one of the following defects: 1) assuming that the path loss factor is known, in the application of MEM-WSNs, the path loss factor is unknown and changes in real time with the dynamic change of the marine environment: 2) the scenario that all nodes move in real time is not considered; and 3) the attenuation loss and the absorption effect of the signal transmission in the underwater environment are not considered. Therefore, it is necessary to propose a new node positioning method to realize a node positioning with high-precision and high-robustness for MEM-WSNs.

SUMMARY

In order to solve the problem of low positioning accuracy caused by the highly dynamic underwater environment and the underwater absorption effect of acoustic communication, a high-precision and high-robustness positioning method for marine environmental monitoring wireless sensor networks (MEM-WSNs) is provided in the application.

The purpose of the application may be achieved by following a technical scheme:

• • a node positioning method for MEM-WSNs, including following steps:
  • establishing a received signal strength model according to a signal propagation loss and an underwater sound propagation absorption effect;
  • establishing an optimization function by a joint estimation of a path loss factor and a target position;
  • re-expressing the optimization function by using a penalty function, and finding an approximate solution of position coordinates of marine nodes by using an interior point method; and
  • transforming an original positioning problem into a generalized trust region sub-problem, and using the approximate solution found by the interior point method as an initial estimated value of the generalized trust region sub-problem, and finding an optimal solution by combining a block coordinate update.

Further, specific steps of establishing the received signal strength model specifically include:

• • according to a wireless signal loss propagation model, a signal strength of information sent by a node received by node i is expressed as:

$L_i^t=L_s^t-L\left(d_0\right)-10\beta {\log }_{10}\frac{x^t-a_i^t}{l_0}-{\alpha }_fx^t-a_i^t+n_i$

• • where L_i^t, represents a power of a target node received by an i^th anchor node at a moment t; L_s^t, represents a transmit power of the target node at the moment t; L(l₀) represents a loss value when a reference distance is l₀: β represents the path loss factor: ∥□∥ represents a second-order norm: n_i represents a signal attenuation noise of the i^th anchor node; α_f represents an absorption coefficient: a_i^t represents a position of the i^th anchor node at the moment t; and x^t represents a position of the target node at the moment t.

Further, the joint estimation of the path loss factor and the target position is as follows:

$F\left({\hat{x}}^t,\beta \right)=\underset{x^t,\beta }{\arg \min }\sum _{i=1}^M\frac{\left(L_i^t-L_0+10\beta {\log }_{10}\frac{l_i^t}{l_0}+c_i^t\right)}{2{\sigma }_i^2}$

• • where x^t represents the position of the target node at the moment t; β represents the path loss factor; L_i^t represents the power of the target node received by the i^th anchor node at the moment t; L₀=L_s^t−L(l₀); L_s^t, represents the transmit power of the target node at the moment t; L(l₀) represents the loss value when the reference distance is l₀; l_i^t=∥x^t−a_i^t∥; c_i^t=α_fl_i^t; a_i^t represents the position of the i^th anchor node at the moment t; and α_f represents the absorption coefficient.

Further, the optimization function by the joint estimation of the path loss factor and the target position is expressed as a following adaptive non-linear sampling (ANLS) problem:

$F\left(\psi \right)=\underset{\mathrm{\psi \dots 0}}{\arg \min }{A\psi -B}^2$ $\mathrm{where}A={\left[\begin{array}{ccc}1& -2{\left(a_1^t\right)}^T& {\mu }_1\ln {\mu }_1\\ ⋮& ⋮& ⋮\\ 1& -2{\left(a_N^t\right)}^T& {\mu }_N\ln {\mu }_N\end{array}\right]}^T,B=\left[\begin{array}{c}{\mu }_1-{a_1^t}^2\\ ⋮\\ {\mu }_N-{a_N^t}^2\end{array}\right];\psi ={\left[{x^t}^2,{\left(x^T\right)}^T,\vartheta \right]}^T;$

a_i^t represents the position of the i^th anchor node at the moment t; x^t represents the position of the target node at the moment

$t;\vartheta =\frac{\beta -\tilde{\beta }}{\tilde{\beta }};$

β represents the path loss factor; {tilde over (β)} represents an estimated parameter;

$\beta =\tilde{\beta }\left(1+\frac{\beta -\tilde{\beta }}{\tilde{\beta }}\right);$

and μ_i=10^{(L0−Lit+α6f)/5 {tilde over (β)}}.

Further, specific steps of re-expressing the optimization function by using the penalty function and finding the approximate solution by using the interior point method include:

• • according to the interior point method, introducing the penalty function to re-describe the problem:

$\Phi \left(\psi ,\theta \right)=\min \left({A\psi -B}^2-\theta \sum _{i=1}^N\ln \left({\psi }_i\right)\right)$

• • where θ>0 represents a penalty factor, θ gradually decreases to zero with iterations, and a feasible region is D={ψ_i>0,i=1, 2, . . . , k};
  • solving Karush-Kuhn-Tucker (KKT) by using Newton iteration method to find the optimal solution of the problem.

Further, specific steps of solving KKT by Newton iteration method to find the optimal solution of the problem include:

• • making z_k=(ψ_k)⁻¹, Z=diag (z₁, z₂ . . . , z_k), then expressing KKT conditions as:

${\psi }^T\left(A^{T}A\psi -A^{T}B\right)-z=0,\Upsilon \mathrm{ZE}-\theta E=0$

• • writing the KKT conditions in a matrix form:

$F\left(\psi ,z\right)=\left[\begin{array}{c}{\psi }^T\left(A^{T}A\psi -A^{T}B\right)-z\\ \Upsilon \mathrm{ZE}-\theta E\end{array}\right]=0$

• • according to F(ψ_k,z_k)+F′(ψ_k,z_k)(d_ψ^k,d_z^k)=0, building a linear system:

$\left[\begin{array}{cc}\nabla {\Phi }_{\mathrm{yy}}^2\left(\psi ,0\right)& -I\\ Z_k& {\Upsilon }_k\end{array}\right]\left[\begin{array}{c}{d}_k^{\psi }\\ d_k^z\end{array}\right]=-\left[\begin{array}{c}\nabla f\left(y_k\right)-z_k\\ {\Upsilon }_k{Z}_{k}E-{\theta }_{j}E\end{array}\right],$

• • where d_k^ψ and d_k^z represent Newton iteration directions;
  • finding an initial point that satisfying the system, and searching linearly along d_k^ψ and d_k^z respectively to find new iterative values; when tolerances of all conditions are established, a final solution ψ* is obtained.

Further, specific steps of transforming the original positioning problem into the generalized trust region sub-problem include:

$F\left(\psi \right)=\arg {\min }_{\mathrm{\psi \dots 0}}{{\stackrel{\_}{\omega }}_i\left(A\psi -B\right)}^{2}s.t.{\left(\psi \right)}^{T}D\psi +2f^T\psi =0\mathrm{where}{\stackrel{\_}{\omega }}_i=1-l_i^t/\sum _{i=1}^N{l}_i^t,D=\left[\begin{array}{ccc}0& 0_{1\times 4}& 0\\ 0& I_4& 0\\ 0& 0_{I\times 4}& 0\end{array}\right],f=\left[\begin{array}{c}-\frac{1}{2}\\ 0_{4\times 1}\end{array}\right];$

I and 0 represent respectively an identity matrix and a zero matrix;

• • introducing a regularization function h(ω_i):

$G\left(\psi ,\stackrel{\_}{\omega }\right)=\arg {\min }_{\mathrm{\psi \dots }-}{{\stackrel{\_}{\omega }}_i\left(A\psi -B\right)}^2+h\left({\stackrel{\_}{\omega }}_i\right)s.t.{\left(\psi \right)}^{T}D\psi +2f^t\psi =0h\left({\stackrel{\_}{\omega }}_i\right)=\sum _{i=1}^N{\delta }^2{\stackrel{\_}{\omega }}_i+\frac{2–\rho }{\rho }{\stackrel{\_}{\omega }}_i^{-\frac{\rho }{2-\rho }}.$

Further, transforming the original positioning problem into two sub-problems:

$\{\begin{array}{c}{\psi }^{k+1}=\underset{\psi }{\arg \min }G\left(\psi ,{\stackrel{\_}{\omega }}^k\right)\\ s.t.{\left(\psi \right)}^{T}D\psi +2f^t\psi =0\end{array}\{\begin{array}{c}{\stackrel{\_}{\omega }}^{k+1}=\underset{\stackrel{\_}{\omega }}{\arg \min }G\left({\psi }^{k+1},\stackrel{\_}{\omega }\right)\\ s.t.{\left(\psi \right)}^{T}D\psi +2f^t\psi =0\end{array}$

• • where k represents a number of the iterations;
  • a first sub-problem is a non-convex function, solved by constructing a block coordinate updating equation based on Lipschitz constant; and
  • a second sub-problem is a strictly convex function, solved in each iteration according to a following formula:

${{\stackrel{\_}{\omega }}_i^k={\left(A\psi -B\right)}^2+{\delta }^2)}^{\frac{\rho -2}{\rho }}.$

Further, specific steps of solving the first sub-problem by constructing the block coordinate updating equation based on Lipschitz constant include.

${\psi }^{k+1}=\arg {\min }_{\psi }\left\{{\nabla }_{\psi }G\left(\hat{\psi },{\stackrel{\_}{\omega }}^k\right),{\stackrel{\_}{\omega }}^{k-1},\psi -{\hat{\psi }}^k\right\}+{\lambda }^k{\psi -{\hat{\psi }}^k}^{2}s.t.{\left(\psi \right)}^{T}D\psi +2f^t\psi =0\mathrm{where}{\hat{\psi }}^k={\psi }^{k-1}{+}^k\left({\psi }^{k-1}-{\psi }^{k-2}\right),{\eta }^k=\frac{1}{12}\sqrt{{\lambda }^{k-1}/{\lambda }^k}$

is a given extrapolation coefficient, λ⁰=2∥k^Tωk∥;

• • introducing a multiplier ∈□, there is a global minimum when the KKT condition is obeyed

$\left({\lambda }^k{I}_3+D\right){\psi }^k=-{\kappa }^T{\stackrel{\_}{\omega }}^{k-1}\left(A{\hat{\psi }}^k-B\right)+{\lambda }^k{\hat{\psi }}^k-f,{\left(\psi \right)}^{T}D\psi +2f^t\mathrm{\psi 3}’0,\dots \max \left(-l^k,{\left(\sqrt{\left(A^{T}A\right)D^2\left(A^{T}A\right)}\right)}^{-1}\right)$

• • obtaining the optimal solution by a following formula:

$\psi \left(\right)={\left({\lambda }^k{I}_3+{*}^{}D\right)}^{-1}\times \left(-A^T{\stackrel{\_}{\omega }}^{k-1}\left(A{\hat{\psi }}^k-B\right)+{\lambda }^k{\hat{\psi }}^k-{*}^{}f\right)$

• • where obtaining , * by following formulas:

${\stackrel{\circ }{a}}^{}=\max \left[-\mathrm{diag}\left({\lambda }^k\right)/\mathrm{diag}\left(D\right),\right],{\psi \left(\right)}^{T}D\psi \left(\right)+2f^T\psi \left(\right)=0$

• • calculating the precise position coordinate x^t of the marine node and the estimated value β={tilde over (β)}(1+ϑ) of the path loss factor after a variable ψ is obtained, where ϑ is the first column of the last row of ψ.

Further, in order to further optimize the estimated value of the path loss factor, calculating an average value of the path loss factor according to the solved position information:

$\stackrel{\_}{\beta }=\frac{\sum _{i=1}^N\left(L_0-L_i-{\alpha }_f{\hat{l}}_i^t\right)}{\sum _{i=1}^N\left(10{\log }_{10}{\hat{l}}_i^t\right)}$

• • where {circumflex over (l)}_i^t represents an estimated value of a distance between an estimated position of an unknown node and the anchor node.

Compared with the prior art, the application has following beneficial effects:

• • the positioning method provided by the application is suitable for a highly dynamic marine environment, and may jointly estimate the node positions and the path loss factor. The real-time movement of all nodes, the path loss and the absorption effect of the underwater communication are considered, which is in line with the real situation when using wireless sensor networks to monitor the marine environment. A positioning method based on the interior point method and the block coordinate update is proposed to jointly estimate the position of marine nodes and the path loss factor, and the positioning with high-precision and high-robustness for ocean target nodes is realized.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a model diagram of a system for marine environmental monitoring wireless sensor networks (MEM-WSNs).

FIG. 2 is a flowchart of a node positioning method for marine environmental monitoring wireless sensor network node according to the present application.

FIG. 3 is a schematic diagram of a motion restriction area for the anchor node.

FIG. 4 is a comparison chart showing estimation errors of positions and of path loss factors corresponding to different anchor node numbers.

FIG. 5 is a comparison chart showing estimation errors of positions and of path loss factors corresponding to different deployment ranges.

FIG. 6 is a comparison chart showing estimation errors of positions and of path loss factors corresponding to different absorption coefficients.

DETAILED DESCRIPTION OF THE EMBODIMENTS

The present application is described in detail with the drawings and specific embodiments. This embodiment is implemented on the premise of the technical scheme of the present application, and the detailed implementation and specific operation process are given, but the protection scope of the present application is not limited to the following embodiments.

Embodiment 1

In order to solve the problem of low positioning accuracy caused by the highly dynamic underwater environment and the underwater absorption effect of acoustic communication, a high-precision and high-robustness positioning method for marine environmental monitoring wireless sensor networks (MEM-WSNs) is provided in the application. This method is suitable for the accurate position of underwater three-dimensional spatial nodes with an unknown path loss factor, and effectively improves the positioning accuracy of MEM-WSNs nodes.

The application is suitable for underwater node positioning requirements in the marine ecological environment, and the scenario is shown in FIG. 1.

As shown in FIG. 2, on the basis of comprehensively considering the complexity and the particularity of underwater target nodes for marine environmental monitoring, and aiming at the characteristics of node positioning of underwater wireless sensor networks, the application discloses a node positioning method for MEM-WSNs, which introduces an estimated parameter to jointly estimate the node position and the environmental path loss factor. The method includes two stages. Firstly, an original non-convex problem is transformed into a non-negative constrained least squares framework, and a two-stage method based on an interior point method (IPM) and a block coordinate update (BCU) is proposed to find the optimal solution of the marine node position. In the first stage, the penalty function is used to re-express the problem and find the approximate solution according to the IPM. In the second stage, the original problem is transformed into a generalized trust region sub-problem (GTRS), and the approximate solution obtained by the IPM is used as the initial estimation in the second stage, then, by combining with the BCU, the accurate estimation values of the marine node position and of the path loss factor are obtained. The node positioning method for MEM-WSNs proposed by the application includes following steps:

• • S1, establishing a received signal strength model according to a signal propagation loss and an underwater sound propagation absorption effect;
  • S2, constructing an optimization function by a joint estimation of a path loss factor and a target position through approximate operations such as Taylor expansion;
  • S3, re-expressing the optimization function by using the penalty function, and finding an approximate solution by using the IPM;
  • S4, transforming the original problem into the GTRS sub-problem; and
  • S5, taking the approximate solution found by the IPM as the initial estimation value of the GTRS sub-problem, and combining with the BCU to find the optimal solution.

In this embodiment, the S1 specifically includes:

In underwater target node positioning, the node positioning method based on RSS has attracted wide attention because of its low cost and no need for time synchronization. However, in the underwater environment, due to high dynamic and complexity of the underwater environment, the path loss factor easily changes with the change of environment, which leads to the decline of positioning accuracy. Due to the influence of water flow, the position of the sensor node and the target node of the buoy changes at every moment, and the movement of the buoy is usually restricted by the anchor chain and confined to a circle with a radius of M=√{square root over (L²−Dep²)}, where L is the length of the anchor chain and Dep is the water depth, as shown in FIG. 3. In this part, the received signal strength model is constructed according to the signal propagation loss and the underwater sound propagation absorption effect, and then the optimization function by the joint estimation is constructed according to the model.

Deploying a total of N anchor nodes and a target node in underwater and surface environments, assuming that the position of the i^th anchor node at the moment t is a_i^t=[a_i1^t,a_i2^t, a_i3^t]^T, T represents a transposition: the position of the target node at the moment t is x^t=[x₁^t, x₂^t,x₃^t]^T.

• • according to a wireless signal loss propagation model, a signal strength of information sent by a node received by node i is expressed as:

$\begin{array}{cc}{L}_i^t=L_s^t-L\left(l_0\right)-10\beta {\log }_{10}\frac{x^t-a_i^t}{l_0}-{\alpha }_fx^t-a_i^t+n_i& \left(1\right)\end{array}$

• • where L_i^t, represents a power of a target node received by an i^th anchor node at a moment t; L_s^t represents a transmit power of the target node at the moment t; L(l₀) represents a loss value when a reference distance is l₀, l₀ is usually 1 m; β represents the path loss factor: ∥□| represents a second-order norm; n_i represents a signal attenuation noise of the i^th anchor node, it is assumed that the noise variance of each anchor node is equal, if it obeys a Gaussian distribution with a mean of zero and a variance of σ_i², it may be expressed as n_i˜G(0,σ_i²): α_f represents an absorption coefficient, obtained by using Thorpe formula through the emission frequency f, that is:

$\begin{array}{cc}{\alpha }_f=0.11\frac{f^2}{1+f^2}+44\frac{f^2}{4100{+}^2}+2.75\times {10}^{-4}{f}^2+0.003& \left(2\right)\end{array}$

• • if the observation vector at the moment t is L^t=[L_i^t]^T, and the path loss factor is unknown, its joint probability density function may be expressed as:

$\begin{array}{cc}p\left(L^t❘x^t,\beta \right)=\prod _{i=1}^N\frac{1}{\sqrt{2{\mathrm{\pi \sigma }}_i^2}}\exp \left\{\frac{L_i-L_0+10\beta {\log }_{10}\frac{l_i^t}{l_0}+c_i^t}{2{\sigma }_i^2}\right\}& \left(3\right)\end{array}$

• • where l_i^t=∥x^t−a_i^t∥, L₀=L_s^t−L(l₀), c_i^t=α_fl_i^t.
  • obtaining the ML estimation by maximizing the probability density function, that is

$\begin{array}{cc}F\left({\hat{x}}^t,\beta \right)=\underset{x^t,\beta }{\arg \min }\sum _{i=1}^M\frac{\left(L_i^t-L_0+10\beta {\log }_{10}\frac{l_i^t}{l_0}+c_i^t\right)}{2{\sigma }_i^2}& \left(4\right)\end{array}$

The above formula is very difficult to solve due to the unknown path loss factor and the non-convexity of the height.

The S2 specifically includes:

• • assuming that the noise n_i is small enough, making l₀=1 m; obtaining the following approximate expression according to formula (1):

$\begin{array}{cc}{l}_i^t\approx {10}^{\frac{L_0-L_i^t+{\alpha }_f}{10\beta }}{10}^{-\frac{{\alpha }_f{l}_i^t}{10\beta }}& \left(5\right)\end{array}$

• • making λ=10^{−(αflit)/10β}, and it may be proved that when the absorption coefficient α_f is small enough, especially in deep water, the absorption term

${\alpha }_f{l}_i^t\frac{10\beta }{\ln 10},$

therefore 0<∥v∥□ 1.

Therefore, for a small absorption coefficient α_f, Taylor first-order expansion may be used to approximate λ.

$\begin{array}{cc}{10}^{-\frac{{\alpha }_f{l}_i^t}{10\beta }}\approx 1-\ln 10·\frac{{\alpha }_f{l}_i^t}{10\beta }& \left(6\right)\end{array}$

• • introducing formula (6) into formula (5), obtaining

$\begin{array}{cc}{10}^{\frac{L_0}{10\beta }}=l_i^t\left({10}^{\frac{l_i^t-{\alpha }_f}{10\beta }}+{10}^{\frac{L_0}{10\beta }}·\ln 10·\frac{{\alpha }_f}{10\beta }\right)& \left(7\right)\end{array}$

• • making

$v={10}^{\frac{L_0}{10\beta }},{\tilde{n}}_i={10}^{\frac{L_i^t-{\alpha }_f}{10\beta }},$

then

$\begin{array}{cc}v=l_i^t\left({\tilde{n}}_i+v·\frac{\ln 10{\alpha }_f}{10\beta }\right)& \left(8\right)\end{array}$

• • obtaining a distance l_i^t;

$\begin{array}{cc}{l}_i^t=\frac{v}{\left({\tilde{n}}_i+v·\frac{\ln 10{\alpha }_f}{10\beta }\right)}& \left(9\right)\end{array}$

• • transforming formula (9)

$\begin{array}{cc}{l}_i^t\approx \frac{\frac{10\beta }{{\alpha }_f\ln 10}\left({\tilde{n}}_i+v·\frac{\ln 10{\alpha }_f}{10\beta }\right)-\frac{10\beta }{{\alpha }_f\ln 10}{\tilde{n}}_i}{\left({\tilde{n}}_i+v·\frac{\ln 10{\alpha }_f}{10\beta }\right)}& \left(10\right)\end{array}$

• • rewriting formula (10) as

$\begin{array}{cc}{l}_i^t\approx \frac{10\beta }{{\alpha }_f\ln 10}-\frac{\frac{10\beta }{{\alpha }_f\ln 10}{\tilde{n}}_i}{\left({\tilde{n}}_i+v·\frac{\ln 10{\alpha }_f}{10\beta }\right)}& \left(11\right)\end{array}$

• • making

$\zeta =\frac{\frac{10\beta }{{\alpha }_f\ln 10}{\tilde{n}}_i}{\left({\tilde{n}}_i+v·\frac{\ln 10{\alpha }_f}{10\beta }\right)},$

dividing the numerator and the denominator by ñ_i, at the same time

$\begin{array}{cc}\zeta =\frac{\frac{10\beta }{{\alpha }_f\ln 10}}{\left(1+\frac{v}{{\tilde{n}}_i}·\frac{\ln 10{\alpha }_f}{10\beta }\right)}& \left(12\right)\end{array}$

• • making

$ϛ=\frac{v_i}{{\tilde{n}}_i}·\frac{\ln 10{\alpha }_f}{10\beta }$

may prove |ζ|□ 1 in most underwater wireless sensor making network environments. Therefore, when the absorption coefficient α_f, is small enough, the first-order Taylor expansion with ζ=0 may be used to approximate formula (12) as follows:

$\begin{array}{cc}\frac{1}{1+\frac{v}{{\tilde{n}}_i}·\frac{\ln 10{\alpha }_f}{10\beta }}\approx 1-\frac{v}{{\tilde{n}}_i}·\frac{\ln 10{\alpha }_f}{10\beta }& \left(13\right)\end{array}$

• • introducing formula (13) into formula (11), obtaining

$\begin{array}{cc}{\left(l_i^t\right)}^2\approx {\left(\frac{v}{{\tilde{n}}_i}\right)}^2={\left({10}^{\frac{L_0-L_i^t+{\alpha }_f}{10\beta }}\right)}^2& \left(14\right)\end{array}$

• • however, because the path loss factor β is unknown, its distance cannot be obtained. Making

$\beta =\tilde{\beta }\left(1+\frac{\beta -\tilde{\beta }}{\tilde{\beta }}\right)\mathrm{and}\vartheta =\frac{\beta -\tilde{\beta }}{\tilde{\beta }}$

by introducing the estimated parameter {tilde over (β)}, where {tilde over (β)}≠0. And formula (14) may be transformed into

$\begin{array}{cc}{\left(l_i^t\right)}^2\approx {10}^{\frac{L_0-L_i^t+{\alpha }_f}{5\tilde{\beta }\left(1+\vartheta \right)}}& \left(15\right)\end{array}$

• • assuming that |ϑ| is relatively small, that is

$\frac{1}{1+\vartheta }\approx 1-\vartheta ,$

then

$\begin{array}{cc}{\left(l_i^t\right)}^2\approx {\mu }_i^{l-\vartheta }\mathrm{where}{\mu }_i={10}^{\frac{L_0-L_i^t+{\alpha }_f}{5\tilde{\beta }}}.& \left(16\right)\end{array}$

Further approximating formula (16) by performing the first-order Taylor expansion on μ_t^−ϑ

$\begin{array}{cc}{\left(l_i^t\right)}^2\approx {\mu }_i\left(1-\vartheta \ln {\mu }_i\right)& \left(17\right)\end{array}$

• • therefore, the LS framework of the distance is as follows:

$\begin{array}{cc}\underset{x^t,\vartheta }{\arg \min }\sum _{i=1}^N{\left[{x^t-a_i^t}^2-{\left(l_i^t\right)}^2\right]}^2& \left(18\right)\end{array}$

• • making ψ=[∥x^t∥², (x^t)^T, ϑ]^T, expanding the formula (18), and then transforming into the adaptive non-linear sampling (ANLS) problem:

$\begin{array}{cc}F\left(\psi \right)=\underset{\mathrm{\psi \dots 0}}{\arg \min }{A\psi -B}^2& \left(19\right)\end{array}$

• • where:

$\begin{array}{cc}A={\left[\begin{array}{ccc}1& -2\left(a_1^t\right)& {\mu }_1\ln {\mu }_1\\ ⋮& ⋮& ⋮\\ 1& -2{\left(a_N^t\right)}^T& {\mu }_N\ln {\mu }_N\end{array}\right]}^T,B=\left[\begin{array}{c}{\mu }_1-{a_1^t}^2\\ ⋮\\ {\mu }_N-{a_N^t}^2\end{array}\right].& \left(20\right)\end{array}$

The S3 specifically includes:

• • according to the IPM, introducing the penalty function to re-express the problem (19)

$\begin{array}{cc}\Phi \left(\psi ,\theta \right)=\min \left({A\psi -B}^2-\theta \sum _{i=1}^N\ln \left({\psi }_i\right)\right)& \left(21\right)\end{array}$

• • where θ>0 represents the penalty factor, θ gradually decreases to zero with an iteration, and the feasible region is D={ψ_i>0,i=1, 2, . . . , k}. The first-order gradient and the second-order gradient of formula (21) are respectively

$\begin{array}{cc}\nabla {\Phi }_y\left(\psi ,\theta \right)=2A^T\left(A\psi -B\right),-{\theta }_j·Y^{-1}·E,\nabla {\Phi }_{\mathrm{yy}}^2\left(\psi ,\theta \right)=2A^{T}A-{\theta }_k·Y^{-2}·E& \left(22\right)\end{array}$

• • where E=[1 1 . . . 1]^T, Υ=diag(ψ₁, ψ₂ . . . ψ_k), the problem may be solved by using Newton iteration method to find the optimal solution: making z_k=(ψ_k)⁻¹, Z=diag(z₁, z₂ . . . , z_k) expressing KKT conditions as:

$\begin{array}{cc}{\psi }^T\left(A^{T}A\psi -A^{T}B\right)-z=0,\mathrm{YZE}-\theta E=0& \left(23\right)\end{array}$

• • writing the KKT conditions in a matrix form:

$\begin{array}{cc}F\left(\psi ,z\right)=\left[\begin{array}{c}{\psi }^T\left(A^{T}A\psi -A^{T}B\right)-z\\ \mathrm{YZE}-\theta E\end{array}\right]=0& \left(24\right)\end{array}$

• • according to F(ψ_k, z_k)+F′(ψ_k, z_k) (d_ψ^k,d_z^k)=0, building a linear system:

$\begin{array}{cc}\left[\begin{array}{cc}\nabla {\Phi }_{\mathrm{yy}}^2\left(\psi ,\theta \right)& -I\\ Z_k& {{\rm Y}}_k\end{array}\right]\left[\begin{array}{c}{d}_k^{\psi }\\ d_k^z\end{array}\right]=-\left[\begin{array}{c}\nabla f\left(y_k\right)-z_k\\ {{\rm Y}}_k{Z}_{k}E-{\theta }_{j}E\end{array}\right],& \left(25\right)\end{array}$

• • where d_k^ψ and d_k^z represent Newton iteration directions.

A great advantage of IP-based method is that it does not need a strictly feasible initial point, but only needs non-negative conditions. Once the initial point satisfying the above system is found, a linear search may be conducted along d_k^ψ to find a new iteration value, which follows ψ_k+1=ψ_k+α_kd_k^ψ, where α_k is the step size. Similarly, z_k+1=z_k+α_kd_k^z. When tolerances of all conditions are established, a final solution aψ* may be obtained.

The S4 specifically includes:

$\begin{array}{cc}F\left(\psi \right)=\underset{\mathrm{\psi \dots 0}}{\arg \min }{{\varpi }_j\left(A\psi -B\right)}^2& \left(26\right)\end{array}$ $s.t.{\left(\psi \right)}^{T}D\psi +2f^T\psi =0$ $\mathrm{where}$ ${\varpi }_i=1-l_i^t/\sum _{i=1}^N{l}_i^t,$ $D=\left[\begin{array}{ccc}0& 0_{1\times 4}& 0\\ 0& I_4& 0\\ 0& 0_{1\times 4}& 0\end{array}\right],$ $f=\left[\begin{array}{c}-\frac{1}{2}\\ 0_{4\times 1}\end{array}\right];$

I and 0 represent respectively an identity matrix and a zero matrix.

Nodes in underwater wireless sensor networks are vulnerable to Byzantine attacks and become fault nodes, which leads to the problem in (26) not being applicable to the actual situation, so a regularization function h(ω_i) is introduced, that is

$\begin{array}{cc}G\left(\psi ,\varpi \right)=\underset{\mathrm{\psi \dots 0}}{\arg \min }{{\varpi }_i\left(A\psi -B\right)}^2+h\left({\varpi }_i\right)& \left(27\right)\end{array}$ $s.t.{\left(\psi \right)}^{T}D\psi +2f^t\psi =0$

• • where:

$\begin{array}{cc}h\left({\varpi }_i\right)=\sum _{i=1}^N{\delta }^2{\varpi }_i+\frac{2-\rho }{\rho }{\varpi }_i^{-\frac{\rho }{2-\rho }}.& \left(28\right)\end{array}$

The S5 specifically includes:

• • transforming the original problem into two sub-problems, that is

$\begin{array}{cc}\{\begin{array}{c}{\psi }^{k+1}=\underset{\psi }{\arg \min }G\left(\psi ,{\varpi }^k\right)\\ s.t.{\left(\psi \right)}^{T}D\psi +2f^t\psi =0\end{array}& \left(29\right)\end{array}$ $\{\begin{array}{c}{\varpi }^{k+1}=\underset{\varpi }{\arg \min }G\left({\psi }^{k+1},\varpi \right)\\ s.t.{\left(\psi \right)}^{T}D\psi +2f^t\psi =0\end{array}$

• • where k represents a number of the iterations.

The second problem in formula (29) is a strictly convex function, so it may be solved according to formula (30) in each iteration

$\begin{array}{cc}{\varpi }_i^k={\left({\left(A\psi -B\right)}^2+{\delta }^2\right)}^{\frac{\rho -2}{\rho }}.& \left(30\right)\end{array}$

The first problem in formula (29) is non-convex, so Lipschitz constant λ is introduced to construct BCU equation.

$\begin{array}{cc}{\psi }^{k+1}=\arg {\min }_{\psi }\left\{{\nabla }_{\psi }G\left(\hat{\psi },{\varpi }^k\right),{\varpi }^{k-1},\psi -{\widehat{\psi }}^k\right\}+{\lambda }^k{\psi -{\widehat{\psi }}^k}^2& \left(31\right)\end{array}$ $s.t.{\left(\psi \right)}^{T}D\psi +2f^t\psi =0$ $\mathrm{where}$ $\begin{array}{cc}{\widehat{\psi }}^k={\psi }^{k-1}+{\eta }^k\left({\psi }^{k-1}-{\psi }^{k-2}\right)& \left(32\right)\end{array}$ ${\eta }^k=\frac{1}{12}\sqrt{{\lambda }^{k-1}/{\lambda }^k}$

is a given extrapolation coefficient, λ⁰=2∥k^Tωk∥.

According to the definition of the Lipschitz constant

$\begin{array}{cc}{\hat{\psi }}^k={\psi }^{k-1}+{\eta }^k\left({\psi }^{k-1}-{\psi }^{k-2}\right)& \left(33\right)\end{array}$

• • literature shows that this method is capable of converging.

According to Theorem 1, the first problem in formula (31) has a global minimum in each iteration. In the process of solving, a multiplier ∈ is introduced, and there is a global minimum when the KKT condition is obeyed

$\begin{array}{cc}\left({\lambda }^k{I}_3+\wp D\right){\psi }^k=-{\kappa }^T{\varpi }^{k-1}\left(A{\hat{\psi }}^k-B\right)+{\lambda }^k{\hat{\psi }}^k-\wp f,& \left(34\right)\end{array}$ ${\left(\psi \right)}^{T}D\psi +2f^t\psi =0,$ $\wp \dots \max \left(-l^k,{\left(\sqrt{\left(A^{T}A\right)D^2\left(A^{T}A\right)}\right)}^{-1}\right)$

Obtaining the optimal solution by formula (35)

$\begin{array}{cc}\psi \left(\wp \right)={\left({\lambda }^k{I}_3+{\wp }^{*}D\right)}^{-1}\times (-A^T{\varpi }^{k-1}\left(A{\hat{\psi }}^k-B\right)+{\lambda }^k{\hat{\psi }}^k-{\wp }^{*}f& \left(35\right)\end{array}$

• • where obtaining , * by formula (36):

$\begin{array}{cc}{\wp }^{å}=\max \left[-\mathrm{diag}\left({\lambda }^k\right)/\mathrm{diag}\left(D\right),\wp \right],& \left(36\right)\end{array}$ ${\psi \left(\wp \right)}^{T}D\psi \left(\wp \right)+2f^T\psi \left(\wp \right)=0$

• • calculating the precise position coordinate x^t of the marine node and the estimated value β={tilde over (β)}(1+ϑ) of the path loss factor after a variable y is obtained, where ϑ is the first column of the last row of ψ, in order to further optimize the estimated value of the path loss factor, calculating an average value of the path loss factor according to the solved position information, that is:

$\begin{array}{cc}\overline{\beta }=\frac{\sum _{i=1}^N\left(L_0-L_i-{\alpha }_f{\hat{l}}_i^t\right)}{\sum _{i=1}^N\left(10{\log }_{10}{\hat{l}}_i^t\right)}& \left(37\right)\end{array}$

• • where {circumflex over (l)}_i^t, represents an estimated value of a distance between an estimated position of an unknown node and the anchor node.

Next the problem model established by the application and the node positioning method for MEM-WSNs are simulated and analyzed through MATLAB R2022a, so as to further introduce the implementation of the example of the application and its advantages in the field of marine environmental monitoring. According to different scenarios, compared with the existing methods WLS, SR-LS, SR-IRLS, PM, RNLA, the minimum root mean square error is used as the evaluation index to evaluate the performance of the proposed method, that is

$\begin{array}{cc}{\mathrm{RMSE}}_x=\sqrt{\frac{1}{\mathrm{MC}}\sum _{i=1}^{\mathrm{MC}}{\left(\widehat{x}-x\right)}^2},& \left(38\right)\end{array}$ ${\mathrm{RMSE}}_{\beta }=\sqrt{\frac{1}{\mathrm{MC}}\sum _{i=1}^{\mathrm{MC}}{\left(\widehat{\beta }-\beta \right)}^2}$

• • where MC is a total number of Monte Carlo experiments, {circumflex over (x)} and {circumflex over (β)} are respectively the estimated values in the i^th MC iteration. In order to simulate the highly dynamic of the underwater real environment, the positions of the anchor nodes and the target node are randomly changed in each simulation. The six positioning methods are compared from following three aspects:
  • 1) different anchor node numbers;
  • 2) different deployment ranges;
  • 3) different absorption coefficients.

FIG. 4 shows the estimation errors of the positions and of the path loss factors corresponding to different anchor node numbers. With the increase of the anchor node number, the information of RSS used for estimation increases, so the estimation error of each method decreases with the increase of the anchor node number. As can be seen from the figure, with the increase of the anchor node number, except SR-LS and WLS, the positioning errors of other methods decrease. For IPM-BCU, WLS and SR-LS, they are more sensitive to the change of the anchor node number, in the process of increasing the anchor node number N from 6 to 15, their positioning accuracies increase by nearly 60%, but of IPM-BCU is higher. However, the positioning accuracies of SR-IRLS, IPM and ANLS are improved by 22%, 20% and 14% respectively. When the number of nodes is low, IPM-BCU outperforms SR-LS and WLS and is close to IPM, RNLA and SR-IRLS. For the estimation error of the path loss factor, the estimation errors of IPM-BCU and ANLS are similar and obviously better than that of WLS and SRLS. Therefore, on the whole, the estimation effect of IPM-BCU is obviously better than other methods.

FIG. 5 shows the estimation errors of the positions and of the path loss factors corresponding to different deployment ranges. Due to the increase of deployment range, the influence of absorption effect on positioning accuracy is gradually increasing, and the position estimation errors of each method are increasing. In the process of increasing the side length of the deployment area from 50 m to 120 m, the estimation errors of SR-LS, WLS, IPM-BCU, IPM, RNLA and SR-IRLS increase by 500%, 500%, 216%, 202%, 200% and 147% respectively. Compared with other positioning methods, the positioning accuracies are still improved by 10%-25%. Although the error increase ratio of IPM-BCU is large, it may be seen from the figure that the error of IPM-BCU is still the smallest. For the estimation error of the path loss factor, when the side length of the deployment area is 60 m, the estimation error of SRLS and WLS decrease, but with the increase of the deployment area, the overall error trend still increases. The final results show that the IPM-BCU method has higher estimation accuracy for the position and the path loss factor, and has greater advantages compared with other methods.

FIG. 6 shows the estimation errors of the positions and the path loss factors corresponding to different absorption coefficients. As can be seen from the figure, the position estimation errors of WLS and SRLS increase with the increase of the absorption coefficient, and the growth ratio reaches 178% and 190% respectively. The other methods are relatively stable, and the fluctuation range of position error is less than 10%, showing a strong robustness. At the same time, it is found that the proposed positioning method is the best for estimating the target position and the closest positioning scheme to CRLB. For a single IPM method, its performance is at a medium level. In the estimation of the path loss factor, when α_f=[0,0.12], the estimation accuracy of IPM-BCU is higher than all methods except SR-IRLS, and when α_f=[0.08,0.12], the performance of IPM-BCU is the best. When α_f=[0.12,0.2], the performance of IPM-BCU is still better than other methods except IPM and RNLA. As can be seen from the figure, the robustness of each method is not ideal for estimating the path loss factor, but the proposed method still maintains a high performance.

To sum up, the node positioning method for MEM-WSNs provided by the application realizes a positioning with high-precision and high-robustness for the marine target nodes, and has the best performance compared with the existing positioning methods in different scenarios.

The preferred embodiments of the present application have been described in detail above. It should be understood that those skilled in the art can make many modifications and changes according to the concept of the present application without creative work. Therefore, any technical scheme that may be obtained by a person skilled in the technical field through logical analysis, reasoning or limited experiments on the basis of the existing technology according to the concept of the present application should be within the protection scope determined by the claims.

Claims

1. A node positioning method for marine environmental monitoring wireless sensor networks comprising: establishing a received signal strength model according to a signal propagation loss and an underwater sound propagation absorption effect;establishing an optimization function by a joint estimation of a path loss factor and a target position;re-expressing the optimization function by using a penalty function, and finding an approximate solution of position coordinates of marine nodes by using an interior point method; andtransforming an original positioning problem into a generalized trust region sub-problem, using the approximate solution of the interior point method as an initial estimated value of the generalized trust region sub-problem, and finding an optimal solution by combining a block coordinate update.

2. The node positioning method for marine environmental monitoring wireless sensor networks according to claim 1, wherein establishing the received signal strength model comprises: according to a wireless signal loss propagation model, expressing a signal strength of information sent by a node received by node i as:

3. The node positioning method for marine environmental monitoring wireless sensor networks according to claim 1, wherein the joint estimation of the path loss factor and the target position is as follows:

4. The node positioning method for marine environmental monitoring wireless sensor networks according to claim 3, wherein the optimization function by the joint estimation of the path loss factor and the target position is expressed as a problem of adaptive non-linear sampling (ANLS):

5. The node positioning method for marine environmental monitoring wireless sensor networks according to claim 4, wherein re-expressing the optimization function by using the penalty function and finding the approximate solution by using the interior point method comprise: introducing the penalty function to re-describe the problem according to the interior point method:

6. The node positioning method for marine environmental monitoring wireless sensor networks according to claim 5, wherein solving KKT by using the Newton iteration method to find the optimal solution of the problem comprises: making zk=(ψk)−1, Z=diag(z1, z2, . . . , zk), then expressing KKT conditions as:

7. The node positioning method for marine environmental monitoring wireless sensor networks according to claim 4, wherein transforming the original positioning problem into the generalized trust region sub-problem comprises:

8. The node positioning method for marine environmental monitoring wireless sensor networks according to claim 5, wherein transforming the original positioning problem into two sub-problems comprises:

9. The node positioning method for marine environmental monitoring wireless sensor networks according to claim 8, wherein solving the first sub-problem by constructing the block coordinate updating equation based on Lipschitz constant comprises:

10. The node positioning method for marine environmental monitoring wireless sensor networks according to claim 9, wherein in order to further optimize the estimated value of the path loss factor, an average value of the path loss factor is calculated according to solved position information: