Privacy position data distribution estimation method based on geographically undistinguishable model

Abstract

The invention discloses a privacy location data distribution estimation method based on a geographically indistinguishable model. The method comprises a brand-new location data disturbance method and a method for carrying out privacy location data distribution estimation at a server side. According to the method, for multiple users, each user holds one piece of position data, under a fixed privacy budget, the position data of each user is disturbed into a set composed of k random positions, and then the set is submitted to a server; the server divides an interested area into square grids with the same size, and then estimates the number of real positions in each grid according to the received disturbance points and the disturbance probability between the grids. According to the method, the whole protection process can be performed at the user side, and the attack of an untrusted third-party server can be resisted; according to the user side disturbance method and the server side distribution estimation method provided by the invention, the availability of a distribution estimation result is effectively improved.

Description

Technical Field

[0001] The present invention belongs to the field of information security technology and relates to a method for estimating the distribution of private location data based on a geographically indistinguishable model. Background Art

[0002] In recent years, with the rapid development of mobile networks, applications based on location data have penetrated into people's daily lives. Many applications need to collect users' location to provide and improve their services. For example, Google Maps allows users to view the busy time histogram of certain places (restaurants, gas stations, etc.) every day of the week, and can see the busyness of the place in each time period; navigation software can display the congestion of each road (indicated by green, yellow, red, dark red) and estimate the time required to pass this road. These services require a large amount of collection of location data of many users, but the location data of users may be sensitive (such as hospitals, workplaces, and home addresses). Moreover, the servers that collect this data may not be trustworthy. If this data is abused, it may cause serious privacy issues. Therefore, how to ensure the service quality of applications while protecting the privacy of users' location is a very meaningful issue.

[0003] In order to reduce the risk of location privacy leakage, Bamba proposed a location anonymity model in 2008, which anonymizes location data through location generalization. However, k-anonymity technology can only provide weak privacy protection, for example, it is vulnerable to background knowledge attacks. Differential privacy (DP) was proposed by Dwork in 2006 as a strict privacy model that can resist background knowledge attacks, and has been widely studied to solve the problem of location privacy risks. It can ensure that adding or deleting any location data in the database can produce indistinguishable results. With the development of technology, there is a problem of untrustworthy third parties in real life, and local differential privacy (LDP) came into being. This model inherits the advantages of differential privacy, and all users perturb the data locally before uploading it to the server, and will not expose the real data to any third party. The location data protection technology under the local differential privacy model has also been studied in depth. As a variant of local differential privacy, Andres proposed the Geo-Indistinguishability (GI) model in 2013. Compared with the traditional local differential privacy model, this model allows two geographical locations with similar Euclidean distances to have a higher probability of generating the same reported location. It allows two locations to be indistinguishable within a radius of r.

[0004] To date, many geo-indistinguishability mechanisms have been proposed to protect users' location privacy in various applications. However, in these mechanisms, the server does not post-process the perturbed data, and the quality of services related to data distribution estimation will be affected to a certain extent. Summary of the invention

[0005] Purpose of the invention: The present invention aims to provide a method for estimating the distribution of private location data based on a geographically indistinguishable model, which effectively improves the availability of distribution estimation results over the user-side perturbation method and the server-side distribution estimation method.

[0006] In order to achieve the above-mentioned object of the invention, the technical solution provided by the present invention is as follows.

[0007] A method for estimating the distribution of private location data based on a geographically indistinguishable model comprises the following steps:

[0008] S1. User location data disturbance: Each user has a geographical location data represented by two-dimensional rectangular coordinates. The original data of all users is represented by the set (x, y) = {(x1, y1), (x2, y2), …, (x n ,y n )}, each user independently according to the probability density function Select a random number r, and then select k points uniformly and randomly from a circle with the real position as the center and r as the radius and submit them to the server;

[0009] S2. Grid division: Before collecting user data, the server first calculates the total noise error and non-uniform error according to the set privacy budget ε, set size k, and the number of participating users n. Then, the plane area of ​​interest is divided into m square grids of equal size and side length 2R according to the uniform grid division method. The Hadamard matrix H of order sets a candidate set for each grid;

[0010] S3. Calculation of disturbance probability matrix: When the grid is dense enough, the coordinates of the center point of the grid are used to replace all other positions in the grid to calculate the disturbance probability. The grid is approximately regarded as a circle with a radius of R. The server calculates the disturbance probability between every two grids and combines all disturbance probabilities into an m-order disturbance probability matrix P.

