CN116883226A - NMF decomposition-based DEM zero watermark method, device and medium - Google Patents

Abstract

The application belongs to the technical field of copyright protection, and particularly discloses a DEM zero watermarking method, device and medium based on NMF decomposition, which are used for solving the problem of poor robustness of DEM data when facing large-area cutting attack, and particularly storing a coefficient matrix generated by the decomposition of the DEM data and final copyright information into a third party mechanism by means of the characteristic of restorability of NMF decomposition technology; and constructing a gradient feature matrix by using HU invariant moment, carrying out exclusive OR with copyright information to obtain zero watermark information, and extracting a base matrix from data to be detected when detecting the copyright information, so as to reversely calculate original data. Experimental results show that the algorithm has stronger robustness in the face of common watermark attack and large-area clipping attack, and the constructed watermark has better uniqueness.

Description

NMF decomposition-based DEM zero watermark method, device and medium

Technical Field

The application belongs to the technical field of copyright protection, and particularly relates to a DEM zero watermark method, device and medium based on NMF decomposition.

Background

The digital elevation model (Digital Elevation Model, DEM) is an important component of national geographic information resources, and has wide application in the fields of geology, homeland, land resource planning and the like. In recent years, as DEM data use requirements increase, data security problems during transmission and storage are also more and more serious. At present, research on embedded watermark technology has achieved a lot of results, but how to design a copyright protection method suitable for DEM data is always an important direction of research in related fields.

To improve the robustness of zero watermark against geometric attacks, wang et al [ Wang C P, wang X Y, chen X J, zhang c.robustzero-watermarking algorithm based on polar complex exponential transform and Logistic mapping [ J ]. Multimedia Tools and Applications,2017,76 (24) ] propose a new zero watermark algorithm based on PCET and Logistic mapping, which improves the robustness of data by using the geometrical invariance of PCET and improves the security of the algorithm by using the sensitivity of Logistic mapping to initial values. Qu Changbo [ Qu Changbo, wu Deyang ] computer applied research on strong robust zero watermark algorithm based on Curvelet-DSVD and visual cryptography [ J ]. 2019,36 (02): 532-537.DOI:10.19734/J. Issn.1001-3695.2017.10.0936 ] et al propose zero watermark algorithm based on Curvelet-DSVD and visual cryptography by partitioning the transformed low frequency band information and performing dual singular value decomposition, constructing a zero watermark feature matrix by comparing the partition with the whole mean, obtaining a share by using the visual cryptography, scrambling the share, and xoring to generate a zero watermark. The algorithm can better represent the curve characteristics of the natural image and further improve the robustness. Liu Mojun [ Liu Mojun, sun Saiyu, qu Haicheng ], et al, a zero watermark algorithm against geometric rotation attack [ J ]. Computer applied research, 2019,36 (09): 2803-2808 ], et al, propose a zero watermark algorithm combining the SIFT algorithm with wavelet transform, correcting the pixel distortion by rotation of the SIFT transform, after determining a safe region with near lossless center, wavelet transforming the region to extract a low frequency region, and comparing the block with the whole mean value to construct a zero watermark. Through experiments, compared with the previous algorithm such as GH torque and the like, the algorithm has a certain improvement. Gong et al [ Gong C, liu J, gong M, et al, bhatti Uzair Aslam, ma J. Robust medium zero-watermarking algorithm based on Residual DenseNet J. IET biomerics, 2022,11 (6) ] based on deep learning, propose a method for enhancing extraction of different features by means of a residual dense network, applying a hopping connection and a new objective function to the low frequency band features after frequency domain transformation, combining the binarized hash vector with encrypted watermark information, and generating a zero watermark. Experiments prove that the algorithm has better robustness against geometric attacks.

DEM data is a topographical representation of an area, and in the process of generating data, data for a large area is often generated simultaneously. But only a portion of the data is often needed during use. When the data is revealed, the thief only needs to cut out the needed part, and discard the remaining large area, which may cause that the detection fails due to the small remaining part when the copyright is traced back to the part, so that the copyright attribution cannot be determined. In addition, if the copyright owner can form zero watermark in the whole area, and after clipping and distributing, the clipped sub-blocks can still judge the attribution of the copyright, so that great cost is reduced.

Disclosure of Invention

The present application has been made to solve the above-mentioned problems occurring in the prior art. Therefore, a method, a device and a medium for zero watermark of DEM based on NMF decomposition are needed to solve the problem of poor robustness of DEM data in large-area clipping attack.

According to a first technical scheme of the application, there is provided a DEM zero watermarking method based on NMF decomposition, the method comprising:

acquiring a copyright image, and scrambling the copyright image to obtain information M of the encrypted copyright image;

extracting an elevation matrix from original DEM data, performing non-overlapping blocking on the elevation matrix, converting each sub-block into a column vector to form a new matrix B, and performing NMF decomposition on the matrix B to obtain a base matrix W and a coefficient matrix H;

calculating the gradient of each grid in the original DEM data to obtain a gradient image C;

non-overlapping blocking is carried out on the gradient image C, at least three HU invariant moments of each small block and the gradient image C are calculated respectively, and the three high-order moment moments are averaged;

comparing the HU invariant moment average value of the gradient image C with the HU invariant moment average value of each small block, and recording as 1 when the whole is larger and as 0 when the whole is smaller to obtain a binary characteristic image D;

and (3) carrying out exclusive OR on the binary characteristic image D and the encrypted copyright image information M to obtain zero watermark information WM, and converting the original point coordinates, the grid row number and the operation time in the original DEM data into a secret key K.

