J. Cent. South Univ. (2016) 23: 1981−1989 DOI: 10.1007/s11771-016-3255-1
An underwater acoustic data compression method based on compressed sensing GUO Xiao-le(郭晓乐)1, 2, YANG Kun-de(杨坤德)1, 2, SHI Yang(史阳)1, 2, DUAN Rui(段睿)1, 2 1. School of Marine Science and Technology, Northwestern Polytechnical University, Xi’an 710072, China; 2. Key Laboratory of Ocean Acoustics and Sensing (Northwestern Polytechnical University), Ministry of Industry and Information Technology, Xi’an 710072, China © Central South University Press and Springer-Verlag Berlin Heidelberg 2016 Abstract: The use of underwater acoustic data has rapidly expanded with the application of multichannel, large-aperture underwater detection arrays. This study presents an underwater acoustic data compression method that is based on compressed sensing. Underwater acoustic signals are transformed into the sparse domain for data storage at a receiving terminal, and the improved orthogonal matching pursuit (IOMP) algorithm is used to reconstruct the original underwater acoustic signals at a data processing terminal. When an increase in sidelobe level occasionally causes a direction of arrival estimation error, the proposed compression method can achieve a 10 times stronger compression for narrowband signals and a 5 times stronger compression for wideband signals than the orthogonal matching pursuit (OMP) algorithm. The IOMP algorithm also reduces the computing time by about 20% more than the original OMP algorithm. The simulation and experimental results are discussed. Key words: compressed sensing; underwater acoustic signal; compression ratio; improved orthogonal matching pursuit (IOMP)
1 Introduction Several areas in signal processing, such as beam forming, array shape estimation, frequency estimation and tracking, detection and Doppler analysis are of particular relevance to the acoustic surveillance task. Those signal processing areas generally apply a large array that leads to the rapid expansion of underwater acoustic data and requires a large data transmission rate. Underwater acoustic communication shows a serious contradiction between the transmission rate and distance. The underwater acoustic signal frequency must be increased to improve the transmission rate. However, an increase in signal frequency will also increase the transmission loss and decrease the transmission distance. Data compression technology [1] is therefore indispensable for underwater acoustic signal transmission. A widely debated issue in the information field, data compression technology has been widely used in radar, communication systems, seismic exploration and many aspects of the national economy. Several compression methods can be used for voice and image, such as K_L transformation [2], DCT transformation [3], Huffman coding [4], and wavelet transformation [5]. Given that these methods take human auditory masking effect or
visual sensitivity as the basic condition of quantitative coding to improve compression performance, they cannot be directly applied in underwater acoustic signal compression [6]. According to the characteristics of underwater acoustic environment and signal, an underwater acoustic data compression method based on compressed sensing [7] is presented in this study. This compression method selects the effective frequency band for data storage at data receiving terminal. After data transmission, the improved orthogonal matching pursuit (IOMP) algorithm is used to reconstruct the original underwater acoustic data at a data processing terminal. The orthogonal matching pursuit (OMP) algorithm [8] is often used when the sparsity of a signal is known. However, sparsity cannot be easily measured in most cases, especially when a large quantity of data is used. Loop iteration cannot be easily stopped, and the requirements for computing speed are difficult to meet. The proposed IOMP algorithm differs from the OMP algorithm in which it synchronously processes several columns of the recovery matrix at each time and chooses the absolute value of the 2-norm difference between the original and the reconstructed signals as the stop condition of the loop iteration. Therefore, the IOMP algorithm is a robust process that can reduce the computing time. To demonstrate the feasibility of the proposed
Foundation item: Project(11174235) supported by the National Natural Science Foundation of China; Project(3102014JC02010301) supported by the Fundamental Research Funds for the Central Universities, China Received date: 2015−06−11; Accepted date: 2015−11−19 Corresponding author: YANG Kun-de, Professor, PhD; Tel: +86−29−88460586; E-mail:
[email protected]
J. Cent. South Univ. (2016) 23: 1981−1989
1982
method, this study uses a single frequency signal, a linear frequency modulation (LFM) signal, and a wideband random signal as simulation signals. Actual underwater acoustic signals are then used to evaluate the algorithm. The experimental data of a 12-element array are used to verify the compression method.
2 Theoretical model
i (k )
The basic theory underlying this study is categorized into underwater acoustic signal compression theory and underwater acoustic signal reconstruction algorithm. Compression theory uses a measurement matrix [7] to sparsely represent the original underwater acoustic signal. The reconstruction algorithm uses the IOMP algorithm to reconstruct the original underwater acoustic signal. 2.1 Underwater acoustic signal compression theory Consider a finite-length, one-dimensional (1D) underwater acoustic signal x, which is viewed as an N×1 column vector in RN with elements x[i] (i=1, 2, …, N). The high-dimensional data can be treated by vectorizing it into a long 1D vector. Any signal in RN can be represented in terms of N×1 vectors φi (i=1, 2, … , N). For simplicity, the basis vectors are assumed as orthonormal in this study. The N×N basis matrix Ψ is formed by stacking of the vectors φi as columns. The signal x can be expressed as follows: x i 1 ai i or x a N
band in the frequency domain. Frequency-domain analysis has become a major technique to analyze underwater acoustic signals. Therefore, the frequency domain is chosen in this study as the sparse domain for underwater acoustic signals. The frequency domain basis vectors φi (i=1, 2, …, N) are calculated as follows:
(1)
where a is the N×1 column vector of the weighting coefficients ai x , i iT x; the angle bracket denotes the inner product operator; the superscript T denotes the transposition; x and a are equivalent representations of the same underwater acoustic signal, with x in the time domain and a in the Ψ domain. The basis of compressed sensing is that the signal is sparse in the Ψ domain. Underwater acoustic signals are sparse in the frequency domain for several reasons: The underwater acoustic signals that are received by the detection equipment are mainly used to detect underwater targets. The detection equipment comprises active and passive sonars. The active sonar detects targets by emitting a sound pulse waveform and accepting the reflection echo of targets. The passive sonar mainly functions as a detector of the physical radiation noise of targets. Physical radiation noise is grouped into three categories. Category I is machinery noise which is emitted by mechanical equipment during its operation. Category II is propeller noise, which is emitted by a propeller during its rotation. Category III is hydrodynamic noise. All of the above-mentioned underwater acoustic signals have a particular frequency
N 1
2π
ei (k ) exp( j N nk ), k 1,
2, , N
(2)
n 0
where φi is the ith column of Ψ, and ei is the ith column of the unit diagonal matrix E. In Eq. (1), a is the sparse coefficient that has K (K≪N) elements that are far greater than zero, where K denotes the sparsity of the underwater acoustic signal. Underwater acoustic signals have little a priori information, so one cannot precisely determine the sparsity K, which indicates the number of points whose amplitude in the frequency domain is greater than zero. Sparsity is also related to the sampling frequency. The sparsity and the number of signal effective points increase along with the sampling frequency. This study sets a threshold to estimate K, and the corresponding threshold is deduced on the basis of the reconstructed effect. The threshold determines the compression ratio (yCR), which is calculated as follows:
yCR
N M
(3)
where N is the length of the original underwater acoustic signal and M is the measurement times or the number of saved data. A high threshold signifies a high yCR but a poor reconstructed effect. Consider the linear measurement process that computes M inner products between x and vector i (i=1, 2, …, N) as yi x , i . By stacking the measurements yi into the M×1 vector y and by stacking the measurement vectors jT as rows into the M×N matrix Φ and substituting them in Eq. (1), we obtain the following equation: y x a a
(4)
where Φ is an M×N matrix; Ψ is an N×N matrix; Θ is an M×N matrix; x is an N×1 vector; a is an N×1 vector; and y is an M×1 vector. The measurement process is nonadaptive, namely, Φ does not depend in any way on signal x. Underwater acoustic signals have a certain regularity or statistical properties, so most of these signals are irrelevant to the Gaussian matrix [9]. The Gaussian matrix may be related to marine environment noise. Therefore, this study chooses the Gaussian matrix as the measurement matrix. This study aims to determine the dimension of Φ. The measurement times M must be determined before
J. Cent. South Univ. (2016) 23: 1981−1989
1983
the original signal x is stably reconstructed. The reconstruction will not be possible if the measurement process damages the information in x. Unfortunately, the measurement process is a linear algebra problem and M≪N, so the unknown values outnumber the equations and hence make the solution ill-posed in general. In this case, the measurement vector y is just a linear combination of the K columns of Θ, whose projection is ai>>0. Therefore, if we know a priori in which K entries of a are bigger than zero, then we can form an M×K system of linear equations to solve these non-zero entries. In this way, the number of equations M now equals or exceeds the number of unknowns K. To ensure that the M×K system is well-conditioned, the following equation must be satisfied for any vector v sharing the same K non-zero entries:
1
v v
2
1
(5)
2
where ε>0 and || ||2 denotes the 2-norm in Eq. (5). The locations of the K non-zero entries in a cannot be determined in practice. Interestingly, a stable inverse for both K-sparse and compressible signals can be obtained when Θ satisfies Eq. (5) for an arbitrary 3K-sparse vector v. This condition is called restricted isometry property (RIP) [10−11]. Incoherence, a related condition, requires that the rows {ϕi} of Φ must not sparsely represent the columns { i} of Ψ. The Gaussian matrix Φ is incoherent with the unit diagonal matrix I. Ψ is computed with the use of Eq. (2). According to RIP, the reconstructed effect becomes perfect for the Gaussian matrix when the measurement times M≥cK×lg(N/K), with c being a small constant [11]. This study estimates the sparsity of underwater acoustic signals, and the number of measurements is calculated as follows: M 2 K lg( N K )
(6)
where Φ, y, and Θ can be calculated according to the M value. The Gaussian matrix, Bernoulli matrix [12], Fourier matrix and Hadamard matrix [13] are always used as measurement matrixes in compressed sensing theory. Different measurement matrixes require various measurement times to accurately reconstruct the original signal. 2.2 Underwater acoustic signal reconstruction algorithm The proposed IOMP algorithm chooses several columns that have the maximum correlation coefficient between the recovery matrix and the residual value and then advances these columns to the front of the recovery matrix. The recovery matrix is then rearranged by
analogy. In this way, the proposed algorithm can reduce the computing time. The IOMP algorithm chooses the absolute value of the 2-norm difference between the original and the reconstructed signals as the stop condition of the loop iteration. The loop iteration will stop when the difference is less than the threshold value ε. Therefore, the IOMP algorithm is robust and can reconstruct the original signal whose sparsity is unknown. The inputs of the IOMP algorithm during the reconstruction include the recovery matrix Θ, the measurement values y, and the estimated sparsity K, whereas the output is the reconstructed signal x. The correlation coefficient is calculated as follows:
u u j | u j rnew , j
j 1, 2, , N
(7)
where θj is the jth column of Θ; rnew is the new residual value obtained by the iteration; uj is the correlation coefficient between θj and rnew. The IOMP algorithm uses Eqs. (8) and (9) to approach the original signal in the frequency domain and to update the residual value, respectively. aˆ T ( T ) 1 y
(8)
rnew y aˆ
(9)
where xˆ rew is the recovery signal in each iteration; aˆ is the original signal in the frequency domain after each iteration; ΘΛ is the recovery matrix after the rearrangement in each iteration. The principle and implementation of the IOMP algorithm are outlined as follows: 1) The IOMP algorithm sets the initial residual value as r0=y, sets the number of choosing columns in each iteration as L=3, and sets the index value as Λ , J=. 2) The algorithm calculates the correlation coefficient u by using Eq. (7) and advances those columns that have the maximum correlation coefficient to the front of the recovery matrix. 3) The algorithm uses J to update the support set Λ. 4) The algorithm calculates xˆrew by using Eq. (8) and updates the residual value by using Eq. (9). 5) The algorithm will stop the iteration when the absolute value of the 2-norm difference between the original and the reconstructed signals is less than ε. Otherwise, the algorithm repeats all procedures from step 2. 6) After the iteration is stopped, the original signal is reconstructed by x α . Figure 1 shows the implementation procedures of the underwater acoustic data compression method based on compressed sensing. The principle and implementation of this
J. Cent. South Univ. (2016) 23: 1981−1989
1984
where x(i) is the original signal; xˆ(i ) is the reconstructed signal; f ori is the normalized frequency of the original signal; f rec is the normalized frequency of the reconstructed signal; X(f) is the frequency domain description of the original signal; and e(f) is the difference between the original and the reconstructed signals in the frequency domain.
3 Single hydrophone signal compression This study chooses the single frequency signal, the LFM signal, and the wideband random signal as the simulation signals and uses actual underwater acoustic signals to verify the compression method. Fig. 1 Process flow diagram of underwater acoustic compression method
compression method are outlined as follows: 1) A threshold is set to estimate the sparsity K of the original signal. 2) M=2Klg(N/K) is taken as the measurement times. 3) M is used to generate the measurement matrix Φ. Φ is then utilized to measure the original signal and to obtain the measurement values with the use of Eq. (4). Most of the underwater acoustic signals are sparse in the frequency domain. In this study, Ψ is calculated with Eq. (2), and the recovery matrix Θ is calculated by Θ=ΦΨ. 5) The IOMP algorithm is used to calculate unknown variables with the use of the least square method [14−16] a=ΘT(ΘΘT)−1y. 6) The original signal is reconstructed by x=Ψa. 2.3 Compression algorithm evaluation indexes The proposed compression method does not code for the underwater acoustic signal, so the method cannot be evaluated by the bits error ratio (BER). This study uses yCR, the cross correlation coefficient (R), the largest normalized frequency domain error (E), and the frequency domain quantization signal noise ratio y(FQSNR) as objective evaluation indexes. A favorable reconstructed effect is obtained when a large R, a small E, and a large yFQSNR are derived.
i 1 xˆ (i) x (i) N N i 1 x 2 (i) i 1 xˆ 2 (i) N
R
(10)
E max f ori f rec yFQSNR
10 lg
(11) N X 2( i 1
N 2 e ( i 1
f)
f)
(12)
3.1 Analog signal simulation The single frequency signal in this study is set as follows: x(t ) cos(4000πt )
(13)
The sample frequency for all analog signals is set to 25 kHz. The waveform of LFM is expressed as follows: f u fl sin 2π f x(t ) l 2T 0, others
t t , 0 t T
(14)
where fu and fl are the upper and lower frequency boundaries of the LFM signal, respectively; and T is the signal duration. In the simulation, fl is 7 kHz, fu is 8 kHz, and T is 1 s. The wideband random signal is generated by filtering of Gaussian white noise. The upper and lower frequency boundaries of the wideband random signal are 5.6 kHz and 7.4 kHz, respectively. Figures 2(a) and (b) show the time domain, frequency domain, and power spectrum comparison diagrams of the original and reconstructed wideband random signals with different yCR. Figure 2(a) shows that the reconstructed signal retains the characteristics of the original signal in the time domain, frequency domain, and power spectrum. Therefore, this signal retains the effective components in the frequency domain when the yCR is 5.78. However, as shown in Fig. 2(b), a poor reconstructed effect is obtained when the yCR is 15.4. Table 1 shows the compression evaluation results of the signals. The cross correlation coefficient between the original and the reconstructed signals of three typical underwater acoustic signals is above 0.99, the E is below 0.1, and the yFQSNR becomes large after 5 to 10 times data compression. Table 1 also shows that the single frequency signal obtains an acceptable reconstructed effect when the yCR is 24, the LFM signal obtains a poor reconstructed effect when the yCR is 13.33, and the wideband random signal obtains a poor reconstructed effect when the yCR is 15.4.
J. Cent. South Univ. (2016) 23: 1981−1989
1985
Fig. 2 Comparison diagrams of original (a1, a3, a5, b1, b3, b5) and reconstructed (a2, a4, a6, b2, b4, b6) wideband random signals: (a1, a2, a3, a4, a5, a6) yCR=5.78; (b1, b2, b3, b4, b5, b6) yCR=15.4
J. Cent. South Univ. (2016) 23: 1981−1989
1986
Table 3 Computing time of OMP and IOMP algorithms
Table 1 Compression results of simulation signals Signal form
yCR
E
yFQSNR
R
Wideband random signal_1
5.78
0.03
28.8
0.9917
Wideband random signal_2
15.4
0.4
−12.65 0.8346
Single frequency signal_1
10
0.03
39.41
0.9916
Single frequency signal_2
24
0.08
22.66
0.9829
LFM signal_1
6.41
0.1
19.15
0.9954
Active sonar signal
5.05 0.973212 0.780983 19.75%
LFM signal_2
13.33
0.42
−10.92 0.8533
Ship-radiated noise
6.17 2.879213 2.179925 24.29%
Unlike the LFM and wideband random signals, the single frequency signal is known to have few efficient points in the frequency domain. Therefore, the proposed compression method can obtain a high compression ratio for the single frequency signal, but it obtains a low compression ratio for the LFM and wideband random signals. 3.2 Actual underwater acoustic signal verification This study chooses two types of actual underwater acoustic signals, namely, active sonar signal and ship-radiated noise, to verify the compression method. The comparison diagrams of the original and the reconstructed signals are shown in Fig. 3. The active sonar signal and ship-radiated noise are accurately reconstructed when the yCR values are 5.05 and 6.17, respectively. Table 2 shows the compression evaluation results of the two signals with different yCR. Table 2 shows that the cross correlation coefficient remains high, and the largest normalized frequency domain error remains low for actual underwater acoustic signals when the yCR is about 6. However, an unacceptable reconstructed effect is obtained when the yCR exceeds 10. Table 2 Compression results of actual signals Signal form
yCR
E
yFQSNR
R
Active sonar signal_1
5.05
0.006
52.41
0.9987
Active sonar signal_2
22.74
0.1
2.73
0.9337
Ship-radiated noise_1
6.17
0.015
20.17
0.9961
Ship-radiated noise_2
14.66
0.4
−16.52
0.8000
3.3 IOMP algorithm advantages The computing time of the OMP and IOMP algorithms are shown in Table 3. Only the loop iterative time is shown in Table 3. The IOMP algorithm reduces the computing time by about 20% after achieving a 5 to 10 times larger compression ratio.
4 Array signal compression The above-mentioned signals are all single hydrophone signals, so an estimation experiment on a
Signal form
yCR
Single frequency signal
10
LFM signal
Computing time/s OMP
IOMP
Saved time/%
0.936528 0.711258 24.05%
6.41 0.867352 0.532194 38.64%
Wideband random signal 5.78 3.223983 2.527797 21.59%
12-element circle array DOA is performed in an anechoic tank to verify the validity of the proposed method for array signals. A transducer transmits a continuous wave (CW) pulse signal and a wideband random signal, and a uniform, 12-element circular array is used to receive the signal. The frequency of CW is 2.0 kHz, the sample frequency is 104 kHz, and the pulse width is 20 ms. The bandwidth of the wideband signal ranges from 500 Hz to 2.0 kHz. The transducer and the array have a depth of 3 m and a range of 10 m, whereas the 12-element array has a radius of 0.5 m. The sound speed is assumed to be at 1500 m/s. 4.1 Narrowband signal experimental results The circular array has a scanning range from 0° to 360°. Figure 4 shows the maximum receiving signal of each element. The amplitudes of elements No. 0 to No. 6 are bigger than those of elements No. 7 to No. 11. Therefore, elements No. 0 to No. 6 are considered the forward direction elements. Element No. 3 is chosen as the center element, and element No. 0 is scanned clockwise from the initial direction. After the experimental data are obtained, this study uses the conventional beamformer (CBF) method and the minimum variance distortionless response (MVDR) method for DOA estimation. Figure 5 shows the DOA estimation diagrams of the narrowband signal with different yCR. The sidelobe level of the beam pattern increases with the increase of yCR. A DOA estimation error is obtained for both CBF and MVDR methods when the yCR is greater than 10. The DOA estimation error increases with the yCR, as shown in Fig. 6(a). To consider the effect of the signal noise ratio (SNR), white Gaussian noise is added to the experimental data. Figure 6(b) shows the DOA estimation error with a different SNR when the yCR is 10. When the SNR is too low, the reconstructed signal obtains a DOA estimation error. By contrast, no DOA estimation error is observed when the SNR exceeds 35 dB. 4.2 Wideband random signal experimental results The wideband random signal is decomposed into narrowband components, which are then processed with
J. Cent. South Univ. (2016) 23: 1981−1989
1987
Fig. 3 Comparison diagrams of original (a1, a3, a5, b1, b3, b5) and reconstructed (a2, a4, a6, b2, b4, b6) signals: (a1, a2, a3, a4, a5, a6) Active sonar signal yCR=5.05; (b1, b2, b3, b4, b5, b6) Ship-radiated noise yCR=6.17
1988
J. Cent. South Univ. (2016) 23: 1981−1989
Fig. 4 Maximum receiving signal of a 12-element array
Fig. 6 DOA estimation error diagram with different yCR and SNR: (a) Different CR; (b) Different SNR
Fig. 5 DOA estimation diagrams of narrowband signal with different yCR: (a) CBF method; (b) MVDR method
the compression method. The final result is calculated by averaging of all narrowband components. Figures 7(a) and (b) show the results of the wideband random signal DOA estimation obtained by the CBF and MVDR methods. The proposed compression method is effective for the wideband random signal when the yCR is about 5. A DOA estimation error is obtained when the yCR is greater than 8.
Fig. 7 DOA estimation diagrams of wideband signal with different yCR: (a) CBF method; (b) MVDR method
J. Cent. South Univ. (2016) 23: 1981−1989
5 Conclusions This study introduces an underwater acoustic data compression method based on compressed sensing. The underwater acoustic data are transformed into the sparse domain for data storage at the receiving terminal, and the IOMP algorithm is used to reconstruct the original underwater acoustic data at the data processing terminal. After the compression method is used, a single frequency or narrowband signal can obtain a high CR, whereas a wideband random signal obtains a low CR. When the correlation coefficient of the reconstructed and original signals is over 0.95, increasing the sidelobe level will seldom cause a DOA estimation error. The proposed method can increase the compression ratios for single frequency and wideband random signals by more than 10 and 5 times, respectively. The method can also alleviate the pressure of the underwater acoustic data transmission, reduce system costs, and significantly improve system efficiency. The IOMP algorithm can decrease the computing time by about 20% compared with the original OMP algorithm. The simulation analyses that are performed in this study do not consider the effect of environment noise. A DOA estimation error is observed when the SNR is very low.
1989
[5]
[6]
[7] [8]
[9]
[10]
[11]
[12]
[13]
[14]
References [1] [2] [3]
[4]
SAYOOD K. Introduction to data compression [M]. United States: Newnes Press, 2006. JAIN A K. A fast Karhunen-Loeve transform for a random processes [M]. Chicago: IEEE Computer Society, 1974, 24: 1023−1029. TURCZA P, DUPLAGA M. Low-Power image compression for wireless capsule endoscopy [C]// IEEE International Workshop on Imaging Systems and Techniques-IST. Krakow, Poland, 2007: 1−4. MO Y B, QIU Y B, LIU J Z, LING Y X. A data compression algorithm baseed on adaptive huffman code for wireless sensor
[15]
[16]
networks [C]// Intelligent Computation Technology and Automation (ICICTA). 2011: 3−6. SHEN Yan-chun, GUAN Yu-jun, WANG Fang, LUN Zhi-xin. The investigation of image compress coding based on wavelet transformation [C]// International Conference on Future Information Technology and Management Engineering. 2010: 324−326. BERGER C R, ZHOU S L, PREISIG J C. Sparse channel estimation for multicarrier underwater acoustic communication: From subspace methods to compressed sensing [J]. IEEE Transactions on Signal Processing, 2010, 58(3): 1708−1721. DONOHO D L. Compressed sensing [J]. IEEE Transactions on Information Theory, 2006, 52(4): 1289−1306. TROPP J A, GILBERT A C. Signal recovery from random measurements via orthogonal matching pursuit [J]. IEEE Transactions on Information Theory, 2007, 53(12): 4655−4666. BOASHASH B, O'SHEA P. A methodology for detection and classification of some underwater acoustic signals using timefrequency analysis techniques [J]. IEEE Transactions on Acoustics, Speech and Signal Processing, 1990, 38(11): 1829−1841. CANDES E J, TAO T. Near-optimal signal recovery from random projections: universal encoding strategies [J]. IEEE Transactions on Information Theory, 2006, 52(12): 5406−5425. DAVENPORT M A, WAKIN M B. Analysis of orthogonal matching pursuit using the restricted isometry property [J]. IEEE Transactions on Information Theory, 2010, 56(9): 4395−4401. ZHANG Ge-sen, JIAO Shu-hong, XU Xiao-li, WANG Lan. Compressed sensing and reconstruction with bernoulli matrices [C]// IEEE International Conference on Information and Automation. 2010: 455−460. (in Chinese) GAN Lu, LI Ke-zhi, LING Cong. Golay meets hadamard: Golaypaired hadamard matrices for fast compressed sensing [C]// IEEE Information Theory Workshop. 2012: 637−641. TANG Gui-lin, QIU Yun-ming. Improved least square method apply in ship performance analysis [C]// International Conference on Advanced Computer Theory and Engineering (ICACTE). 2010: 594− 596. MARQUARDT D W. An algorithm for last-squares estimation of nonlinear parameters [J]. Journal of the Society for Industrial and Applied Mathematics, 1963, 11(2): 431−441. LI Ling-zhi, ZOU Bei-ji, ZHU Cheng-zhang. Improved nonconvex optimization model for low-rank matrix recovery [J]. Journal of Central South University, 2015, 22(3): 984−991. (Edited by YANG Hua)