[0011] S4, data distribution estimation: When the server receives the perturbation data submitted by all users, it first counts the frequency of the perturbation points in each grid, multiplies it with the Hadamard matrix to get the frequency of the perturbation points in each candidate set, represents the frequency of all candidate sets with vector p(C), and then solves the linear equation system Get an estimate of the true data frequency within all grids.

[0012] In step (S1), the user end perturbs the original data including the following process:

[0013] Assume that each user holds a geographic location data represented by a two-dimensional rectangular coordinate system, and user i holds private location data represented by (x i ,y i ), the original data of all users can be represented by the set (x,y)={(x1,y1),(x2,y2),…,(x n ,y n )} means that all users independently follow the probability density function Sampling random number r, random number r satisfies gamma distribution, and then independently extracting k random numbers θ1, θ2, …, θ from a uniform distribution from 0 to 2π with the true position as the center k , and then independently draw k random numbers l1, l2, ..., l from a uniform distribution between 0 and 1 k , all the random numbers required by user i in the data perturbation are sampled, and then the set {(x i1 ,y i1 ),(x i2 ,y i2 ),…,(x ik ,y ik )} is submitted to the server, that is, k points are uniformly randomly selected within a circle with a radius of r and submitted to the server, where the data perturbation is as follows:

[0014]

[0015]

[0016] In the formula, r i Represents the random number sampled by user i according to the probability density P(r).

[0017] The perturbation data submitted by each user is a set, k represents the size of the set, and all users use the same k.

[0018] Furthermore, in step (S2), the server divides the plane area into grids and confirms the candidate set of each grid, including the following process:

[0019] Assume that the server divides the plane area of ​​interest into m square grids of equal size and side length 2R, and then calculates the total noise error and non-uniform error according to the set privacy budget ε, set size k, and number of participating users n using the uniform grid division (UG) method. The specific calculation process is as follows:

[0020] The standard deviation between the actual number of locations and the estimated number of locations in each grid is in is a quantity that depends mainly on the value of ε, p(C x )and They refer to the actual perturbation probability of the candidate set of grid x and our estimated perturbation probability, respectively. It can be regarded as a function of m b = b (m), p x→y refers to the perturbation probability from grid x to grid y; because the value of a mainly depends on ε, we can omit it here. Assuming r is the ratio of the query area to the total area size, the total noise error is

[0021] The non-uniformity error is proportional to the number of data points in the grid that are located at the query boundary. Considering the query area as a square, the proportion of each edge is Then on the query boundary there is grids, assuming that N = nk is the number of all disturbance points, then the number of points on the boundary can be considered to be

[0022] The next task is to find the m that minimizes the sum of the two types of errors. The sum of the errors can be considered as a function of m. The optimal m can be found using the following optimization method:

[0023]

[0024] use The i+1th row of the Hadamard matrix H defines the candidate set of the i-th grid. The index of the column corresponding to all the 1 elements in the i+1th row of the matrix H (the index value is less than or equal to m) represents the candidate set C of the i-th grid. i .

[0025] For ease of understanding, take m=3 as an example. The Hadamard matrix of order is Then the candidate sets of grids 1, 2, and 3 are C1 = {1,3}, C2 = {1,2}, and C3 = {1} respectively.

[0026] It is further pointed out that in step (S3), the server side calculates the disturbance probability between all grids including the following process:

[0027] (S31) When the grid is dense enough, the coordinates of the center point of the grid can be used to replace all other positions to calculate the disturbance probability. The grid is approximately regarded as a circle with a radius of R. The probability of data in each grid being disturbed to the same grid or any other grid can be approximately calculated, and all disturbance probabilities are combined into an m-order disturbance probability matrix P. The specific calculation method of the disturbance probability is as follows:

[0028] (S32) The data before and after the disturbance are not in the same grid: Assume that the Euclidean distance between the center points of the two grids before and after the disturbance is d. When the random radius size r∈[0,dR], the random circle has no overlap with the other grid, so the disturbance probability is 0; when r∈(dR,d+R), the random circle intersects with the other grid, so the disturbance probability is Where S represents the area of ​​the intersection of two circles whose center distance is d and whose radius is R respectively. When r∈[d+R,+∞), the random circle includes another grid, so the perturbation probability is So the probability of the true position being disturbed from one grid to another can be approximately expressed as:

[0029] The data before and after the disturbance are in the same grid: When r∈[0,R], the disturbance point must be in the original grid, so the disturbance probability is When r∈(R,+∞), the perturbation probability is So the probability that the true position and the perturbation position are in the same grid is

[0030] In step (S4), after receiving the disturbance data, the server estimates the distribution of the original data, including the following process:

[0031] When the server receives the perturbation data submitted by all users, it first counts the frequency of the perturbation points in each grid, and multiplies it with the Hadamard matrix to get the frequency of the perturbation points in each candidate set. The frequencies of all candidate sets are represented by vector p(C) to solve the linear equation system Get the estimated value of the real data frequency in all grids; the specific process is as follows:

[0032] The perturbation data submitted by n users are represented by the set {{z 11 ,z 12 ,…,z 1k},{z 21 ,z 22 ,…,z 2k},…,{z n1 ,z n2 ,…,z nk}} means that after receiving the disturbance data, the server first counts the number of disturbance points in each grid. The fast Walsh-Hadamard transform can efficiently find the number of disturbance points in the candidate set of all grids, and the set {F1,F2,…,F m}, then the frequency of the disturbance point in the candidate set of the grid x can be expressed as p(C)={p(C1),p(C2),…,p(C m )}; Assume that the frequency of the real point in the grid x is p x , the perturbation probability from grid x to grid y is p x→y ,According to the approximate value of the inter-grid disturbance probability we obtained, it can be considered that According to this relationship, after we change all -1 in the Hadamard matrix H to 0, the linear equation system The solution p is the estimate of the frequency of the true data in all grids.

[0033] Beneficial effects: The method for estimating the distribution of private location data based on a geographically indistinguishable model described in the present invention first perturbs a user's real data into a set of a certain size by randomly sampling within a circle with the user's real location as the center and a random number as the radius, and submits the set to a server, thereby protecting the user's privacy; on the server side, the present invention proposes a method for estimating the distribution of geographic location data, firstly divides a plane area into square grids of equal size and then calculates a perturbation probability matrix, then uses the properties of the Hadamard matrix to achieve efficient counting, and then constructs a linear equation system to solve the estimated values ​​of the frequency of real data in all grids; thanks to the new perturbation method and post-processing method proposed in the present invention, the effectiveness and accuracy of various services based on geographic location data and the utility of the overall data are improved, while protecting the privacy of users. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 It is a working schematic diagram of the user end and the server end in the present invention;

[0035] Figure 2 It is the overall flow chart of the present invention. DETAILED DESCRIPTION

[0036] The above scheme is further described below in conjunction with specific embodiments. It should be understood that these embodiments are used to illustrate the present invention and are not limited to the scope of the present invention. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the protection scope of the present invention.

[0037] First, for the problem of estimating the distribution of geographic location data under privacy protection, the present invention provides a method for estimating the distribution of privacy location data based on a geographically indistinguishable model. This method is based on a new location perturbation algorithm, which perturbs the single location data of each user into a set of fixed size and then uploads it to the server. Specifically, it refers to a method of random sampling within a circle with the user's real location as the center and a radius of a random number to perturb a user's real data into a set of fixed size. A geographically indistinguishable mechanism; a new location data distribution estimation method is designed on the server side. First, the plane area is divided into square grids of equal size and the perturbation probability matrix is ​​obtained. Then, efficient counting is achieved by relying on the properties of the Hadamard matrix, and then a linear equation system is constructed to solve the estimated value of the frequency of real data in all grids. Efficient counting is achieved by relying on the properties of the Hadamard matrix, and then a linear equation system is constructed. The solution of the equation system is the result we need.

[0038] Specifically, the steps of implementing the method for estimating the distribution of private location data of the present invention are as follows:

[0039] Step 1: User location data disturbance

[0040] Each user holds a geographic location data represented by a two-dimensional rectangular coordinate. Each user independently selects a random number r, and then uniformly randomly selects k points in a circle with its own real coordinate as the center and r as the radius and submits them to the server.