Further, the copyrighted image is scrambled by the following formula:

where mod is used to extract the modulus in the formula, N is the row width of square copyright data, (x ', y') is the processed coordinates, and (x, y) is the coordinates in the original information.

Further, the performing NMF decomposition on the matrix B to obtain a base matrix W and a coefficient matrix H includes:

based on the matrix B, randomly initializing to respectively obtain a base matrix W and a coefficient matrix H, wherein W is more than 0 and H is more than or equal to 0.

The base matrix W is iteratively updated according to an update rule shown in the following equation (5):

wherein i represents the row number of the original matrix, j represents the column number of the original matrix, and B is approximately equal to W and H after multiple iterations;

normalizing the rows and columns of the base matrix W;

iteratively updating the coefficient matrix H through the base matrix W according to an updating rule shown in a formula (5);

and comparing the calculated error with a first threshold value, if the calculated error is smaller than or equal to the first threshold value, finishing operation, outputting an operation result, and if the calculated error is larger than the threshold value, iteratively updating the base matrix W according to an updating rule shown in the formula (5) until the calculated error is smaller than or equal to the first threshold value.

Further, the update rule is determined by:

taking the Euclidean distance square as an objective function, extracting a matrix B of original DEM data, calculating a reconstructed matrix WH, and judging whether iteration is continued or not by calculating a distance value between the matrix B and the matrix WH and comparing the distance value with a second threshold value, wherein a Euclidean distance square formula is shown in a formula (2):

D(W,H)＝||B-WH|| ² ＝∑ _ij (B _ij -(WH) _ij ) ² (2)

and respectively deriving W and H on two sides of the formula (2), wherein the derived result is shown as the formula (3):

wherein T is a rank-converting operator, and k is a decomposition dimension;

by means of the gradient descent method, an update rule is deduced as shown in formula (4):

order theThe update method shown in the formula (4) is simplified from addition to multiplication, and as shown in the formula (5), the update rules of the base matrix W and the coefficient matrix H are obtained respectively.

Further, the method also includes zero watermark extraction.

Further, the zero watermark extraction includes:

reading DEM data A' to be detected, and obtaining zero watermark information WM, a coefficient matrix H and a secret key K;

judging whether the DEM data A' is cut or not according to the original point coordinates and the grid number in the secret key; if the cutting is performed, the step 3 is performed, and if the cutting is not performed, the step 4 is performed;

under the condition that DEM data A ' is cut, 0 is complemented to a cutting part according to key information, then the data A ' is converted into column vectors to be extracted into a matrix I, corresponding column vectors are obtained from H according to the matrix I, H ' is formed, and W=I [ H ' '] ^-1 Calculating a base matrix W, multiplying W by a coefficient matrix H, and obtaining reconstruction data B';

under the condition that DEM data A 'is not cut, performing invariant moment feature extraction on the data A' to be detected or the reconstruction data B 'to obtain a feature image D';

and exclusive-or the characteristic image D ' and the zero watermark information WM, and performing Arnold scrambling on the characteristic image D ' to obtain the extracted information M ' of the encrypted copyright image.

Further, the information M of the encrypted copyright image is compared with the information M' of the encrypted copyright image, and copyright authentication is performed.

Further, the clipping part is complemented with 0 according to the key information, then the data a ' is converted into a column vector to be extracted into a matrix I, the corresponding column vector is obtained from H according to I, and H ' is formed, and the data a ' is represented by the formula w=i [ H ' '] ^-1 Calculating a base matrix W, multiplying W by a coefficient matrix H, and obtaining reconstruction data B', wherein the method comprises the following steps:

the basic model of NMF decomposition is shown as formula (6):

B _m*n ＝W _m*r *H _r*n (6)

wherein m, n are the dimensions of the non-negative matrix B, and r is the dimension of NMF decomposition;

b _i is the element of the ith column in the non-negative matrix B, h _i For the element in the ith column in the coefficient matrix H, there are:

b _i ＝W _m*r *h _i i∈[1,n] (7)

when the value of the sub-block to be cut is set to 0 after the cutting attack, there are:

B _m*n ＝[b ₁ ,…b _s ,0,…,0,b _t ,…,b _n ] (8)

namely:

B′ _k ＝[b ₁ ,…b _s ,b _t ,…,b _n ] (9)

wherein k represents the number of column vectors and the number of the remaining complete sub-blocks of the original data matrix after the clipping attack, when k is more than or equal to r, namely the number of the matrix sub-blocks after the clipping attack is more than or equal to a dimension value r, the complete base matrix can be recovered, and the inverse operation reconstruction of NMF decomposition is realized;

in the case of B ' being known, both sides are simultaneously multiplied by [ H ', based on B ' =w×h ' '] ^-1 The method comprises the following steps:

W＝B′*[H′] ^-1 (10) In the case of the known coefficient matrix H, the original matrix is reconstructed according to the above procedure.

According to a second aspect of the present application, there is provided an NMF decomposition-based DEM zero watermarking method apparatus, the apparatus comprising a processor configured to:

acquiring a copyright image, and scrambling the copyright image to obtain information M of the encrypted copyright image;

extracting an elevation matrix from original DEM data, performing non-overlapping blocking on the elevation matrix, converting each sub-block into a column vector to form a new matrix B, and performing NMF decomposition on the matrix B to obtain a base matrix W and a coefficient matrix H;

calculating the gradient of each grid in the original DEM data to obtain a gradient image C;

non-overlapping blocking is carried out on the gradient image C, at least three HU invariant moments of each small block and the gradient image C are calculated respectively, and the three high-order moment moments are averaged;

