Mixed domain image watermarking algorithm based on multi-correlation statistical modeling

Abstract

The invention discloses a mixed domain image watermarking algorithm based on multi-correlation statistical modeling, and the algorithm comprises the steps: carrying out the secondary non-subsampled wavelet transformation of a carrier image, carrying out the subtraction of a high-frequency sub-band of a first scale and a high-frequency sub-band of a second scale in the same direction, and selecting a wavelet domain difference sub-band with the highest energy as a watermark carrier; performing block processing on the high-energy wavelet domain difference sub-band, selecting a high-entropy block to perform polar-harmonic Fourier decomposition, modifying a UDWT-PHFMs amplitude coefficient of the high-entropy block, and further implementing polar-harmonic Fourier reconstruction and inverse non-subsampled wavelet transform to obtain a watermark-containing image; deriving a probability density function of a vector Sankan-Nel binary Pareto, performing statistical modeling by using a partitioned polar-harmonic Fourier amplitude coefficient of a wavelet domain high-frequency sub-band to obtain corresponding shape parameters, scale parameters and correlation parameters, and performing parameter estimation by using a maximum likelihood method; and according to a maximum likelihood decision criterion, constructing a digital watermark decoder based on an SNBP statistical model, and extracting watermark information bits.

Description

Technical Field

[0001] The present invention belongs to the technical field of copyright protection and content authentication of digital image works, and relates to a digital image watermarking method based on a statistical model, and in particular to a mixed domain image watermarking algorithm based on multi-correlation statistical modeling. Background Art

[0002] With the rapid development of digital technology and network technology, digital multimedia products and their applications have increased significantly, which not only brings many development opportunities to various industries, but also makes human life more convenient. However, the ease and convenience of digital information copying and dissemination has led to various information security issues, especially the copyright issues and content authentication issues of digital image works. As an important supplement to traditional encryption technology, digital watermarking has significantly enhanced the security and reliability of digital information transmission. In particular, digital image watermarking technology can effectively protect the copyright of digital image works and ensure the authenticity of the content. This technology refers to the sender embedding recognizable watermark information into the digital image work before transmitting it, and it will not reduce the visual perception quality of the image work. After receiving the image work, the receiver can accurately extract the watermark information through a specific algorithm, and then determine the source and authenticity of the digital image work based on the watermark information.

[0003] For digital image watermarking technology, imperceptibility and robustness are key indicators to measure its performance. How to achieve the best balance between imperceptibility and robustness is a major challenge. In recent years, scholars at home and abroad have tried to propose the idea of ​​transform domain multiplicative watermarking based on statistical models, which provides a possible solution to effectively solve the problem of a good balance between imperceptibility and robustness. However, the existing statistical model image watermarking schemes are still immature, and there are a series of common problems such as not considering the vulnerability of the single transform coefficient modification strategy, the weak modeling ability of the adopted statistical model, and the abnormal lack of high-performance image watermark detection algorithms. Summary of the invention

[0004] In order to solve the above technical problems existing in the existing related technologies, the present invention proposes a mixed domain image watermarking algorithm based on multi-correlation statistical modeling.

[0005] The technical solution of the present invention is: a mixed domain image watermarking algorithm based on multi-correlation statistical modeling, which is carried out in the following steps:

[0006] Convention: I represents the carrier image; F k (k=1,2,3) represents the k-th high-frequency subband in the wavelet domain scale direction; S k (k=1,2,3) represents the k high frequency subband in the wavelet domain scale direction; D k(k=1,2,3) represents the wavelet domain difference subband in direction k; x,y represent the horizontal and vertical coordinates of the subband coefficients; L represents the number of wavelet domain difference coefficient blocks with larger entropy values; n represents the order of the extreme harmonic Fourier moment; m represents the repetition degree of the extreme harmonic Fourier moment; B l (l=1,2,...,L) represents the UDWT-PHFMs amplitude coefficient set used for watermark embedding; x i represents the original UDWT-PHFMs amplitude coefficient; y i represents the corresponding watermarked UDWT-PHFMs amplitude coefficient; fx i and fy i Respectively represent the original scale two wavelet coefficients and their watermark wavelet coefficients; sx i and sy i They represent the original scale wavelet coefficient and its watermark wavelet coefficient respectively; λ represents the positive weighting factor; w l represents the watermark position to be embedded in the lth wavelet domain difference coefficient block; I′ represents the watermarked image; x1 represents the UDWT-PHFMs amplitude coefficient vector of the wavelet domain coefficient block of scale one; x2 represents the UDWT-PHFMs amplitude coefficient vector of the wavelet domain coefficient block of scale two; Θ represents the model parameter set; θ represents the shape parameter; α1 and α2 represent the scale parameters of scale one and scale two respectively; α0 represents the correlation parameter between scale one and scale two; y represents the UDWT-PHFMs amplitude coefficient vector of the wavelet domain interpolation high entropy block containing the watermark; N represents the number of UDWT-PHFMs amplitude coefficients of the wavelet domain interpolation high entropy block containing the watermark; y 1i and 2i represent the UDWT-PHFMs amplitude coefficients of scale 1 and scale 2 wavelet coefficient blocks respectively; w l ' represents the watermark bit extracted from the lth watermarked wavelet domain difference coefficient block;

[0007] a. Initial Setup

[0008] Input carrier image I and initialize settings;

[0009] b. Watermark Embedding

[0010] b.1 Perform a two-level non-subsampled discrete wavelet transform (UDWT) on the carrier image I to obtain three scale two high frequency subbands F k (k = 1, 2, 3) and 3 scale-high frequency subbands S k (k=1,2,3);

[0011] b.2 Based on the correlation between the coefficients of the scales, the difference between the high-frequency subband coefficients of scale 2 and scale 1 in the same direction is calculated to generate the difference subband D in the wavelet domain. k(k=1,2,3), the calculation method of the difference subband coefficient is as follows:

[0012] D k [x,y]=F k [x,y]-S k [x,y];

[0013] b.3 Select the difference subband with the largest energy as the target subband to embed the digital watermark;

[0014] b.4 Divide the target subband into non-overlapping and equal-sized wavelet domain difference coefficient blocks, calculate the entropy value of the wavelet domain difference coefficient blocks, and then arrange the wavelet domain difference coefficient blocks in descending order according to the entropy value;

[0015] b.5 Calculate the nth-order extreme harmonic Fourier moments (PHFMs) of the first L wavelet domain difference coefficient blocks with larger entropy values, and select the UDWT-PHFMs amplitude coefficient set B in the region of |m|=1,2,3,5&n≤5 l (l=1,2,...,L) is used to hide digital watermark information;

[0016] b.6 Use the multiplicative strategy and modify B l The UDWT-PHFMs amplitude coefficient in is used to embed 1-bit watermark information. The watermark embedding expression is as follows:

[0017]

[0018] b.7 Reconstruct the unmodified and modified UDWT-PHFMs amplitude coefficients by the n-th order extreme harmonic Fourier moment, and map the high entropy wavelet domain difference coefficient block containing the watermark information back to the original position to obtain the wavelet domain high frequency subband containing the watermark;

[0019] b.8 performing a two-stage inverse non-subsampled discrete wavelet transform on the high frequency subband containing the watermark together with the unmodified wavelet domain subband to generate an image I′ containing the watermark;

[0020] c. Sankaran-Nair Bivariate Pareto (SNBP) statistical modeling

[0021] c.1 Perform a two-level non-subsampled discrete wavelet transform on the watermarked image I′ to obtain three scale-two high-frequency subbands and three scale-one high-frequency subbands;

[0022] c.2 Divide all the above scale 2 and scale 1 high frequency subbands into non-overlapping and equal-sized wavelet domain coefficient blocks;

[0023] c.3 For each wavelet domain coefficient block, calculate its n-th order extreme harmonic Fourier moment, and select the UDWT-PHFMs amplitude coefficients in the region of |m|=1,2,3,5&n≤5 to form the training samples;