[0041] Assuming that each user holds a geographical location data represented by a two-dimensional rectangular coordinate system, the original data of all users can be represented by the set (x, y) = {(x1, y1), (x2, y2), …, (x n ,y n )} means that all users independently follow the probability density function Select a random number r. It can be seen that the random number r satisfies the gamma distribution, so this sampling process is easy to implement; then select k points uniformly and randomly from the circle with the real location as the center and r as the radius and submit them to the server. The specific perturbation process is as follows: For example, user i holds private location data (x i ,y i ), first sample a random number r according to the probability density P(r) i , and then independently draw k random numbers θ1, θ2, …, θ from a uniform distribution between 0 and 2π k , and then independently draw k random numbers l1, l2, ..., l from a uniform distribution between 0 and 1 k , all the random numbers required by user i in the data perturbation are sampled, and then the set {(x i1 ,y i1 ),(x i2 ,yi2 ),…,(x ik ,y ik )} is submitted to the server, where:

[0042]

[0043] In the formula, r i Represents the random number sampled by user i according to the probability density P(r).

[0044] Step 2: Mesh division

[0045] Before collecting user data, the server first divides the plane area of ​​interest into m square grids of equal size and side length 2R according to the set privacy budget ε, set size k, and number of participating users n. The Hadamard matrix H sets a candidate set for each grid.

[0046] First, assume that the server divides the plane area of ​​interest into m square grids of equal size and side length 2R, and then calculates the total noise error and non-uniform error according to the set privacy budget ε, set size k, and number of participating users n using the uniform grid division (UG) method. The specific process is as follows:

[0047] The standard deviation between the actual number of locations and the estimated number of locations in each grid is in is a quantity that depends mainly on the value of ε, p(C x )and They refer to the actual perturbation probability of the candidate set of grid x and our estimated perturbation probability, respectively. It can be regarded as a function of m b = b (m), p x→y refers to the perturbation probability from grid x to grid y; because the value of a mainly depends on ε, we can omit it here. Assuming r is the ratio of the query area to the total area size, the total noise error is

[0048] The non-uniformity error is proportional to the number of data points in the grid that are located at the query boundary. Considering the query area as a square, the proportion of each edge is Then on the query boundary there is grids, assuming that N = nk is the number of all disturbance points, then the number of points on the boundary can be considered to be

[0049] The next task is to find the m that minimizes the sum of the two types of errors. The sum of the errors can be considered as a function of m. The optimal m can be found using the following optimization method:

[0050]

[0051] use The i+1th row of the Hadamard matrix H defines the candidate set of the i-th grid. The index of the column corresponding to all the 1 elements in the i+1th row of the matrix H (the index value is less than or equal to m) represents the candidate set C of the i-th grid. i , taking m=3 as an example, The Hadamard matrix of order is Then the candidate sets of grids 1, 2, and 3 are C1 = {1,3}, C2 = {1,2}, and C3 = {1} respectively.

[0052] Step 3: Calculation of perturbation probability matrix

[0053] When the grid is dense enough, the coordinates of the grid center point can be used to replace all other positions in the grid to calculate the disturbance probability. The server calculates the disturbance probability between every two grids (including between itself) and combines all disturbance probabilities into an m-order disturbance probability matrix P.

[0054] When the grid is dense enough, the coordinates of the center point of the grid can be used to replace all other positions to calculate the disturbance probability. The grid is approximately regarded as a circle with a radius of R. The probability of data in each grid being disturbed to the same grid or any other grid can be approximately calculated, and all disturbance probabilities are combined into an m-order disturbance probability matrix P. The specific calculation method of the disturbance probability is as follows:

[0055] The data before and after the disturbance are not in the same grid: Assume that the Euclidean distance between the center points of the two grids before and after the disturbance is d. When the random radius size r∈[0,dR], the random circle has no overlap with the other grid, so the disturbance probability is 0; when r∈(dR,d+R), the random circle intersects with the other grid, so the disturbance probability is Where S represents the area of ​​the intersection of two circles whose center distance is d and whose radius is R respectively. When r∈[d+R,+∞), the random circle includes another grid, so the perturbation probability is So the probability of the true position being disturbed from one grid to another can be approximately expressed as:

[0056] The data before and after the disturbance are in the same grid: When r∈[0,R], the disturbance point must be in the original grid, so the disturbance probability is When r∈(R,+∞), the perturbation probability is So the probability that the true position and the perturbation position are in the same grid is

[0057] Step 4: Data distribution estimation