comparing the HU invariant moment average value of the gradient image C with the HU invariant moment average value of each small block, and recording as 1 when the whole is larger and as 0 when the whole is smaller to obtain a binary characteristic image D;

and (3) carrying out exclusive OR on the binary characteristic image D and the encrypted copyright image information M to obtain zero watermark information WM, and converting the original point coordinates, the grid row number and the operation time in the original DEM data into a secret key K.

Further, the processor is further configured to: scrambling the copyrighted image by the following formula:

where mod is used to extract the modulus in the formula, N is the row width of square copyright data, (x ', y') is the processed coordinates, and (x, y) is the coordinates in the original information.

Further, the processor is further configured to:

based on the matrix B, randomly initializing to respectively obtain a base matrix W and a coefficient matrix H, wherein W is more than 0 and H is more than or equal to 0.

The base matrix W is iteratively updated according to an update rule shown in the following equation (5):

wherein i represents the row number of the original matrix, j represents the column number of the original matrix, and B is approximately equal to W and H after multiple iterations;

normalizing the rows and columns of the base matrix W;

iteratively updating the coefficient matrix H through the base matrix W according to an updating rule shown in a formula (5);

and comparing the calculated error with a first threshold value, if the calculated error is smaller than or equal to the first threshold value, finishing operation, outputting an operation result, and if the calculated error is larger than the threshold value, iteratively updating the base matrix W according to an updating rule shown in the formula (5) until the calculated error is smaller than or equal to the first threshold value.

Further, the processor is further configured to: taking the Euclidean distance square as an objective function, extracting a matrix B of original DEM data, calculating a reconstructed matrix WH, and judging whether iteration is continued or not by calculating a distance value between the matrix B and the matrix WH and comparing the distance value with a second threshold value, wherein a Euclidean distance square formula is shown in a formula (2):

D(W,H)＝||B-WH|| ² ＝∑ _ij (B _ij -(WH) _ij ) ² (2)

and respectively deriving W and H on two sides of the formula (2), wherein the derived result is shown as the formula (3):

wherein T is a rank-converting operator, and k is a decomposition dimension;

by means of the gradient descent method, an update rule is deduced as shown in formula (4):

order theThe update method shown in the formula (4) is simplified from addition to multiplication, and as shown in the formula (5), the update rules of the base matrix W and the coefficient matrix H are obtained respectively.

Further, the processor is further configured to: and (5) extracting the zero watermark.

Further, the processor is further configured to: reading DEM data A' to be detected, and obtaining zero watermark information WM, a coefficient matrix H and a secret key K;

judging whether the DEM data A' is cut or not according to the original point coordinates and the grid number in the secret key; if the cutting is performed, the step 3 is performed, and if the cutting is not performed, the step 4 is performed;

under the condition that DEM data A ' is cut, 0 is complemented to a cutting part according to key information, then the data A ' is converted into column vectors to be extracted into a matrix I, corresponding column vectors are obtained from H according to the matrix I, H ' is formed, and W=I [ H ' '] ^-1 Calculating a base matrix W, multiplying W by a coefficient matrix H, and obtaining reconstruction data B';

under the condition that DEM data A 'is not cut, performing invariant moment feature extraction on the data A' to be detected or the reconstruction data B 'to obtain a feature image D';

and exclusive-or the characteristic image D ' and the zero watermark information WM, and performing Arnold scrambling on the characteristic image D ' to obtain the extracted information M ' of the encrypted copyright image.

Further, the processor is further configured to compare the information M of the encrypted copyrighted image with the information M' of the encrypted copyrighted image for copyright authentication.

Further, the processor is further configured to: the basic model of NMF decomposition is shown as formula (6):

B _m*n ＝W _m*r *H _r*n (6)

wherein m, n are the dimensions of the non-negative matrix B, and r is the dimension of NMF decomposition;

b _i is the element of the ith column in the non-negative matrix B, h _i For the element in the ith column in the coefficient matrix H, there are:

b _i ＝W _m*r *h _i i∈[1,n] (7)

when the value of the sub-block to be cut is set to 0 after the cutting attack, there are:

B _m*n ＝[b ₁ ,…b _s ,0,…,0,b _t ,…,b _n ] (8)

namely:

B′ _k ＝[b ₁ ,…b _s ,b _t ,…,b _n ] (9)

wherein k represents the number of column vectors and the number of the remaining complete sub-blocks of the original data matrix after the clipping attack, when k is more than or equal to r, namely the number of the matrix sub-blocks after the clipping attack is more than or equal to a dimension value r, the complete base matrix can be recovered, and the inverse operation reconstruction of NMF decomposition is realized;

in the case of B ' being known, both sides are simultaneously multiplied by [ H ', based on B ' =w×h ' '] ^-1 The method comprises the following steps:

W＝B′*[H′] ^-1 (10)

in the case of the known coefficient matrix H, the original matrix is reconstructed according to the above procedure.

According to a third aspect of the present application, there is provided a readable storage medium storing one or more programs executable by one or more processors to implement the method as described above.

The application has at least the following beneficial effects:

the application stores the coefficient matrix generated by decomposition and the last copyright information into a third party mechanism by means of the characteristic of restorability of NMF decomposition technology, and extracts the base matrix from the data to be detected when detecting the copyright information, thereby reversely calculating the original data. The robustness against other geometric attacks is improved by extracting HU invariant moment from the topography factors. Experimental results show that the algorithm has stronger robustness in the face of common watermark attack and large-area clipping attack, and the built watermark has higher uniqueness.

Drawings