[0024] c.4 According to the one-to-one relationship between the scale 2 and scale 1 wavelet domain coefficient blocks, the corresponding two groups of UDWT-PHFMs amplitude coefficients (i.e., training samples) are input as two variables into the SNBP probability density function for statistical modeling. The probability density function of SNBP is expressed as follows:

[0025]

[0026] c.5 The maximum likelihood estimation (MLE) method is used to estimate the parameters of the SNBP model. The parameter set can be expressed as:

[0027] Θ={θ,α0,α1,α2};

[0028] d. Decoder construction and watermark extraction

[0029] d.1 Perform a two-level non-subsampled discrete wavelet transform on the watermarked image I′ to obtain three scale-two high-frequency subbands and three scale-one high-frequency subbands;

[0030] d.2 According to the correlation between scale coefficients, the wavelet domain difference subband is calculated, and the wavelet domain difference subband with the largest energy is selected, and the subband is divided into blocks and the high entropy block is determined;

[0031] d.3 Select the scale 2 and scale 1 wavelet coefficient blocks corresponding to the wavelet domain interpolation high entropy block, calculate its n-th order extreme harmonic Fourier moment, and select the UDWT-PHFMs amplitude coefficients in the region of |m|=1,2,3,5&n≤5 to form a test set;

[0032] d.4 According to the maximum likelihood decision (ML) criterion, construct the watermark decoder expression:

[0033]

[0034] d.5 The binary watermark information bits contained in the wavelet domain interpolation high entropy block are extracted as follows:

[0035]

[0036] The invention firstly implements a two-level non-subsampled wavelet transform on the carrier image, and makes a difference between the high-frequency subbands of scale one and scale two in the same direction, thereby obtaining a wavelet domain difference subband, and at the same time selects the wavelet domain difference subband with the highest energy as the watermark carrier; secondly, the high-energy wavelet domain difference subband is processed in blocks, and a high entropy block is selected for extreme harmonic Fourier moment decomposition, and at the same time, a multiplicative embedding strategy is used to modify the high entropy block UDWT-PHFMs amplitude coefficient, and then extreme harmonic Fourier moment reconstruction and inverse non-subsampled wavelet transform are implemented, thereby obtaining a watermarked image; then, according to the correlation between scales in the wavelet domain, the probability density function of the vector Sankaran-Nair binary Pareto (SNBP) is derived, and the block extreme harmonic Fourier moment amplitude coefficients of the high-frequency subband in the wavelet domain are used for statistical modeling, thereby obtaining the corresponding shape parameters, scale parameters and correlation parameters, and at the same time, the maximum likelihood method is used for parameter estimation; finally, according to the maximum likelihood decision criterion, a digital watermark decoder based on the SNBP statistical model is constructed, and the watermark information bits are extracted. Experimental results show that the algorithm of the present invention combines multiple strong correlations between the amplitude coefficients of robust UDWT-PHFMs, constructs a hybrid domain maximum likelihood watermark decoder by establishing a vector Sankaran-Nair Binary Pareto (SNBP) statistical model, and achieves a good balance between imperceptibility and robustness.

[0037] Compared with the prior art, the present invention has the following advantages:

[0038] First, a robust carrier based on the amplitude coefficient of the sub-band harmonic Fourier moment in the wavelet domain difference is constructed and used for watermark embedding and statistical modeling, which improves the robustness of the algorithm.

[0039] Second, combining the multiple strong correlations of UDWT-PHFMs amplitude coefficients between scales in the wavelet domain, a vector Sankaran-Nair Bivariate Pareto (SNBP) statistical model with stronger descriptive ability was established;

[0040] Thirdly, a hybrid domain maximum likelihood watermark decoder based on the vector SNBP statistical model is derived, which achieves a good balance between imperceptibility and robustness. BRIEF DESCRIPTION OF THE DRAWINGS

[0041] Figure 1 This is a result diagram for verifying the robustness of the UDWT-PHFMs amplitude coefficient according to an embodiment of the present invention.