[0058] When the server receives the perturbation data submitted by all users, it first counts the frequency of the perturbation points in each grid, and multiplies it with the Hadamard matrix to get the frequency of the perturbation points in each candidate set. The frequencies of all candidate sets are represented by vector p(C) to solve the linear equation system Get an estimate of the true data frequency within all grids.

[0059] When the server receives the perturbation data submitted by all users, it first counts the frequency of the perturbation points in each grid, and multiplies it with the Hadamard matrix to get the frequency of the perturbation points in each candidate set. The frequencies of all candidate sets are represented by vector p(C) to solve the linear equation system Get the estimated value of the real data frequency in all grids; the specific process is as follows:

[0060] The perturbation data submitted by n users are represented by the set {{z 11 ,z 12 ,…,z 1k},{z 21 ,z 22 ,…,z 2k},…,{z n1 ,z n2 ,…,z nk}} means that after receiving the disturbance data, the server first counts the number of disturbance points in each grid. The fast Walsh-Hadamard transform can efficiently find the number of disturbance points in the candidate set of all grids, and the set {F1,F2,…,F m}, then the frequency of the disturbance point in the candidate set of the grid x can be expressed as p(C)={p(C1),p(C2),…,p(C m )}; Assume that the frequency of the real point in the grid x is p x , the perturbation probability from grid x to grid y is p x→y ,According to the approximate value of the inter-grid disturbance probability we obtained, it can be considered that According to this relationship, after we change all -1 in the Hadamard matrix H to 0, the linear equation system The solution p is the estimate of the frequency of the true data in all grids.

[0061] The following are the experimental results of the present invention. The data sets used in the experiment are Gowalla, T-Drive and Foursquare. The Gowalla data set records 6,442,890 location check-in data of 196,591 users in the United States from February 2009 to October 2010, and each data corresponds to a user number and a longitude and latitude coordinate; the T-Drive data set records about 110,000,000 trajectory data of 17,000 taxis in Beijing from February 2008 to January 2010; the Foursquare data set records 90,048,627 location check-in data of 2,733,324 users in New York. The present invention randomly extracts 100,000 data from each data set for the experiment, and divides them into 8 privacy levels in the experiment, with privacy budgets of 0.1, 0.3, 0.5, 0.8, 1.0, 1.5, 2.0, and 3.0 respectively. In the experiment, each user independently perturbs his own data and then sends it to the server. The measure of the mechanism utility is the square of the L2 norm. and earth moving distance (EMD). In the experiment, the scheme of the present invention is implemented separately with and without post-processing to prove that the data distribution estimation method proposed by us on the server side is beneficial to improving data validity.

[0062] Table 1 Affected by parameter ε

[0063]

[0064]

[0065] Table 2 EMD affected by parameter ε

[0066]

[0067] It can be seen from Tables 1 and 2 that except for the Gowalla dataset with a privacy budget of 0.1, the post-processing method proposed in the present invention significantly improves the utility of the overall data and makes the data distribution estimation results closer to the actual situation.

[0068] The present invention can protect the user's original location data from being obtained by attackers, resist attackers with arbitrary background knowledge and prevent privacy attacks from untrusted third parties. At the same time, it solves the problem that the existing geographical indistinguishability mechanism cannot post-process the data, thereby improving the utility of privacy data while ensuring the security of privacy data.

[0069] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. For those skilled in the art, the present invention may have various modifications and variations. Any modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present invention shall be included in the protection scope of the present invention.

Claims

1. A method for estimating the distribution of private location data based on a geographically indistinguishable model, characterized in that: Includes steps: S1. User location data disturbance: Each user has a geographical location data represented by two-dimensional rectangular coordinates. The original data of all users is represented by the set (x, y) = {(x1, y1), (x2, y2), …, (x n ,y n )}, each user independently according to the probability density function Select a random number r, and then select k points uniformly and randomly from a circle with the real position as the center and r as the radius and submit them to the server; S2. Grid division: Before collecting user data, the server first calculates the total noise error and non-uniform error according to the set privacy budget ε, set size k, and the number of participating users n. Then, the plane area of ​​interest is divided into m square grids of equal size and side length 2R according to the uniform grid division method. The Hadamard matrix H of order sets a candidate set for each grid; S3. Calculation of disturbance probability matrix: When the grid is dense enough, the coordinates of the center point of the grid are used to replace all other positions in the grid to calculate the disturbance probability. The grid is approximately regarded as a circle with a radius of R. The server calculates the disturbance probability between every two grids and combines all disturbance probabilities into an m-order disturbance probability matrix P. S4, data distribution estimation: When the server receives the perturbation data submitted by all users, it first counts the frequency of the perturbation points in each grid, multiplies it with the Hadamard matrix to get the frequency of the perturbation points in each candidate set, represents the frequency of all candidate sets with vector p(C), and then solves the linear equation system Get an estimate of the true data frequency within all grids.