Fig. 1 shows a flowchart of a DEM zero watermarking method based on NMF decomposition according to an embodiment of the application;

FIG. 2 shows an Arnold scrambling transformation effect diagram in accordance with an embodiment of the present application;

fig. 3 shows an experimental data diagram according to an embodiment of the present application, in which (a) represents original data and (b) represents a copyrighted image;

fig. 4 shows the clipping effect of different intensity clipping and the recovery effect of NMF transformation and watermark extraction effect after facing such an attack, according to an embodiment of the application;

FIG. 5 shows translation attack experimental results according to an embodiment of the present application;

FIG. 6 shows the results of a rotation attack experiment according to an embodiment of the present application;

FIG. 7 shows scaling attack experimental results according to an embodiment of the present application;

FIG. 8 illustrates the results of an altitude translation attack experiment according to an embodiment of the present application;

FIG. 9 shows Gaussian noise attack experimental results according to an embodiment of the application;

fig. 10 shows the results of a salt and pepper noise attack experiment according to an embodiment of the present application;

FIG. 11 shows a graph of comparative test data in which (a) represents raw data 1 and (b) represents comparative data 2, according to an embodiment of the present application;

FIG. 12 shows a uniqueness verification comparison graph in accordance with an embodiment of the application;

FIG. 13 shows a plot of BER values versus crop attack for different algorithms according to embodiments of the present application;

FIG. 14 shows NC value versus graph for different algorithms facing a clipping attack, according to an embodiment of the application.

Detailed Description

The present application will be described in detail below with reference to the drawings and detailed description to enable those skilled in the art to better understand the technical scheme of the present application. Embodiments of the present application will be described in further detail below with reference to the drawings and specific examples, but not by way of limitation. The order in which the steps are described herein by way of example should not be construed as limiting if there is no necessity for a relationship between each other, and it should be understood by those skilled in the art that the steps may be sequentially modified without disrupting the logic of each other so that the overall process is not realized.

Aiming at the problems of illegal acquisition, tampering and copyright attribution possibly existing in the sharing process of high-precision DEM data and the copyright traceability requirement of the data after the large-space clipping under specific conditions, the embodiment of the application provides a DEM zero watermarking method based on NMF decomposition. The method adopts the idea of data reconstruction, and utilizes the characteristic of restorability of NMF decomposition, and utilizes the characteristic of unchanged relative position relation of DEM data in the face of attack, thereby realizing the embedding and detection of copyright information under the condition of not changing the data precision.

The specific flow chart of the method is shown in fig. 1, and the basic idea is as follows: firstly, NMF is used for decomposing and processing original data to obtain a coefficient matrix and a base matrix, secondly, the gradient of the original data is calculated, the gradient is divided into non-overlapping blocks, HU invariant moment processing is respectively carried out on the whole gradient and the blocks, the sizes of the whole gradient and the blocks are compared to form a binary matrix, the binary matrix is exclusive-or with encrypted copyright information to generate a zero watermark, and the zero watermark and the coefficient matrix are stored in a third party mechanism together.

The NMF decomposition-based DEM zero watermarking method specifically comprises the following step one and step two.

Step one, watermark information generation.

Step 1.1: in order to facilitate exclusive or with DEM data, raster images are selected as copyright data.

Step 1.2: in order to improve the security and confidentiality of the copyright information, arnold scrambling is selected to be used, and the copyright image is scrambled to achieve the encryption effect.

In some embodiments, the format of the copyright information selected in the present document is raster data, and in order to improve the security of the copyright information, the copyright information is scrambled.

Arnold scrambling is one of common image scrambling algorithms, has periodicity, and can change data to be encrypted from ordered to unordered by processing (x, y) coordinates in original information through the algorithm, so that the scrambling effect is achieved. And because the copyright data needs to be restored after watermark extraction, the data can be changed from unordered to ordered by continuing processing. The scrambling formula is shown in (1):

wherein mod is used to extract a modulus value in the formula, N is a row-column width of square copyright data, and (x ', y') is a processed coordinate, in this embodiment, a 64×64 image is used as a copyright image, and the period is 48, and the specific flow is shown in (a), (b), (c), and (d) in fig. 2.

Step 1.3: and acquiring information M of the encrypted copyright image to obtain the row number r and the period T of the encrypted copyright image.

Step two, watermark information embedding.

Step 2.1: and extracting an elevation matrix from the original DEM data, and carrying out 16-by-16 non-overlapping blocking on the matrix, wherein each sub-block is converted into a column vector to form a new matrix B. And (3) performing NMF decomposition on the B to obtain a base matrix W and a coefficient matrix H.

In some embodiments, NMF decomposition is applied to DEM data, the main solution idea being: and introducing an objective function for calculating errors between the matrix extracted from the original DEM data and the matrix reconstructed and restored by NMF, iterating continuously, and stopping after the threshold is reached by given threshold constraint to obtain an NMF decomposition result.

The Euclidean distance square is used as an objective function, a matrix B of original DEM data is extracted, a reconstructed matrix WH is calculated, and whether iteration is continued or not is judged by comparing a distance value with a threshold value through calculating the distance between the two matrices. The formula of the square of the euclidean distance is shown in formula (2), and the value of the square of the euclidean distance is 0 if and only if B- (WH) =0, at this time, B and w×h are infinitely close.

D(W,H)＝||B-WH|| ² ＝∑ _ij (B _ij -(WH) _ij ) ² (2)

W and H are respectively derived on two sides of the formula (2), and the derived result is shown as the formula (3).

Where T is the rank-converting operator and k is the decomposition dimension. By means of the gradient descent method, an update rule is deduced as shown in formula (4).