[0042] Figure 2 This is a graph of the Sankaran-Nair Bivariate Pareto (SNBP) PDF fitting results according to an embodiment of the present invention.

[0043] Figure 3 This is a watermarked result image of a grayscale image hiding a 1024-bit watermark according to an embodiment of the present invention.

[0044] Figure 4This is a grayscale image in an embodiment of the present invention containing a 1024-bit watermark image and a 20-fold difference result image of the original image.

[0045] Figure 5 This is a graph showing the 1024-bit watermark extraction results under various attacks according to an embodiment of the present invention.

[0046] Figure 6 The following is a flow chart of watermark embedding according to an embodiment of the present invention.

[0047] Figure 7 This is a flow chart of watermark extraction according to an embodiment of the present invention. DETAILED DESCRIPTION

[0048] The method of the present invention comprises three stages: multiplicative digital watermark embedding, vector Sankaran-Nair binary Pareto modeling, maximum likelihood watermark decoder construction and watermark extraction. Figure 6 , Figure 7 As shown, follow the steps below:

[0049] Convention: I represents the carrier image; F k (k=1,2,3) represents the k-th high-frequency subband in the wavelet domain scale direction; S k (k=1,2,3) represents the k high frequency subband in the wavelet domain scale direction; D k (k=1,2,3) represents the wavelet domain difference subband in direction k; x,y represent the horizontal and vertical coordinates of the subband coefficients; L represents the number of wavelet domain difference coefficient blocks with larger entropy values; n represents the order of the extreme harmonic Fourier moment; m represents the repetition degree of the extreme harmonic Fourier moment; B l (l=1,2,...,L) represents the UDWT-PHFMs amplitude coefficient set used for watermark embedding; x i represents the original UDWT-PHFMs amplitude coefficient; y i represents the corresponding watermarked UDWT-PHFMs amplitude coefficient; fx i and fy i Respectively represent the original scale two wavelet coefficients and their watermark wavelet coefficients; sx i and sy i They represent the original scale wavelet coefficient and its watermark wavelet coefficient respectively; λ represents the positive weighting factor; w lrepresents the watermark position to be embedded in the lth wavelet domain difference coefficient block; I′ represents the watermarked image; x1 represents the UDWT-PHFMs amplitude coefficient vector of the wavelet domain coefficient block of scale one; x2 represents the UDWT-PHFMs amplitude coefficient vector of the wavelet domain coefficient block of scale two; Θ represents the model parameter set; θ represents the shape parameter; α1 and α2 represent the scale parameters of scale one and scale two respectively; α0 represents the correlation parameter between scale one and scale two; y represents the UDWT-PHFMs amplitude coefficient vector of the wavelet domain interpolation high entropy block containing the watermark; N represents the number of UDWT-PHFMs amplitude coefficients of the wavelet domain interpolation high entropy block containing the watermark; y 1i and 2i represent the UDWT-PHFMs amplitude coefficients of scale 1 and scale 2 wavelet coefficient blocks respectively; w l ' represents the watermark bit extracted from the lth watermarked wavelet domain difference coefficient block;

[0050] a. Initial Setup

[0051] Input carrier image I and initialize settings;

[0052] b. Watermark Embedding

[0053] b.1 Perform a two-level non-subsampled discrete wavelet transform (UDWT) on the carrier image I to obtain three scale two high frequency subbands F k (k = 1, 2, 3) and 3 scale-high frequency subbands S k (k=1,2,3);

[0054] b.2 Based on the correlation between the coefficients of the scales, the difference between the high-frequency subband coefficients of scale 2 and scale 1 in the same direction is calculated to generate the difference subband D in the wavelet domain. k (k=1,2,3), the calculation method of the difference subband coefficient is as follows:

[0055] D k [x,y]=F k [x,y]-S k [x,y];

[0056] b.3 Select the difference subband with the largest energy as the target subband to embed the digital watermark;