2. The method for estimating the distribution of private location data based on a geographically indistinguishable model according to claim 1, characterized in that: In step (S1), the user end perturbs the original data including the following process: Assume that user i holds private location data represented by (x i ,y i ), first sample a random number r according to the probability density P(r) i , and then independently draw k random numbers θ1, θ2, …, θ from a uniform distribution between 0 and 2π k , and then independently draw k random numbers l1, l2, ..., l from a uniform distribution between 0 and 1 k , all the random numbers required by user i in the data perturbation are sampled, and then the set {(x i1 ,y i1 ),(x i2 ,y i2 ),…,(x ik ,y ik )} is submitted to the server, where: In the formula, r i Represents the random number sampled by user i according to the probability density P(r).

3. The method for estimating the distribution of private location data based on a geographically indistinguishable model according to claim 1, characterized in that: In step (S2), the server divides the plane area into grids and confirms the candidate set of each grid, including the following process: The standard deviation between the actual number of locations and the estimated number of locations in each grid is in p(C x )and They refer to the actual perturbation probability and the estimated perturbation probability of the candidate set of grid x, respectively. is a function of m b = b(m), p x→y refers to the perturbation probability from grid x to grid y; assuming r is the ratio of the query area to the total area size, the total noise error is The non-uniformity error is proportional to the number of data points in the grid that are located at the query boundary. Considering the query area as a square, the proportion of each edge is Then on the query boundary there is grids, assuming that N = nk is the number of all disturbance points, then the number of points on the boundary is The next task is to find the m that minimizes the sum of the two types of errors, considering the sum of the errors as a function of m The optimal m is obtained by the following optimization method: use The i+1th row of the Hadamard matrix H of order 1 defines the candidate set of the i-th grid, and the index of the column corresponding to all the 1 elements in the i+1th row of the matrix H represents the candidate set C of the i-th grid. i .

4. The method for estimating the distribution of private location data based on a geographically indistinguishable model according to claim 1, characterized in that: In step (S3), the server side calculates the disturbance probability between all grids, including the following two cases: (S31) The data before and after the disturbance are not in the same grid Assuming that the Euclidean distance between the center points of the two grids before and after the disturbance is d, when the random radius size r∈[0,dR], the random circle has no overlap with the other grid, so the disturbance probability is 0; when r∈(dR,d+R), the random circle intersects with the other grid, so the disturbance probability is Where S represents the area of ​​the intersection of two circles whose center distance is d and whose radius is R respectively: When r∈[d+R,+∞), the random circle includes another grid, and the perturbation probability is So the probability of the true position being disturbed from one grid to another is approximately expressed as: (S32) The data before and after the disturbance are in the same grid When r∈[0,R], the disturbance point must be within the original grid, so the disturbance probability is When r∈(R,+∞), the perturbation probability is So the probability that the true position and the perturbation position are in the same grid is:

5. The method for estimating the distribution of private location data based on a geographically indistinguishable model according to claim 1 is characterized in that: in step (S4), after receiving the disturbance data, the server estimates the distribution of the original data, including the following process: The perturbation data submitted by n users are represented by the set {{z 11 ,z 12 ,…,z 1k },{z 21 ,z 22 ,…,z 2k },…,{z n1 ,z n2 ,…,z nk }} means that after receiving the disturbance data, the server first counts the number of disturbance points in each grid. The fast Walsh-Hadamard transform can efficiently find the number of disturbance points in the candidate set of all grids, using the set {F1,F2,…,F m }, then the frequency of the disturbance point in the candidate set of the grid x is expressed as p(C)={p(C1),p(C2),…,p(C m )}; Assume that the frequency of the real points in the grid x is p x , the perturbation probability from grid x to grid y is px →y ,According to the approximate value of the inter-grid disturbance probability, it is believed that According to this relationship, after changing all -1 in the Hadamard matrix H to 0, the linear equation system The solution p is the estimated value of the real data frequency in all grids.

Images