To make the formula non-negative, letThe update rule shown in equation (4) can be simplified from addition to multiplication, and as shown in (5), the update rule of the base matrix W and the coefficient matrix H can be obtained, respectively.

Where i represents the number of rows of the original matrix and j represents the number of columns of the original matrix. After several iterations b≡w will be H.

The NMF decomposition method is as shown above, after the corresponding updating rule is determined by adopting the Euclidean distance square method, different thresholds are set according to different application scenes, and the thresholds are used as termination conditions of NMF decomposition iteration, or the iteration times can be regulated, and the automatic stop can be carried out after the iteration times are reached.

In summary, the implementation steps of the NMF decomposition method are as follows:

(1) Randomly initializing to obtain a base matrix W and a coefficient matrix H respectively, wherein W is more than 0 and H is more than or equal to 0.

(2) And iteratively updating the base matrix W according to the deduced updating rule.

(3) Normalizing rows and columns of the base matrix W, e.g.As shown.

(4) The coefficient matrix H is iteratively updated by the basis matrix W according to the update law.

(5) Comparing the calculated error with a threshold value, judging whether iteration is ended, if the condition is met, finishing operation, outputting an operation result, and if the condition is not met, turning to the step (2), and repeating calculation.

Step 2.2: and calculating the gradient of each grid of the original data to obtain a gradient image C.

Step 2.3: and (3) carrying out 16×16 non-overlapping block division on the gradient image C, respectively calculating seven HU invariant moments of each small block and the whole gradient image C, and averaging the three high-order moment moments.

With respect to the generation of gradient invariant moment, the attack on DEM data usually does not change the height relationship of adjacent grids, so that the magnitude relationship between adjacent gradients does not change after the gradient is generated. The gradient of each grid in eight directions is calculated and averaged to obtain the gradient value of the grid. In order to further highlight the characteristics, the 5 th, 6 th and 7 th invariant moment in the HU invariant moment formula are selected, three invariant moment are calculated respectively in calculation, and finally, the invariant moment characteristic values of the grid are obtained through averaging.

Step 2.4: comparing the HU moment average value of the gradient image C with the HU moment average value of each small block, and marking as 1 when the whole is larger, and marking as 0 when the whole is smaller, thereby obtaining a binary characteristic image D.

Step 2.5: the obtained binarization characteristic D is exclusive-or with the binarization copyright information M to obtain zero watermark information WM, origin coordinates in DEM data, grid row and column numbers and operation time are converted into keys, the keys are recorded as K, and WM, coefficient matrix H and keys K are registered in a third party organization together.

In some embodiments, the method further comprises a step three, zero watermark extraction. The specific steps of zero watermark extraction are as follows:

step 1: reading DEM data A' to be detected, and obtaining information WM registered in a third party, a coefficient matrix H and a secret key K;

step 2: and judging whether the DEM data A' is cut or not by reading the original point coordinates and the grid number in the secret key. If the cutting is performed, the step 3 is performed, and if the cutting is not performed, the step 4 is performed;

step 3: according to the key information, 0 is complemented to the clipping part, then the data A ' is segmented according to 16 x 16, converted into column vectors and extracted into a matrix I, the corresponding column vectors are obtained from H according to the I, H ' is formed, and the data A ' is divided into blocks according to the formula W=I [ H ' '] ^-1 Calculating a base matrix W, multiplying W by H stored in a third party, and obtaining reconstruction data B';

step 4: extracting invariant moment characteristics of the data A ' to be detected or the reconstruction data B ', and obtaining a characteristic image D ' in the same step II;

step 5: and D ' is exclusive-or ' with WM, arnold scrambling is carried out on the D ' to obtain M ', and the M is compared with the M ', so that copyright authentication is carried out.

In order to verify that the zero watermark method after NMF decomposition processing has better robustness in the face of large-space clipping attack. In this embodiment, the robustness and uniqueness experiments are performed by selecting data respectively, and because the researches on DEM zero watermarks are less, the present gray image algorithm is selected and compared in this embodiment, so as to embody the advantage of the method in the strong clipping resistant direction.

Robustness experiments:

in the embodiment, DEM data in regions of 31 degrees of north latitude and 103 degrees of east longitude are selected as original carrier data, the data name is ASTMTMTY 003_N31E103.GIF, a two-dimensional display diagram is shown in (a) of fig. 3, and the data format is tif; meanwhile, the present section uses a meaningful binary image as copyright information, and its size is 64×64, as shown in fig. 3 (b).

In order to verify the robustness of the present method, the present embodiment designs an attack experiment, constructs the completed DEM data by attacking the zero watermark, extracts the copyright information after the attack, compares the copyright image with the extracted copyright information, and calculates the NC value to determine the robustness of the experiment (the NC value results each retain two decimal places).

(1) Strong cutting attack resistance experiment

Various regular and irregular clipping attacks are mainly enumerated, and robustness is verified through experiments. Fig. 4 shows the clipping effect of clipping with different intensities, the recovery effect of NMF transformation and the watermark extraction effect after facing such an attack, and NC values obtained by reading the zero watermark copyright image and the original copyright image under the corresponding clipping attack are shown in table 1.

TABLE 1 NC values for extracting copyright information under different attack intensities

As can be seen from fig. 4, comparing the original data with the inverse NMF-based generated data after different intensity cropping attacks at different locations facing the row, column, edge, upper left corner, etc., no distinction can be seen by the naked eye. And by extracting the zero watermark from the DEM data generated by inverse transformation, the extracted zero watermark information can successfully distinguish logo information of the image, and the logo information cannot be distinguished from the original copyrighted image by naked eyes.