[0057] b.4 Divide the target subband into non-overlapping and equal-sized wavelet domain difference coefficient blocks, calculate the entropy value of the wavelet domain difference coefficient blocks, and then arrange the wavelet domain difference coefficient blocks in descending order according to the entropy value;

[0058] b.5 Calculate the nth-order extreme harmonic Fourier moments (PHFMs) of the first L wavelet domain difference coefficient blocks with larger entropy values, and select the UDWT-PHFMs amplitude coefficient set B in the region of |m|=1,2,3,5&n≤5 l (l=1,2,...,L) is used to hide digital watermark information;

[0059] b.6 Use the multiplicative strategy and modify B l The UDWT-PHFMs amplitude coefficient in is used to embed 1-bit watermark information. The watermark embedding expression is as follows:

[0060]

[0061] b.7 Reconstruct the unmodified and modified UDWT-PHFMs amplitude coefficients by the n-th order extreme harmonic Fourier moment, and map the high entropy wavelet domain difference coefficient block containing the watermark information back to the original position to obtain the wavelet domain high frequency subband containing the watermark;

[0062] b.8 performing a two-stage inverse non-subsampled discrete wavelet transform on the high frequency subband containing the watermark together with the unmodified wavelet domain subband to generate an image I′ containing the watermark;

[0063] c. Sankaran-Nair Bivariate Pareto (SNBP) statistical modeling

[0064] c.1 Perform a two-level non-subsampled discrete wavelet transform on the watermarked image I′ to obtain three scale-two high-frequency subbands and three scale-one high-frequency subbands;

[0065] c.2 Divide all the above scale 2 and scale 1 high frequency subbands into non-overlapping and equal-sized wavelet domain coefficient blocks;

[0066] c.3 For each wavelet domain coefficient block, calculate its n-th order extreme harmonic Fourier moment, and select the UDWT-PHFMs amplitude coefficients in the region of |m|=1,2,3,5&n≤5 to form the training samples;

[0067] c.4 According to the one-to-one relationship between the scale 2 and scale 1 wavelet domain coefficient blocks, the corresponding two groups of UDWT-PHFMs amplitude coefficients (i.e., training samples) are input as two variables into the SNBP probability density function for statistical modeling. The probability density function of SNBP is expressed as follows:

[0068]

[0069] c.5 The maximum likelihood estimation (MLE) method is used to estimate the parameters of the SNBP model. The parameter set can be expressed as:

[0070] Θ={θ,α0,α1,α2};

[0071] d. Decoder construction and watermark extraction

[0072] d.1 Perform a two-level non-subsampled discrete wavelet transform on the watermarked image I′ to obtain three scale-two high-frequency subbands and three scale-one high-frequency subbands;

[0073] d.2 According to the correlation between scale coefficients, the wavelet domain difference subband is calculated, and the wavelet domain difference subband with the largest energy is selected, and the subband is divided into blocks and the high entropy block is determined;

[0074] d.3 Select the scale 2 and scale 1 wavelet coefficient blocks corresponding to the wavelet domain interpolation high entropy block, calculate its n-th order extreme harmonic Fourier moment, and select the UDWT-PHFMs amplitude coefficients in the region of |m|=1,2,3,5&n≤5 to form a test set;

[0075] d.4 According to the maximum likelihood decision (ML) criterion, construct the watermark decoder expression:

[0076]

[0077] d.5 The binary watermark information bits contained in the wavelet domain interpolation high entropy block are extracted as follows:

[0078]

[0079] Experimental test and parameter setting:

[0080] The experimental environment is MATLAB R2011a, and the grayscale images are all 512×512. Download address:

[0081] http: / / decsai.ugr.es / cvg / dbimagenes / index.php.

[0082] Figure 1 This is a result diagram for verifying the robustness of the UDWT-PHFMs amplitude coefficient according to an embodiment of the present invention.

[0083] Figure 1 (a) No attack; (b) JPEG 30; (c) JPEG 70; (d) AWGN 10; (e) AWGN 30; (f) Median 3×3; (g) Median 7×7; (h) Salt 0.03; (i) Salt 0.1; (j)Gamma1.5; (k)Gamma0.75; (l)Gaussian3×3.

[0084] Figure 2 This is a graph of the Sankaran-Nair Bivariate Pareto (SNBP) PDF fitting results according to an embodiment of the present invention.

[0085] Figure 2 (a) Joint probability density distribution diagram; (b) SNBP distribution diagram; Contour plot of SNBP distribution (solid line) fitting the bivariate empirical PDF (dashed line).

[0086] Figure 3 This is a watermarked result image of a grayscale image hiding a 1024-bit watermark according to an embodiment of the present invention.

[0087] Figure 3 (a) original image Lena; (b) original image Baboon; (c) original image Peppers; (d) watermarked image Lena; (e) watermarked image Baboon; (f) watermarked image Peppers.

[0088] Figure 4 This is a grayscale image in an embodiment of the present invention containing a 1024-bit watermark image and a 20-fold difference result of the original image.

[0089] Figure 4 (a) Lena - 20 times difference image; (b) Baboon - 20 times difference image; (c) Peppers - 20 times difference image.

[0090] Figure 5 This is a graph showing the 1024-bit watermark extraction results under various attacks according to an embodiment of the present invention.

[0091] Figure 5 (a) JPEG compression; (b) Gaussian noise; (c) median filter; (d) salt and pepper noise; (e) Gaussian filter; (f) rotation.

[0092] Figure 5 Comparative literature used: M Amini, MO Ahmad, MNS Swamy. Arobust multibitmultiplicative watermark decoder using vector-based hidden Markov model inwavelet domain. IEEE Transactions on Circuits&Systems for Video Technology, 2018, 28(2): 402-413.

[0093] Xiang-yang Wang,Jing Tian,Jia-lin Tian,Pan-pan Niu,Hong-yingYang.Statistical image watermarking using local RHFMs magnitudes and Betaexponential distribution.Journal ofVisual Communication and ImageRepresentation,2021,77:103123。

Claims