From table 1, it can be known that the NC value calculated by the original copyright information and the zero watermark information extracted by the data seed generated by the inverse NMF transform under different attack intensities is also close to 1 because NMF can recover the data before decomposition almost without loss. Therefore, the watermark algorithm has better robustness in the face of different positions and different strengths of cutting attacks, and also meets the requirement of better robustness of zero watermark in the face of large-area strong cutting.

(2) Translation attack resistance experiment

In the three-dimensional grid zero-watermarking method, the original point coordinates of the header file are used as keys to determine the position, and the watermark is damaged when the data is translated. And (3) performing translation attack on the data from an x axis or a y axis according to the experimental plan, extracting copyright information of the attacked data at a change interval of 100m, and calculating an NC value to judge the anti-attack effect. The attack level, zero watermark extraction effect and NC values are shown in fig. 5.

As can be seen from fig. 5, the NC value of the algorithm is 1.00 regardless of the translation of the data, because the algorithm does not need to use the data in the header file as an index when it is determined that no clipping is performed or clipping is resumed, and therefore the grid width and elevation values remain unchanged when the coordinates change, i.e., the gradient values of each grid remain unchanged and the gradient relationship between grids remains unchanged, so it can be concluded that the method is completely resistant to translation attacks.

(3) Anti-rotation attack experiment

The design experiment tests the robustness of the algorithm in the present section against rotation attack, and the specific design is as follows: and (3) taking the origin in the DEM data header file as a center point, rotating the DEM data with the zero watermark constructed clockwise by 30 degrees, 60 degrees and 90 degrees. After attack, the zero watermark copyright information is extracted, and NC value is calculated, as shown in figure 6.

As can be seen from fig. 6, the NC value calculated by the method herein is always 1.00 when the rotation attack is faced with different angles, because the HU invariant moment used herein has rotation invariance and is binarized by comparing the value of the block HU moment with the value of the mean HU moment, although the relationship between the block and the whole HU moment value is unchanged when the rotation attack is faced with large angle transformation, and experiments prove that the method can resist better when the rotation attack is faced with.

(4) Anti-scaling attack experiment

The design experiment tests the robustness of the method for shortening and prolonging the grid length during attack, extracts copyright information in the attacked data, and specifically designs the method as follows, the grid width is respectively scaled to be 0.5 times of the original width, and is respectively scaled to be 2 times and 5 times of the original width. The result of comparing the post-attack extraction effect with the calculated NC value is shown in fig. 7.

As can be seen from fig. 7, the NC value extracted after the scaling attack is always kept at 1.00, because the mesh is only used as a positioning index when the algorithm is designed to resist the cropping attack, and the width of the mesh is not used as a feature element participating in calculation when other attacks are resisted. Therefore, when the scaling attack is faced, even though the gradients are changed, the magnitudes among the gradients are not changed, and the magnitudes among HU invariant moments of the gradients calculated correspondingly are unchanged, so that the scaling cannot influence zero watermark information.

(5) Experiment for resisting elevation translation attack

The design experiment tests the robustness of the method facing the elevation translation attack, and extracts zero watermark information from the data after the attack, and the specific design is as follows, the elevation value is respectively increased by 50m, 100m and reduced by 50m. The result of comparing the post-attack extraction effect with the calculated NC value is shown in fig. 8.

As can be seen from fig. 8, after the scaling attack, the NC value obtained by calculation is always 1.00, and similar to the translational attack, the overall deviation of the elevation occurs, but the corresponding gradient is transformed, so that the calculated HU moment and the compared binary are not transformed, and the method can be completely resisted when facing the elevation translational attack.

(6) Anti-noise attack experiment

Experimental tests were designed to test the robustness of the methods herein in the face of noise attacks. Noise attacks are commonly known as pretzel noise and gaussian noise, which are used herein to attack DEM data, respectively. When Gaussian noise attack is designed, the mean value is set to be 0, and the standard deviation is set to be 0.1; when designing salt and pepper noise, due to the difference between DEM data and gray images, 0 elevation value and the maximum value in the data are designed to be taken into experiments as white noise and black noise. As can be seen from fig. 9, the calculated NC value gradually decreases as the gaussian noise attack intensity gradually increases. It can be seen that by experiment, the noise in time reaches 0.05, the nc value is still within the threshold value of 0.75, and the same conclusion can be drawn from fig. 10, so that the method can resist noise attack within a certain intensity.

Uniqueness analysis:

the design experiment verifies whether the Wen Ling watermarking method meets the uniqueness requirement when being applied to DEM data copyright protection. The carrier data in the robustness experiment is selected as the original data 1 to be constructed, the comparison data 2 is common DEM test data SRTM3_V4_90m.GIF, and the two comparison data are rough, as shown in (a) and (b) in fig. 11, but the region where the data are located and the terrain exhibited by the data are different from those of the original data.

And respectively constructing zero watermarks for the two groups of data, calculating NC values, and comparing the NC values of the two groups of data to judge whether the uniqueness of the zero watermarks meets the standard or not, wherein the comparison result is shown in figure 12. Since the zero watermark is constructed from the original data, the NC value of the copyright information extracted from the original data is necessarily 1.00. As can be seen from the table, the NC value of the comparison data 2 is 0.53, which is far smaller than the empirical threshold value of 0.75, and the extracted copyright information presents messy codes and cannot be recognized. Therefore, the method has better uniqueness and meets the requirements of DEM copyright protection.

Comparison experiment:

in order to further embody the superiority of the algorithm in the aspect of resisting strong clipping attack, taking an original picture as an example, because the research on the grid DEM data is less, a zero watermark algorithm on gray image data is selected, namely, compared with the documents [19] and [20], the document 19] adopts NSST to decompose a host image; then randomly extracting a sub-graph from the decomposed low-frequency approximation components by using a Logistic chaotic system, and dividing the sub-graph into mutually non-overlapping sub-blocks; and finally, performing QR decomposition on each sub-block, and constructing a zero watermark by judging the size relation between the l1 norm of the 1 st row element vector in each sub-block R matrix and the average value of the l1 norms of the 1 st row element vectors of all sub-block R matrices. Document [20] after normalization processing, firstly obtaining low-frequency components of an image through NSST and DCT transformation, and then adopting NMF to decompose and extract feature vectors as feature information for constructing a zero watermark; when processing the copyright image, a spread spectrum technology and a chaotic scrambling technology are adopted to improve the safety of the algorithm. As particularly shown in fig. 13 and 14.

When facing a large-area clipping attack, the larger the NC value calculated by the zero watermark is, the smaller the BER value is, and the stronger the clipping attack resistance of the corresponding algorithm is. As can be seen from fig. 13 and 14, when the cutting attack is performed on a small area, all three can be well resisted, and the robustness is very high, but as the cutting area increases, NC values of the document [19] and the document [20] gradually decrease, BER values gradually increase, the ability to resist the cutting attack becomes weak, and the robustness is reduced. From this, it can be derived that the algorithm herein is very robust and has better stability compared to literature [19] and literature [20 ].

Document [19] is [19] Han Shaocheng, zhang Zhaoning, zhang Yujin ] a robust zero-watermark algorithm based on non-downsampled shear wave transformation and QR decomposition [ J ]. Opto-electronic laser, 2012,23 (10): 1957-1964.

Document [20] refers to [20]Sun S,Zeng Y,Zhang H,et al.A Robust Image Perceptual Hashing Scheme Based on Improved Normalization Algorithm[J ]. Advanced Materials Research,2012,1864 (532-533).

In summary, the application designs a zero watermark method aiming at strong clipping attack by analyzing the problems of the digital watermark technology in the DEM data application process. The NMF decomposition technology is utilized to reduce the cutting part through the residual part, so that the effect of strong cutting resistance is achieved. And the purpose of resisting the attacks such as rotation and the like is achieved by calculating the invariant moment of the gradient HU. Experiments show that the method has strong robustness for large-area clipping attacks and also has good robustness for other attacks. The uniqueness experiment shows that the method can effectively identify copyright attribution when facing different data, and has a certain use prospect and popularization value.

Furthermore, although exemplary embodiments have been described herein, the scope thereof includes any and all embodiments having equivalent elements, modifications, omissions, combinations (e.g., of the various embodiments across), adaptations or alterations as pertains to the present application. The elements in the claims are to be construed broadly based on the language employed in the claims and are not limited to examples described in the present specification or during the practice of the application, which examples are to be construed as non-exclusive. It is intended, therefore, that the specification and examples be considered as exemplary only, with a true scope and spirit being indicated by the following claims and their full scope of equivalents.

The above description is intended to be illustrative and not restrictive. For example, the above-described examples (or one or more aspects thereof) may be used in combination with each other. For example, other embodiments may be used by those of ordinary skill in the art upon reading the above description. In addition, in the above detailed description, various features may be grouped together to streamline the application. This is not to be interpreted as an intention that the features of the claimed application are essential to any of the claims. Rather, inventive subject matter may lie in less than all features of a particular inventive embodiment. Thus, the following claims are hereby incorporated into the detailed description as examples or embodiments, with each claim standing on its own as a separate embodiment, and it is contemplated that these embodiments may be combined with one another in various combinations or permutations. The scope of the application should be determined with reference to the appended claims, along with the full scope of equivalents to which such claims are entitled.

Claims (10)

1. An NMF decomposition-based DEM zero watermarking method, the method comprising:

acquiring a copyright image, and scrambling the copyright image to obtain information M of the encrypted copyright image;

extracting an elevation matrix from original DEM data, performing non-overlapping blocking on the elevation matrix, converting each sub-block into a column vector to form a new matrix B, and performing NMF decomposition on the matrix B to obtain a base matrix W and a coefficient matrix H;

calculating the gradient of each grid in the original DEM data to obtain a gradient image C;

non-overlapping blocking is carried out on the gradient image C, at least three HU invariant moments of each small block and the gradient image C are calculated respectively, and the three high-order moment moments are averaged;

comparing the HU invariant moment average value of the gradient image C with the HU invariant moment average value of each small block, and recording as 1 when the whole is larger and as 0 when the whole is smaller to obtain a binary characteristic image D;

and (3) carrying out exclusive OR on the binary characteristic image D and the encrypted copyright image information M to obtain zero watermark information WM, and converting the original point coordinates, the grid row number and the operation time in the original DEM data into a secret key K.

2. The method of claim 1, wherein the copyrighted image is scrambled by the formula:

where mod is used to extract the modulus in the formula, N is the row width of square copyright data, (x ', y') is the processed coordinates, and (x, y) is the coordinates in the original information.

3. The method according to claim 1, wherein the performing NMF decomposition on the matrix B to obtain the base matrix W and the coefficient matrix H includes:

based on the matrix B, randomly initializing to respectively obtain a base matrix W and a coefficient matrix H, wherein W is more than 0 and H is more than or equal to 0.

The base matrix W is iteratively updated according to the update rule shown in the above equation (5):

wherein i represents the row number of the original matrix, j represents the column number of the original matrix, and B is approximately equal to W and H after multiple iterations;

normalizing the rows and columns of the base matrix W;

iteratively updating the coefficient matrix H through the base matrix W according to an updating rule shown in a formula (5);

and comparing the calculated error with a first threshold value, if the calculated error is smaller than or equal to the first threshold value, finishing operation, outputting an operation result, and if the calculated error is larger than the threshold value, iteratively updating the base matrix W according to an updating rule shown in the formula (5) until the calculated error is smaller than or equal to the first threshold value.

4. A method according to claim 3, characterized in that the update law is determined by:

taking the Euclidean distance square as an objective function, extracting a matrix B of original DEM data, calculating a reconstructed matrix WH, and judging whether iteration is continued or not by calculating a distance value between the matrix B and the matrix WH and comparing the distance value with a second threshold value, wherein a Euclidean distance square formula is shown in a formula (2):

D(W,H)＝||B-WH|| ² ＝∑ _ij (B _ij -(WH) _ij ) ² (2)

and respectively deriving W and H on two sides of the formula (2), wherein the derived result is shown as the formula (3):

wherein T is a rank-converting operator, and k is a decomposition dimension;

by means of the gradient descent method, an update rule is deduced as shown in formula (4):

order theThe update method shown in the formula (4) is simplified from addition to multiplication, and as shown in the formula (5), the update rules of the base matrix W and the coefficient matrix H are obtained respectively.

5. The method of claim 1, further comprising zero watermark extraction.

6. The method of claim 5, wherein the zero watermark extraction comprises:

reading DEM data A' to be detected, and obtaining zero watermark information WM, a coefficient matrix H and a secret key K;

judging whether the DEM data A' is cut or not according to the original point coordinates and the grid number in the secret key; if the cutting is performed, the step 3 is performed, and if the cutting is not performed, the step 4 is performed;

under the condition that DEM data A ' is cut, 0 is complemented to a cutting part according to key information, then the data A ' is converted into column vectors to be extracted into a matrix I, corresponding column vectors are obtained from H according to the matrix I, H ' is formed, and W=I [ H ' '] ^-1 Calculating a base matrix W, multiplying W by a coefficient matrix H, and obtaining reconstruction data B';

under the condition that DEM data A 'is not cut, performing invariant moment feature extraction on the data A' to be detected or the reconstruction data B 'to obtain a feature image D';

and exclusive-or the characteristic image D ' and the zero watermark information WM, and performing Arnold scrambling on the characteristic image D ' to obtain the extracted information M ' of the encrypted copyright image.

7. The method according to claim 6, wherein the copyright authentication is performed by comparing the information M of the encrypted copyright image with the information M' of the encrypted copyright image.

8. The method of claim 6, wherein the clipping portion is complemented with 0 according to the key information, and the data a 'is converted into a column vector to be extracted into a matrix I, and the corresponding column vector is obtained from H according to I, and H' is formed by the formula w=i [ H ''] ^-1 Calculating a base matrix W, multiplying W by a coefficient matrix H, and obtaining reconstruction data B', wherein the method comprises the following steps:

the basic model of NMF decomposition is shown as formula (6):

B _m*n ＝W _m*r *H _r*n (6)

wherein m, n are the dimensions of the non-negative matrix B, and r is the dimension of NMF decomposition;

b _i is the element of the ith column in the non-negative matrix B, h _i For the element in the ith column in the coefficient matrix H, there are:

b _i ＝W _m*r *h _i i∈[1,n] (7)

when the value of the sub-block to be cut is set to 0 after the cutting attack, there are:

B _m*n ＝[b ₁ ,…b _s ,0,…,0,b _t ,…,b _n ] (8)

namely:

B′ _k ＝[b ₁ ,…b _s ,b _t ,…,b _n ] (9)

wherein k represents the number of column vectors and the number of the remaining complete sub-blocks of the original data matrix after the clipping attack, when k is more than or equal to r, namely the number of the matrix sub-blocks after the clipping attack is more than or equal to a dimension value r, the complete base matrix can be recovered, and the inverse operation reconstruction of NMF decomposition is realized;

in the case of B ' being known, both sides are simultaneously multiplied by [ H ', based on B ' =w×h ' '] ^-1 The method comprises the following steps:

W＝B′*[H′] ^-1 (10)

in the case of the known coefficient matrix H, the original matrix is reconstructed according to the above procedure.

9. An NMF decomposition-based DEM zero watermarking method apparatus, characterized in that the apparatus comprises a processor configured to:

acquiring a copyright image, and scrambling the copyright image to obtain information M of the encrypted copyright image;

extracting an elevation matrix from original DEM data, performing non-overlapping blocking on the elevation matrix, converting each sub-block into a column vector to form a new matrix B, and performing NMF decomposition on the matrix B to obtain a base matrix W and a coefficient matrix H;

calculating the gradient of each grid in the original DEM data to obtain a gradient image C;

non-overlapping blocking is carried out on the gradient image C, at least three HU invariant moments of each small block and the gradient image C are calculated respectively, and the three high-order moment moments are averaged;

comparing the HU invariant moment average value of the gradient image C with the HU invariant moment average value of each small block, and recording as 1 when the whole is larger and as 0 when the whole is smaller to obtain a binary characteristic image D;

and (3) carrying out exclusive OR on the binary characteristic image D and the encrypted copyright image information M to obtain zero watermark information WM, and converting the original point coordinates, the grid row number and the operation time in the original DEM data into a secret key K.

10. A readable storage medium storing one or more programs executable by one or more processors to implement the method of any of claims 1-7.