1. A mixed domain image watermarking algorithm based on multi-correlation statistical modeling, characterized by Follow these steps: Convention: I represents the carrier image; F k (k=1,2,3) represents the k-th high-frequency subband in the wavelet domain scale direction; S k (k=1,2,3) represents the k high frequency subband in the wavelet domain scale direction; D k (k=1,2,3) represents the wavelet domain difference subband in direction k; x,y represent the horizontal and vertical coordinates of the subband coefficients; L represents the number of wavelet domain difference coefficient blocks with larger entropy values; n represents the order of the extreme harmonic Fourier moment; m represents the repetition degree of the extreme harmonic Fourier moment; B l (l=1,2,...,L) represents the UDWT-PHFMs amplitude coefficient set used for watermark embedding; x i represents the original UDWT-PHFMs amplitude coefficient; y i represents the corresponding watermarked UDWT-PHFMs amplitude coefficient; fx i and fy i Respectively represent the original scale two wavelet coefficients and their watermark wavelet coefficients; sx i and sy i They represent the original scale wavelet coefficient and its watermark wavelet coefficient respectively; λ represents the positive weighting factor; w l represents the watermark position to be embedded in the lth wavelet domain difference coefficient block; I′ represents the watermarked image; x1 represents the UDWT-PHFMs amplitude coefficient vector of the wavelet domain coefficient block of scale one; x2 represents the UDWT-PHFMs amplitude coefficient vector of the wavelet domain coefficient block of scale two; Θ represents the model parameter set; θ represents the shape parameter; α1 and α2 represent the scale parameters of scale one and scale two respectively; α0 represents the correlation parameter between scale one and scale two; y represents the UDWT-PHFMs amplitude coefficient vector of the wavelet domain interpolation high entropy block containing the watermark; N represents the number of UDWT-PHFMs amplitude coefficients of the wavelet domain interpolation high entropy block containing the watermark; y 1i and 2i represent the UDWT-PHFMs amplitude coefficients of scale 1 and scale 2 wavelet coefficient blocks respectively; w l ' represents the watermark bit extracted from the lth watermarked wavelet domain difference coefficient block; a. Initial Setup Input carrier image I and initialize settings; b. Watermark Embedding b.1 Perform a two-level non-subsampled discrete wavelet transform on the carrier image I to obtain three scale two-high frequency subbands F k (k = 1, 2, 3) and 3 scale-high frequency subbands S k (k=1,2,3); b.2 Based on the correlation between the coefficients of the scales, the difference between the high-frequency subband coefficients of scale 2 and scale 1 in the same direction is calculated to generate the difference subband D in the wavelet domain. k (k=1,2,3), the calculation method of the difference subband coefficient is as follows: D k [x,y]=F k [x,y]-S k [x,y]; b.3 Select the difference subband with the largest energy as the target subband to embed the digital watermark; b.4 Divide the target subband into non-overlapping and equal-sized wavelet domain difference coefficient blocks, calculate the entropy value of the wavelet domain difference coefficient blocks, and then arrange the wavelet domain difference coefficient blocks in descending order according to the entropy value; b.5 Calculate the nth-order extreme harmonic Fourier moment of the first L wavelet domain difference coefficient blocks with larger entropy values, and select the UDWT-PHFMs amplitude coefficient set B in the region of |m|=1,2,3,5&n≤5 l (l=1,2,...,L) is used to hide digital watermark information; b.6 Use the multiplicative strategy and modify B l The UDWT-PHFMs amplitude coefficient in is used to embed 1-bit watermark information. The watermark embedding expression is as follows: b.7 Reconstruct the unmodified and modified UDWT-PHFMs amplitude coefficients by the n-th order extreme harmonic Fourier moment, and map the high entropy wavelet domain difference coefficient block containing the watermark information back to the original position to obtain the wavelet domain high frequency subband containing the watermark; b.8 performing a two-stage inverse non-subsampled discrete wavelet transform on the high-frequency subband containing the watermark together with the unmodified wavelet domain subband to generate an image I′ containing the watermark; c. Sankaran-Nair Bivariate Pareto Statistical Modeling c.1 Perform a two-level non-subsampled discrete wavelet transform on the watermarked image I′ to obtain three scale-two high-frequency subbands and three scale-one high-frequency subbands; c.2 Divide all the above scale 2 and scale 1 high frequency subbands into non-overlapping and equal-sized wavelet domain coefficient blocks; c.3 For each wavelet domain coefficient block, calculate its n-th order extreme harmonic Fourier moment, and select the UDWT-PHFMs amplitude coefficients in the region of m|=1,2,3,5&n≤5 to form the training samples; c.4 According to the one-to-one relationship between the scale 2 and scale 1 wavelet domain coefficient blocks, the corresponding two groups of UDWT-PHFMs amplitude coefficients are input as two variables into the SNBP probability density function for statistical modeling. The probability density function of SNBP is expressed as follows: c.5 The maximum likelihood estimation method is used to estimate the parameters of the SNBP model. The parameter set can be expressed as: Θ={θ,α0,α1,α2}; d. Decoder construction and watermark extraction d.1 Perform a two-level non-subsampled discrete wavelet transform on the watermarked image I′ to obtain three scale-two high-frequency subbands and three scale-one high-frequency subbands; d.2 According to the correlation between scale coefficients, the wavelet domain difference subband is calculated, and the wavelet domain difference subband with the largest energy is selected, and the subband is divided into blocks and the high entropy block is determined; d.3 Select the scale 2 and scale 1 wavelet coefficient blocks corresponding to the wavelet domain interpolation high entropy block, calculate its n-th order extreme harmonic Fourier moment, and select the UDWT-PHFMs amplitude coefficients in the region of m|=1,2,3,5&n≤5 to form a test set; d.4 According to the maximum likelihood decision criterion, construct the watermark decoder expression: d.5 The binary watermark information bits contained in the wavelet domain interpolation high entropy block are extracted as follows:

Images